texts: How to create a map to PNG file with specific box coordinates by Mapnik?

texts:
I recently manually built a tile server with this instruction ->
https://switch2osm.org/serving-tiles/manually-building-a-tile-server-20-04-lts/.
I downloaded location data for Russia from geofabrik.de and then loaded to database by osm2pgsql.
Now all working nicely - Apache, mod_tile, renderd, Mapnik and PostgreSQL with PostGIS.
Now I want to create a program on Python for generating a PNG file with specific region (box). This region will be defined by coordinates (left and right corners of box).
I write this code:
import mapnik

mapnik_xml = "openstreetmap-carto/mapnik.xml"
map_output = "region_map_mapnikXml.png"

# Create a map object
m = mapnik.Map(600,300)
mapnik.load_map(m, mapnik_xml)
bbox = mapnik.Box2d(45.23, 41.34, 46.23, 43.23)
m.zoom_to_box(bbox)
mapnik.render_to_file(m, map_output)
print(f"Rendered image to {map_output}")

But after executing this code created only blank file with background.
I don't figure out why?
Can you help me please.
Additionaly, if I run this code:
import mapnik

mapnik_xml = "openstreetmap-carto/mapnik.xml"
map_output = "world_map_mapnikXml.png"

# Create a map object
m = mapnik.Map(600,300)
mapnik.load_map(m, mapnik_xml)
m.zoom_all()
mapnik.render_to_file(m, map_output)
print(f"Rendered image to {map_output}")

in results, generated PNG file - world map with shape of mainlands and countries.
Why is this happening?


texts: Reversing orientation of polygon rings in GDAL/C#?

texts:
I have a shapefile containing only one polygon. I can visualize this file in QGIS without a problem. However when opening this file with C#/GDAL, I cannot create a SqlGeography out of that polygon feature.

I figured out that this is because this polygon has wrong orientation. I wrote code to reverse the geometry's orientation but I feel my inadequate knowledge of GDAL/C# is making me miss a simpler way to achieve orientation reversal.

How do I reverse the orientation of (rings of?) a polygon using C# bindings of GDAL?

Below is the code I wrote to get the orientation reversed. This code doesn't give me the desired result. The line where I add a reversedRing to reversedGeometry throws an exception about adding a ring that is not Linear Geometry.

private static Feature ReverseOrientationForGeometryOf(Feature feature)
{
    //Step - 1: Get the original feature's geometry.
    var geometry = feature.GetGeometryRef();

    //Step - 2: Create a blank geometry.
    var reversedGeometry = new Geometry(geometry.GetGeometryType());
    reversedGeometry.AssignSpatialReference(geometry.GetSpatialReference());

    //Step - 3: Iterate over all the geometries in the original geoemtry.
    for (var geomIndex = 0; geomIndex < geometry.GetGeometryCount(); geomIndex++)
    {
        //The first geometry, i.e. the 0th geometry is outer ring. 1st and onwards are inner rings.

        //Get the ring at this index.
        var ring = geometry.GetGeometryRef(geomIndex);

        //Reverse the ring.
        var reversedRing = Reverse(ring);

        //Add revered ring to the blank geommetry created above in Step - 2.
        reversedGeometry.AddGeometry(reversedRing); //!!THROWS EXCEPTION
    }

    //Step - 4: Now set the reversed geom as the geometry of the input feature.
    feature.SetGeometry(reversedGeometry);
    return feature;
}

private static Geometry Reverse(Geometry ring)
{
    var reversedRing = new Geometry(ring.GetGeometryType());
    reversedRing.AssignSpatialReference(ring.GetSpatialReference());
    var point = new double[3];
    for (var pointIndex = ring.GetPointCount() - 1; pointIndex >= 0; pointIndex--)
    {
        ring.GetPoint(pointIndex, point);
        reversedRing.AddPoint_2D(point[0], point[1]);
    }
    return reversedRing;
}


texts:
I have a model which iterates through feature classes in a geodatabase.
For each feature class, it runs the check Geometry tool.
I then use the "Get Count" tool to count the number of rows in the output table.
I then want to add an expression where if the number of rows is greater than 0 it should keep running the rest of the model, otherwise it should go on to the next feature class.

I thought this could be achieved by using the Calculate Value tool using the following code.

Expression

CountMe(%Row Count%)

Code Block

 def CountMe(n):
     import arcpy
     if n > 0:
      return "true"
     else:
      return "false"

Where %Row Count% is the output from the "Get Count" tool.

This returns True and False values, and I have the output from this set as a precondition to the next steps in the model.

However, I assumed that only "true" data would go through, but all the data seems to go through i.e. the model will push through an empty table from the Check Geometry tool all the way to the end.

I've had a read online and I think that using the Calculate Value tool is the wrong way to go about this. Am I right in assuming that what I need is a python script that will give me 2 outputs (True and False) and then I can connect the "True" as a precondition to the rest of the model, and leave the "False" disconnected, so the model stops and the iterator goes on to the next feature class?

If so, would anybody be able to give me some help with this?


texts: How to extract slope from topographic raster map?

texts:
I am working on a project for a machine learning class wherein I extract features from a topographic raster map and then classify the traversability of the terrain. One of the features I need to extract is slope. My initial thought was that an approximation of slope at any given point could be extracted from a map by calculating the distance between contour lines. However it seems I took for granted the complexity of extracting the contour lines.

In attempting to solve this problem I found myself trying to learn a number of GIS tools. However I still have not been successful which is why I am here. I have installed gdal and used gdal_contour to generate what I believe to be a DEM. From there I loaded it into qgis along with the original raster map. From this I see that when gdal created the DEM it used only one of the two colors that designate contour lines. Worse yet it used the less frequent color so the resulting contour lines are highly segmented.

This leaves me with a couple of questions. The first is whether it is possible to designate the colors gdal uses to extract contour lines and generate the DEM. The second is whether there is a simpler method for extracting the desired data.


texts:
I made this tutorial. I used Landsat 8 bands (1-7), which I atmospheric corrected with DOS and converted to 8-bit. After that I layer stacked the pre processed band (1-7), clipped a little part with a polygon shape of the whole image and used it as input for "Exact Large-Scale Mean Shift segmentation step 1".

But after step 4 I just get a very bad segmented shapefile. The tool just segmented a little percentage of the objects in my Landsat 8 clip. Anyone has a idea how to fix this?


texts:
I need to create an cluster and outlier analysis for each field in my shapefile.

How can I create an iterator that changes the input value in each run?


texts: Exporting Contours to DWG with Elevation Data

texts:
I do not have the 3D Analyst extension. I called one of my old professors and he told me to use features to 3D tool and then convert it to a point file. I need to export the contours to a CAD file that contains the elevation data.

I need a work around that does not require the 3D analyst extension, and also is there any reason why my professor said to convert it to a point file?

They are asking me for contours not points.


texts: Why doesn't Leaflet fill the height of its parent container?

texts:
I have a ` with a specific height, inside of it I want to put my leaflet map. But the problem is that the map doesn't fill 100% of the height of its container. Why does my map not show up?

HTML is as follows:

CSS:

.col-sm-9{
  height: 500px;
}

#map{
  height: 100%;
}

And JS:

var map;
map = L.map( 'map', {
  center: [33.960057, -6.916462],
  minZoom: 2,
  zoom: 14,
  zoomControl: true
});
L.tileLayer('http://{s}.google.com/vt/lyrs=s&x={x}&y={y}&z={z}',{
  maxZoom: 20,
  subdomains:['mt0','mt1','mt2','mt3']
}).addTo( map );


texts: Choosing one certain point (out of many) for nearby polygon based on coordinates and other similar attributes

texts:
I am working on a project in which I have to find a certain lightning strike that caused an area to burn. For now, I have 1300 fire areas with approx. 27,000 lightning strikes in my QGIS project map. Both datasets have coordinates and fire/stroke dates.

For each fire area there 'should be' a lightning strike which sparked a fire (my task). However, as it is right now there are multiple lightning strikes within and nearby each fire area scattered over time (from 2007 to 2017). I have to be able to find the exact lightning strike based on the coordinates and from the dates of fire and lightning strikes.

Additional: I need it for every fire area (1300 areas) so yes automated as manual selection would take way too long. Further, I have included all lightning strikes of a fire area within a buffer of 1 km from the outer parts of the fire area. So in terms of coordinates, we are not talking about changes of more than .1 or .2 The temporal uncertainty can be bigger, as the lightning dataset varies largely within itself. Let's say 3 to 4 weeks of uncertainty. Oh, and output is less important, preferably a shapefile.

Is this possible with QGIS and, if so, how?


texts: How to lock LayoutItem layers in QGIS 3.0 using python?

texts:
I am using the python console to interact with QGIS 3.0 and have accessed the layout manager using:

lomgr = QgsProject.instance().layoutManager()
lo = lomgr.layouts()[0]
myMapItem = lo.itemById('myMapItem')

I need to lock the layers programmatically as I would if I were using the LayoutManager UI (e.g., select the item, go to Item Properties tab, and check/uncheck the 'Lock Layers' checkbox).  I can't find any way of doing this from the documentation (https://python.qgis.org/api/classQgsLayoutItem.html).


texts: Distance between two consecutive points where condition applies for both points

texts:
Is it possible to run through a data.frame and calculate the distance between the points if both point fulfil a condition?

In the enclosed example file i would want to calculate the the distance between two points where values in the actsum column equals or exceeds 0.26.

E.g. the only case where a distance should be calculated is the points in line 6 and 7

Would pointDistance {raster} in be usefull for this?


texts: Opening a shapefile with no extension in QGIS

texts: I am trying to open an open license shapefile in QGIS (this is the link: https://pub.uni-bielefeld.de/record/2674065), I am not sure if the information is not there or if the content is not complete; in the package there is only one file without extension...I am a begginer so probably I am doing something wrong.


texts: Geoserver's mosaic created out of 3 pictures, shows only 1 picture

texts:
I am trying to create mosaic in geoserver out of 3 jpg pictures.Every picture is accompanied by appropriate .jgw and prj file.
IMG_3293.jgw

0.0300000
0.00000000
0.00000000
-0.0300000
19.80492000
45.30164700

IMG_3298.jgw

0.0300000
0.00000000
0.00000000
-0.0300000
19.80416700
45.30286500

IMG_3299.jgw

0.0300000
0.00000000
0.00000000
-0.0300000
19.80431000
45.30321300

And i created .prj files for each picture as:
IMG_3293.prj

PROJCS["WGS 84 / Pseudo-Mercator",
  GEOGCS["WGS 84",
    DATUM["World Geodetic System 1984",
      SPHEROID["WGS 84", 6378137.0, 298.257223563, AUTHORITY["EPSG","7030"]],
      AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]],
    UNIT["degree", 0.017453292519943295],
    AXIS["Geodetic longitude", EAST],
    AXIS["Geodetic latitude", NORTH],
    AUTHORITY["EPSG","4326"]],
  PROJECTION["Popular Visualisation Pseudo Mercator", AUTHORITY["EPSG","1024"]],
  PARAMETER["semi_minor", 6378137.0],
  PARAMETER["latitude_of_origin", 0.0],
  PARAMETER["central_meridian", 0.0],
  PARAMETER["scale_factor", 1.0],
  PARAMETER["false_easting", 0.0],
  PARAMETER["false_northing", 0.0],
  UNIT["m", 1.0],
  AXIS["Easting", EAST],
  AXIS["Northing", NORTH],
  AUTHORITY["EPSG","3857"]]

And when i publish layer and go to preview it only thing i am shown is IMG_3298.jpg

I calculated that 1m is represented with 38 pixels and as i am using EPSG 3857 as native srs its default unit is 1m, pixel size is about 0.03m.
What might be the problem?


texts: Using "Select by Location" in ArcGIS ModelBuilder?

texts:
I need some assistance for my workflow in ArcGIS ModelBuilder. I have two datasets:

Points (unique ID field: "bldg_id")
Parcels

I am trying to select Parcels within a certain distance (e.g., within 100 m) from the Points feature. There are several entities within the Points feature. The selected Parcels may sometimes overlap with each other for certain Points. I then want to output the selected Parcels into a new feature class.

My model has been successful up to a point, and I need assistance.

I am able to output each Parcels into a new shapefile based on the "bldg_id" field (%VALUE%) in the Points feature. In other words, if there are 10 unique "bldg_id" values, I get 10 unique parcel outputs called "Output_BldgID_%VALUE%." However, I want to merge all these individual outputs into one single output with a new field ID that identifies the corresponding BldgID (%VALUE%).

I am using the "MERGE" tool, but I don't know how to generate the unique field ID. My full dataset is several hundred points, so dealing with 100+ feature classes would be unwieldy. As a result, I want one single feature class. Any suggestions?


texts: Using orthophotos as basemap in CartoDB?

texts:
I have seen a java example of how to use a raster overlay in CartoDB, but what about replacing the basemap entirely with orthophotos?

In the web GUI for selecting basemap there is an option for the users imagery, but I can't seem to get it working and when I try to upload a TIFF to my datasets I get a "file too large" error.  I think I am supposed to create a Mapbox account and upload the TIFFs there and then point to them from CartoDB?

I guess I am looking for a step-by-step walk through of how to accomplish going from an orthophoto mosaic in ArcMap to exporting it out for use in Mapbox, putting it there, and then pointing to it in CartoDB (which seems to be the easy part).


texts: I am trying to calculate NDWI from Landsat-8 Band3(range:-113 to 11785) and Band5(range -89 to 12960) but the range of the resultant NDWI image is not valid (range: -215.181 to 1151.02). What is the possible workaround?


texts: How To Understand Geological Directions?

texts:
I was reading a "Guide to Some Minerals and Rocks in Indiana," which was published in the 1950s, and in it there were these bewildering directions:

"Bluff along Mud Pine Creek 1.5 miles west of Rainsville, SW 1/ 4NE1/ 4sec. 29, T. 23 N., R. 8 W."

It is on p.38 of this pdf: https://scholarworks.iu.edu/dspace/bitstream/handle/2022/235/C04.pdf?sequence=1

How may those directions be made understandable?


texts:
I want to query the soil data via https://sdmdataaccess.nrcs.usda.gov/WebServiceHelp.aspx.

Input is the polygon coordinates.
Result: features (polygon), key, soil code, soil name, total acres, % of field, PI ...

How can I do it?


texts:
I want to automaticaly insert/update points via cron in cartodb.
I am using SQL API.

How can I generate the_geom row via sql api in WKB format without web-interface?

Like this:

0101000020E61000004260E5D0224F464078D503E621074440


texts:
I was reading a travel diary of an ancestor of mine who describes a cruise from Hamburg to Natal (South Africa) in 1860. At one point a rock island named Trinidad (supposedly only inhabited by animals) is mentioned in the text and I am curious about its location to trace the route of the ship as exactly as possible.

Does anybody know how to find this out? Since according to the text the ship is heading to the "cap" in a straight way four days later the island could be somewhat close to the coast of western South Africa but might also be further to the north - but probably still close to the African shore.


texts: ESRI World Light Gray Base Map JavaScript - Prevent zoom out no more than the full extent

texts:
It is my first attempt to use JavaScript and the ESRI products.
I am using the World Light Gray Base map for my first web GIS application.

I have a question about the map full extent method.
I would like to set the full extent size so users can't zoom out no more than my full extent setting. So the map always fully covers the map window. At the moment, if I zoom out the World Light Gray Base Map, the full size of the map is getting smaller so it covers a small part of my map window.

Could anyone helps me how to set up the maximum zoom out size so my map window is always covered by the world light gray map?

Best,


texts:
I have a dataset with parcels with a unique lotplan - lets say it has 10,000 records. I have a dataset which has owner info which has one or many owners for each lotplan lets say it has 11,0000 records. For example this would be an apartment complex or similar where there are many owners.

We need to join the parcels to owner by lotplan and where a single lotplan has multiple owners the lot plan spatial object should be repeated in the dataset for each.

So we should end up with a combined dataset of 11,000 records.

How can I do this in Mapbasic?

In QGIS I'd use Virtual Layers (Joining multiple records to single feature using QGIS?)

The brute force method maybe to rebuild the parcel table based on the owner info. So row by row go through the owner table by lotplan, find the spatial record in the parcel dataset, append it to a new table, fill in the fields based on the records in both datasets and then do the next one in the owner table. This will ensure one spatial object for each record in the owner dataset.

I hope I don't have to do the above and there's an easier way...


texts:
I made several point shapefile layers. The properties only differ in their assigned colours. Each attrubute table consists of a single field "id" and each record is null.

I now wish to assign a name (that the colour corresponds to) throughout each table. Can I do that in the id column, or is the id field something special that should be left alone? (I eventually plan to combine the tables then add other fields, I just want to start properly.)


texts: Google Satellite imagery gone wild with QGIS 2.18

texts:
I'm using the OpenLayers plugin (version 1.4.3) for my satellite imagery.  When I'm at 1:565 or less (the number on the right is bigger) everything is fine.  See image below.

When I'm at 1:564 or greater the satellite image is incorrectly aligned / scaled.  In this case, I change the number is the bottom readout (circled in red) to avoid moving the map in any way.

I am using an API key.  And can zoom into 1:282 before the image goes away (what appears to be the plugin's way of letting you know the imagery isn't supported at that scale.

I haven't noticed this behavior with Bing imagery, or Google Hybrid, but can't zoom in as far either either case.


texts:
I hasve an XML document originated from this site. It contains location and data for resturants and bars. How can I import those locations theirs data and meta-data to arcmap?

I am still relatively new to arcgis.


texts: Symbology of Outlines are Overlapping Instead of Layers

texts:
I have three layers, three symbologies. If I create a group layer, the symbols are rendered differently than if I create a merged layer of all three layers.

Group layer:

Merged Layer:

Is there a way to get the merged layer symbology to look like the group layer symbology, where the symbol's outlines are layered, instead of combining. I understand that with merged layers, that the rules of layering may not apply, but I want to get rid of the noise from the outlines, while keeping them.

I am using ArcMap 10.6.1


texts: Merging shapefile with dataframe and plotting using R

texts:
I am using R.
I have a shapefile that gives me the provincial boundaries for a state and I want to use another dataset to combine the aggregate total enrolment in schools by 'district' to the shapefile. Since it is one-to-one merge, I thought a simple cbind should work but the file that is created is a list and I cannot map the density of enrolment onto the shapefile using the tmap function.
combined<- cbind(punjab, emis_data, by="ADM2_EN")

The other alternative is to use the merge command with aggregate function but that gives me all NAs in the total enrolment column.
Why might this be happening?


texts: Making multiple convex hulls in a GeoDataFrame based on a categorical attribute in GeoPandas

texts:
Suppose I have a GeoDataFrame as follows:
    | cat    | geom
----+--------+-------------
  0 | 0      | POINT(0, 0)
  1 | 0      | POINT(0, 1)
  2 | 0      | POINT(1, 1)
  3 | 0      | POINT(1, 0)
  4 | 1      | POINT(5, 5)
  5 | 1      | POINT(5, 6)
  6 | 1      | POINT(6, 6)
  7 | 1      | POINT(6, 5)

I want to make a GeoDataFrame with convex hulls for the point geometries for each category, while retaining the cat column to indicate which category the hull corresponds to. Can I do this in one line? I was thinking something like:
import geopandas as gpd
gdf = ...
gpd.tools.collect(gdf.groupby(gdf["cat"])[gdf.geometry.name]).convexhull

I know that this isn't right but I'm not sure about the grouping syntax required to achieve such a result.


texts: Don't route through bollards, concrete blocks, etc with osm2po

texts:
Is it possible to prevent road segments being used when these are blocked by bollards, concrete blocks or the like?

I have one example where there the road Solhøjvej is blocked with a so called concrete block as you can see here

http://www.openstreetmap.org/?lat=57.010248&lon=9.936238&zoom=18&layers=M

Here is a photo of the "block":

https://maps.google.com/?ll=57.010343,9.934331&spn=0.000929,0.00284&t=m&layer=c&cbll=57.010252,9.934346&panoid=SrWxGTFwoDbDZeDgXHl8vg&cbp=12,10.73,,0,10.62&z=19

I’m exporting the Openstreetmap to PostgreSQL using osm2po and performs routing, but this does not know of these concrete blocks or bollards or stuff blocking the road.

I’m aware that this is perhaps not a simple task, to detect and split road segments when such obstacles occur, but that would be very useful. Am I just missing an option?


texts: How to preserve GPX track point (trkpt) info when performing ogr2ogr with sql ST_union operation

texts:
I am performing GDAL ogr2org operation with ST_union(Geometry) on tracks from GPX file to split those at intersection points. Track points have elevation and depth info. In GPX track point entry looks like this:

    <trkpt lat="51.1012332738" lon="20.916904">
      <ele>130</ele>
      <extensions>
        <depth>12.3</depth>
      </extensions>
    </trkpt>

The command being executed is this:

ogr2ogr -nlt MULTILINESTRING -dsco GPX_USE_EXTENSIONS=YES -lco FORCE_GPX_TRACK=YES  -f GPX "%output_gpx%" "%input_gpx%" -dialect sqlite -sql "select ST_union(Geometry) from tracks"

The problem is that resulting GPX file looses depth and elevation data.
The question is how to keep depth and elevation data in resulting GPX file after ST_union(Geometry)?

Sample GPX file: DOWNLOAD


texts:
I imported an OSM file into a Postgres DB using the Osmosis tool. Now, I need to extract some data but I can't find any schema description (I'm really a newbie in GIS).

Where can I find it?

Do you know if there is a tutorial to learn how to query an Osmosis Postgres DB?

I need to extract the names of all streets available within a given bounding box.

How can I do that?


texts: Snapping existing raster to another existing raster using ArcGIS for Desktop?

texts:
I have got two existing rasters that have the same cell size (5x5m). The corner of cells from raster A are located exactly at the centres of cells in raster B. I want to place raster B exactly on top of my raster A.

I have tried these options:

the tool 'resample' on raster B, while maintaining the same cell size and snapping to raster A.
the tools 'raster to points' and consequently 'points to raster' on raster B, again, maintaining the grid size and snapping to raster A.
the tool 'copy raster', snapping raster B-copy to raster A.

All options have the same result: raster B is still exactly where it lay before, with the corners of cells in B at the centre of cells in A.

Whenever Google this problem, I end up at 'environment settings' when making the grid for the first time - but I do not have the data on which the grid was based.
Do you have any suggestions on how I could solve this, without a need for the original data?

I am working in ArcGIS 10.2.1. I've got a basic license, but might be able to do an upgrade to Standard or Advanced if necessary.

Edit: the raster files are .asc files.


texts: GEE skipEmptyTiles has no effect: export produces a bunch of empty tiles

texts:
I am trying to get Google Earth Engine to export many polygons of an Image (thousands all around the world, training data for an ML model) into .tif files on Google Cloud Storage, using a single GEE export task, because given how many polygons I want to export, having 1 export task per polygon results in the whole thing taking days.
As a prototype, I'm trying only 2 polygons, and exporting simply some pre-existing Landsat data.
The code below tries to do that by clipping the Image to the union of the polygons, and exporting a region which is the bounding rectangle of them all, but with skipEmptyTiles: true.
Unfortunately, it does not have the required effect - despite skipEmptyTiles, the code produces a bunch of empty tiles (identical completely black .tif files). The non-empty tiles look correct.
Is this a bug in GEE? Am I doing something wrong?
var landsat = ee.Image('LANDSAT/LC08/C01/T1_TOA/LC08_123032_20140515')
                  .select(['B4', 'B3', 'B2']);
// Brighten a bit
landsat = landsat.add(landsat).add(landsat).add(landsat);

var r1 = ee.Geometry.Polygon([[116.256, 40.515], [116.343, 40.544], [116.313, 40.632], [116.146, 40.632], [116.135, 40.565]]);
var r2 = ee.Geometry.Polygon([[116.819,40.333], [116.926,40.362], [116.852,40.460], [116.662,40.456], [116.668,40.335]]);

var union = r1.union(r2);
var bounds = union.bounds();

// Produces 28 .tif's totaling 32MB and tiling the region 4x7.
// Despite skipEmptyTiles, 20 of the tiles are completely black.
Export.image.toCloudStorage({
  image: landsat.clip(union),
  bucket: 'ekirpichov-gee-test',
  description: 'union',
  scale: 100,
  region: union.bounds(),
  skipEmptyTiles: true
});

// Debugging
Map.centerObject(landsat);
Map.addLayer(r1, {color: 'FF0000'});
Map.addLayer(r2, {color: '00FF00'});
Map.addLayer(landsat.clip(union), {}, 'landsat');


texts: Error of Post request of the WFS-T on the https protocol

texts:
Tell me why the following request does not work for me via the http protocol:
https://test.geoportal.samgau.com/geoserver/AkmolaRegion/wfs
Body:
<GetFeature xmlns="http://www.opengis.net/wfs" service="WFS" version="1.1.0" outputFormat="application/json" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://www.opengis.net/wfs http://schemas.opengis.net/wfs/1.1.0/wfs.xsd"><Query typeName="AkmolaRegion:001_1132201_ads_lin" srsName="EPSG:3857" xmlns:AkmolaRegion=" "/></GetFeature>

But if the protocol is changed, then everything works:
http://test.geoportal.samgau.com/geoserver/AkmolaRegion/wfs


texts: Where to find information about LinkedGeoData ontology?

texts: LinkedGeoData is an ontology with information derived from openstreetmap. I want to query linkedgeodata to get administrative areas and points inside those areas, but I would like to know a bit more about the ontology before starting. I already ran a query to get all the predicates, what I want is real documentation about each predicate.


texts: Using Google Transit Layer on Google My Maps

texts:
I've played with this a fair amount and have determined there isn't a simple way to render Google Transit Layer onto a map I've created through Google My Maps. These are the solutions I've come up with which all lead to dead ends:

Use the Google Maps API to render the Google Transit Layer and my Google My Maps map. The issue with this is I've read My Maps isn't accessibly through the API, and even if it was, I wouldn't be able to use the Google Maps phone app to view the two maps, I'd have to navigate to a page I'm hosting.
I've found 3 stores of public GTFS feeds, none of which have the cities I'm looking into (but Google does). If they did have the cities I wanted, I did find Google's GTFS to KML development tools in Python. Then I would have been able to import the KML to create a layer on My Map, but no such luck.
Use a 3rd party service such as Maps Channel to composite the two maps together. I can't seem to get this service to work, but it stands to reason that if it has access to My Maps that the API should allow access, but who knows.

If someone comes upon a solution to this which doesn't require manually creating a GTFS or KML file for a city or something equally labor intensive, let me know.


texts: Selecting columns in from an attribute table in ArcPy using a list (rather than slices or indices)?

texts:
I have seen this answer to Calling specific columns of attribute table using ArcPy?, but it uses a slice to select columns.

I would like to select them using a list.

Example code that doesn't work:

import arcpy
select_cols = ['FID', 'Shape', 'OBJECTID'] # the list of columns
    for field in field_list[select_cols]:
        print(field.name)

  Traceback (most recent call last):
        File "", line 1, in
      TypeError: list indices must be integers or slices, not list

A different attempt

field_list2 = field_list[select_cols]

  Traceback (most recent call last):
        File "", line 1, in
      TypeError: list indices must be integers or slices, not list


texts: How do I get the distance meter value between two geometries in PostGIS?

texts:
I have a simple question about calculating distances in PostGIS.

I would like to get the distance between two geometries. I am using this sid : 4269 in meters
what i am doing now is this : ST_Distance((a.geom,b.geom)) FROM ... but I am getting result in degrees. I think that I should work with geography but how can I cast a geom to a geography? I tried with (a.geom::geography) but i am getting an error.

What can I do to get my result in meters?

Thank You


texts: How to present three spatial variables in one map?

texts:
I want to assign 3 category of variables in a single map.
In other words, I have to produce a cartographic visualization of the damage including a legend listing:

damaged surface area in the residential areas (in km2),
number of living units and number of people affected, and
damage in terms of euroes.

I just can present two of them by "graduated colors" and "graduated symbols".
I cannot add third variable, in a meaningful way.

How I can present all the 3 variable in a same snapshot simultaneously?


texts: Creating gradient-filled rectangle in PyQGIS layout?

texts:
I have created a print layout with some of the layers of my QGIS Project using PyQGIS. I am trying to add a gradient filled rectangle to the layout. In the GUI, I would use the "Add Shape" --> "Add Rectangle" tool, drag the shape to the desired location, and use the "Item Properties" to fill with a gradient. However, I am having difficulties reproducing the shape programmatically.

To create the Layout, I use the following code:

project = QgsProject.instance()
manager = project.layoutManager()
layout = QgsPrintLayout(project)
layoutName = "Some_name"

# The next step allows the script to run multiple times without naming issues
layouts_list = manager.printLayouts()
for layout in layouts_list:
    if layout.name() == layoutName:
        manager.removeLayout(layout)

layout.initializeDefaults()
layout.setName(layoutName)
manager.addLayout(layout)

After having added my layers, I tried to add a polygon:

# Make a polygon from coordinates (in layout)
# The coordinates are in the layout units
xmin = 20
xmax = 40
ymin = 20
ymax = 60
polygon = QPolygonF()
polygon.append(QPointF(xmin, ymin))
polygon.append(QPointF(xmax, ymin))
polygon.append(QPointF(xmax, ymax))
polygon.append(QPointF(xmin, ymax))

# Create the polygon from nodes
polygonItem = QgsLayoutItemPolygon(polygon, layoutname)

# Add to the layout
layout.addItem(polygonItem)

# Style properties
props = {}
props["color"] = "green"
props["style"] = "solid"
props["style_border"] = "solid"
props["color_border"] = "black"
props["width_border"] = "0.26"
props["joinstyle"] = "miter"

# Create symbolfiller from properties
symbol = QgsFillSymbol.createSimple(props)
polygonItem.setSymbol(symbol)

This approach allows me to create a filled polygon. I cannot find the method to fill the polygon with a gradient, and I haven't yet found the list of options that can be set in the parameter dictionary. I also tried to replace the QgsFillSymbol with QgsGradientFillSymbolLayer, but to no avail.


texts: Find single line strings in MultiLineStrings

texts:
I have a layer of line geometries in QGIS (version 2/Dufour) and want them all to be singlepart linestrings. However, some of them are multipart geometries (MultiLineStrings in a GeoJSON export). The MultiLineStrings sometimes have small gaps, but not always. Sometimes the geometry end points are in the same position.

How can I highlight/select the individual geometries so that I can easily find the gaps?

Once I can see where the gaps are I'd like to treat them manually, e.g. close the gap manually and merge the geometries.


texts: GeoServer-Mviewer : point layers over surface layers

texts:
I am working on an interactive map designed under Mviewer package and using GeoServer to save symbols and layers. I have a specific issue on it: I managed to generate all layers correctly in my interactive map, but I would like to control their overlaying.

In fact, I would like to decide that point layers must always be over linear layers and zonal layers, as they recover the rest of the layers when you click on them. Is there a specific command I could use in, either GeoServer, or my HTML/XML files which control my files for Mviewer? I know there is the command z-index in HTML and if possible, apply it between all my different layers in the map.


texts: Build RemapValue table name in Python script and then access the Remap Table

texts:
I have a bunch of rasters that need to be reclassified. I am not sure the best way to create RemapValue tables, but I made them in Excel. I then dumped them in the Python script and they are now hidden objects.
SiteA_remap = arcpy.sa.RemapValue([[1,100], [2,200]])
SiteB_remap = arcpy.sa.RemapValue([[1,10], [2,20]])

I can see them if I type SiteA_remap; the table will show up in the Python console.
But if try to build SiteA_remap like in the code below, my_remap = 'SiteA_remap1' and this will not bring up the remap table, just a string.
How do I call this Spatial Analyst object (remapTable) without typing the name in the console? Because if I type my_remap = SiteA_remap  it works and my_remap is a now a remapTable. How do I build the name and call it in my code?
import os, arcpy, arcpy.sa
from pathlib import Path

inputRasters = r"D:\FOO\BAR"

#list to hold rasters
fooList = []

for path, dirs, files in arcpy.da.Walk(inputRasters, datatype = "RasterDataset"):
    for f in files:
        fooPath = os.path.join(path, f)
        fooList.append(fooPath)

        thatName = Path(f).stem
        fooName = thatName.split("_",1)[0]

        inPath = os.path.dirname(fooPath)
        inDir = os.path.basename(inPath)

        #my_remap = os.path.join(inDirName + "_remap")
        my_remap = inDirName + "_remap"
        # do more stuff


texts:
How is it possible to make the map zoom to a specific point like an address that the user wants to find?

Codebehind:

private void btnSearch_Click(object sender, RoutedEventArgs e){
    var queryTask = new QueryTask();
    queryTask.Url = /*URL*/;
    queryTask.ExecuteCompleted += new EventHandler<QueryEventArgs>(queryTask_ExecuteCompleted);

    string[] addressInput = txtBoxAddress.Text.Split(' '); //Here goes the address in

    var query = new Query();
    query.ReturnGeometry = true;
    query.OutFields.Add("*");

    if (addressInput[0] != null)
        query.Where = String.Format(/*SQL-Statement*/, addressInput[0], addressInput[1]);

    queryTask.ExecuteAsync(query);
}

void queryTask_ExecuteCompleted(object sender, QueryEventArgs e){
    //throw new NotImplementedException();

    var graphicsLayer = (ESRI.ArcGIS.Client.GraphicsLayer)Map.Layers["graphic"]; //A little circle marks the address on the layer
    graphicsLayer.ClearGraphics();

    foreach (var item in e.FeatureSet){
        item.Symbol = null;
        graphicsLayer.Graphics.Add(item);
    }
}

Code in XAML:

<esri:GraphicsLayer ID="graphic" x:Name="graphic">
    <esri:GraphicsLayer.Renderer>
        <esri:SimpleRenderer>
            <esri:SimpleMarkerSymbol Color="Cyan"></esri:SimpleMarkerSymbol>
        </esri:SimpleRenderer>
     </esri:GraphicsLayer.Renderer>
 </esri:GraphicsLayer>

I'm developing in visual studio with C# (WPF)


texts:
I'm not sure if I've misunderstood how MBTiles and SlippyMap tilenames go together or if there is a problem when I'm doing an export with TileMill.

The issue is that a tile showing the same region - having the same bounding box - is having different tile_column/tile_row in MBtiles and in x/y in OSM.

For example the tile in OSM haves the z/x/y values of 12/2143/1406. But in my MBTiles file it haves the zoom_level/tile_column/tile_row values 12/2143/2689.

Any idea what is wrong with my export? Or am i misunderstanding naming in MBTiles?


texts: Do 1KB Shapefiles with bigger .shp.xml contain any usable data?

texts: I received this folder from a county's GIS department.  Am I right in telling them that there is nothing in it?  They seem convinced that everything is included, but all I see is that there is Metadata associated with the Shapefile that isn't contained within the Shapefile itself.  Am I doing something wrong or is the fault on their end?


texts: What is the formula for calculating world coordinates for a given LatLng in Google maps?

texts:
I want to calculate a given LatLng to the world coordinate .

I know I can use the internal projection of Google maps like this:

var map=new google.maps.Map(...);
map.getProjection().fromLatLngToPoint(latlng);

But I want the formula to calculate, rather then the result.


texts: ArcGIS Model Builder perform split operation on each selected feature and delete identical on collection

texts:
I have created a model to select parcels and create segments for each parcel. If two parcels are sharing the same segment, delete duplicate. My model is attached.

My question is: Before appending or deleting identicals, I need to perform the action on a whole. The model I created is executing for each feature. Please advice..!


texts: WMS layer doesn't show after closing and reopening the browser (Geoserver)

texts:
I am adding a WMS service to ArcGIS Map Viewer using the service URL and by providing my username and password. After I sign in, I am able to choose the layers I want to add. When I do this initially, the layers I just added render in the map window. However, after saving the map and closing the browser, next time I reopen the browser and the map I saved, the layers no longer render. The legend doesn't render either as shown in the picture below.

I used Google Chrome and Internet Explorer, and the same thing happens for both browsers. This makes me believe that the issue is with the geoserver configurations. I also don't believe it's a CORS issue, because CORS is enabled in web.xml, and I also tried disabling it in Chrome, but that did not help.
Has anyone experienced this and was able to find a solution?


texts: ComboBox Field not updating when layer is selected

texts:
I try to develop my plugin. The problem is that when the layer is selected in upper comboBox (QgsMapsLayerComboBox), in lower two boxes (Latitude, Longitude: usual ComboBoxes) shall appear fields correspoding with selected layers.

The problem is, that when I switch the layer in upper field, I can only choose the fields which are from the other layer...

On this picture, the fields "ID", "LONGITUDE", "LATITUDE" belongs to previously selected layer "coordinates" visible on 1st photo.

Here's my code of def_run of my plugin:

def run(self):

    self.dlg.layerCombo2.clear()
    self.dlg.layerCombo3.clear()

    layers = self.iface.legendInterface().layers()
    layer_list = []
    for layer in layers:
        layer_list.append(layer.name())
    self.dlg.layerCombo1.addItems(layer_list)
    self.dlg.layerCombo1.setFilters(QgsMapLayerProxyModel.PointLayer)

    def field_select():
        self.dlg.layerCombo2.clear()
        self.dlg.layerCombo3.clear()
        selectedLayerIndex = self.dlg.layerCombo1.currentIndex()
        selectedLayer = layers[selectedLayerIndex]
        fields = [field.name() for field in selectedLayer.pendingFields()]
        self.dlg.layerCombo2.addItems(fields)
        self.dlg.layerCombo3.addItems(fields)

    # This connects the function to the layer combobox when changed
    self.dlg.layerCombo1.currentIndexChanged.connect(field_select)
    self.dlg.show()

`


texts: I'm looking to programmatically convert UTM to Web Mercator (a.k.a.  Spherical Mercator, EPSG:3857, EPSG:3785, OpenLayers:900913). I'm looking for the algorithm, or something that runs on iOS (Swift, Objective-C or even C). The GDAL C library is not an option to due its huge size (this will be running on a mobile device).


texts: How to speed up map display performance?

texts:
I was wondering if there is any trick I can use to speed up zooming in/ out or panning through my maps in Arcmap. I usually work on county-wide studies in which I have lots of base layers, imagery and annotation which slow me down very much. I try to keep unnecessary layers off but can't think of anything else. I am using ArcGIS Desktop so I suppose I can not use caches.
I have not tried it yet but I guess I can use Dice (Data Management) to increase my speed. Any thought/ comment?


texts: Is elevation (height relative to sea level) constant, despite sea level rising?

texts:
Is a 10,000' peak still a 10,000' peak if sea levels have risen a foot? Is it a 9,999' peak now, or is elevation just a constant?

Seems odd that a patch of land a few feet above "sea level" could sit underwater.


texts:
I am currently working with ArcGIS 10.2.2 and need a script to create a tool (if possible). I have USGS 7.5' Quad Maps as my base and a series of UTM coordinates representing sites. I would like the tool to accept an Easting and Northing. Then zoom the map to the area and put a cross-hair in the location. Is this possible by using arcpy scripting?

If so, I am open to suggestions on where to start or even how to complete this task.


texts:
I have a map of trace gas (NO2) distribution over big cities and other anthropogenic objects. The purpose is to estimate the area of distribution.  I attached the sample image to illustrate what I do.

Is there some tool in ArcGIS to carry out such an operation?

I have just one multi-stage idea: exclude low values pixels---- convert new image into polygons----estimate the area of each polygon. Maybe there is some more elegant solution?


texts: How do I get the layer by name using PyQGIS 3?


texts: Selecting WMS layers from GeoServer in OpenLayers application?

texts:
I am making a weather application and I am almost there!

I can show my GeoServer WMS-layers in my browser using OpenLayers 5. What I want to accomplish is that the user can switch between different types of 'weather'. (The user can choose if it wants to see the wind, the rain or solar for example.)

How can I achieve this?

At the moment I made GeoServer group layers for all the different weather types. With a map baselayer and then the weather specific layers on top of that.

But how do I make the user choose what to see?

With some sort of menu in the corner.

I presume I have to do this in OpenLayers 5 and not GeoServer, but surely one of both has to have the ability to do this. I can't really find it myself.


texts: I am completely new to this and have never worked with shapefiles. What I'm trying to find out is whether it is possible, using the TIGER data in conjunction with the ZCTA files, to determine the physical boundaries of a ZCTA and once that step is complete, to generate a list of random geographic points (LAT/LON pairs) that fall within that ZCTA boundary.


texts: d.rast.edit not yet implemented in the WxGUI

texts:
Am using GRASS to modify a raster and tried using d.rast.edit but am getting the following error when using it in the WxGUI

Command 'd.rast.edit' not yet implemented in the WxGUI. Try adding it as a command layer instead

But am not sure how to add it as a command layer. I looked for information here
d.rast.edit error when attempting to edit a raster

but am not sure if it is the same issue. I also looked for information in the manual
https://grass.osgeo.org/grass76/manuals/d.rast.edit.html

But whenever i type d.rast.edit in the console in get the same error.

Does anyone knows how to address this? or an alternative methodology to modify a raster in GRASS?

Am using GRASS 7.4.2 in Windows


texts: Why Sentinel-2 bands have different spatial resolution

texts:
Sentinel-2 has 13 bands with 3 distinct resolutions.

Band 1 60m
Band 2 10m
Band 3 10m
Band 4 10m
Band 5 20m
Band 6 20m
Band 7 20m
Band 8 10m
Band 8A 20m
Band 9 60m
Band 10 60m
Band 11 20m
Band 12 20m

This page just gives how they separate the spectrum.
But does not elaborate on how these bands have different resolution (the technology behind it).

I am curious about the advantage of this variation in spatial resolution (since Landsat 8 has only 30m)?


texts: Running process in ModelBuilder before loading shapefiles?

texts: I have a model which I would like to check if shapefiles exist or are up to date and download them if they are not before processing them. I wrote the script to do the checking. In the script, I use arcpy.SetParameter() to pass in the file location of my downloaded shapefiles once checking is done and I pass them as derived outputs parameters in the script's properties. My problem is that I am not able to force ModelBuilder to run the script before trying to do anything else. The output shapefiles which I use for processing are none existent until the script is run so I get errors telling me that fields do not exist and such (because of processing steps down the line). Simply put the Model Thinks that I want it to run without input shapefiles when I just want it to wait for them to be obtained. I have tried preconditions but that has not worked, maybe I did it wrong. The answer must be in ModelBuilder, or if it is a python script, it has to be one that I can add to the model, not export my whole model as a script.


texts: How to replace/update authorization key for an ArcGIS Enterprise Geodatabase

texts:
I need to switch the authorization key of an ArcGIS Enterprise Geodatabase (in Postgresql). I am on 10.3 and based on the scarce documentation I have found, ESRI is moving tools that allow such administration to arcpy and the GUI, however, the particular tool that is on the docs in pre 10.3 era, sdesetup is not available with ArcGIS installation (Server or Desktop). So basically I am wondering how to switch keys.

I have tried the "Enable Enterprise Geodatabase" tool with no avail saying that it is already authorized.


texts: High number of unique user defined views

texts:
We're just at the stage of starting to developing a GIS application accessible through a web front end. We're thinking of using the boundless stack which includes GeoServer and PostGIS.
One of our key requirements is having bespoke combinations of objects visible to each user. These can combinations of any type of object, which I expect would all exist across several layers in a mapping system.
My question is, how best to manage this high number of unique views? There can be a very high number of objects that make up one view, while another might have quite a small number, so I'd expect some form of caching will be necessary for performance and drawing them "on the fly" will be too slow. But my understanding of caching is that it will only be beneficial or practical if the same view is seen by number of users.
We have around 100,000 different objects (about half are vectors and half points) to give you a clue of the scale of data and some users will have access to all of it while others only a section.
The definitions of what objects are available to each user are already stored as queries in a relational database, so can potentially pass them over as a list of object ids if that's the right way to do things.
I have two ideas about how to solve this:

Break these definitions down into tiles at a certain zoom level, so all the unique user access will be a list of tiles. This will have the disadvantage that the smallest section we can filter by is the size of a tile and other objects we might not have wanted them to have access to are shown.
Have a general low detail view that everyone sees and some user interaction is performed to see which part they want to see in detail and this can be drawn "on the fly". This will mean only a small number of objects is drawn at the time so performance isn't too bad but will also mean they can't see everything they have access to at one time.

Or a combination of the two.

I'm really after a pointer in the right direction as am pretty new to developing G.I.S. software rather than a complete solution. I just want to check whether I'm on the right track.
Thanks very much.
edit: changed client to users


texts: What is the right SQL query for this case?

texts:
I want to use input string in part of the SQL query, and the code has been posted below.

import arcpy
import os

in_features = arcpy.GetParameterAsText(0)
target_county = arcpy.GetParameterAsText(1)
out_path = os.path.dirname(in_features)
out_name= "%s.shp" % target_county
where_clause=""""CNTY_NM" = %s""" % target_county
arcpy.FeatureClassToFeatureClass_conversion(in_features, out_path, out_name, where_clause)

When I run the script, it gives the error.

I don't know what is wrong with the code. Any suggestions?

I am using ArcGIS 10.3.1 Desktop Advanced License.


texts: Removing time from date field labeling in VBScript

texts:
I'm not at all experienced with writing VBScript for labeling purposes. Here is my code so far:

Function FindLabel ( [Date_of_Birth] )
  if not isnull ( [Date_of_Birth] ) then
         FindLabel = "DOB: "  &  [Date_of_Birth]
            FormatDateTime [Date_of_Birth],????
  else
         FindLabel = ""
  end if
End Function

Basically, I want to label features that have a DOB with the DOB from the field, and if that field is null I don't want it to label. That part works great, but I noticed that some of the dates also had a time tacked on to the end. I would like to format that part out if possible. The line that reads FormatDateTime [Date_of_Birth],???? is what I found searching online, but nothing I've entered after that comma works. I'm curious if I'm on the right track or not with how I'm trying to format the date.

Thanks,
Chase


texts:
I'm trying to convert closed polylines into polygons using ogr2ogr from dxf to mssql. Here It says that polylines are converted into POLYGON if they are closed, but I always get the type LINESTRING instead of POLYGON in the geometry type column. I also tried to convert DXF into Shapefile and into a MySQL database with same results. Converting Shapefiles with polygons into databases works well, so I think the problem is in the gdal dxf driver.

Can someone explain to me what I am doing wrong?

Here is the command I'm running

ogr2ogr -update -append -f MSSQL "MSSQL:server=.\SQLEXPRESS;database=testdb;trusted_connection=yes" test.dxf


texts: Spatial Join to count points in polygons in lat/lon in PostGIS?

texts:
I have a table, called pets, with lats and lng coordinates and variables and a table that was created by uploading a shapefile, called city, with columns left, right, top, bottom, lat, lng, and  geom.

    lat    lng   cats  dogs
   -59.4   42.1   2     0
   -58.1   44.3   1     1
   -59.7   43.6   0     3
   -59.8   42.0   2     1

I'm wondering how I would construct a query that will sum up the number of cats and dogs within each boundary. Usually I would do a spatial join in QGIS but my dataset is too large and it crashes before successfully joining the two.

I am pretty familiar with QGIS and Postgres, but very new to PostGIS. I've constructed a query but unsuccessful in getting the answer I needed:

select sum(pets.dogs)
from pets, city
WHERE ST_Within(pets.dogs, city.geom);

How do I do this in PostGIS?


texts: Getting point from known points to specific distance of polyline using ArcObjects?

texts:
I want to know the coordinates of a polyline from specific point and distance.
I have attached image to clear the specification.
say, I have a points B(2,3). I want to know the coordinates of A and C from 2 feet distance of B point. I have tried this:

public IPoint GetPointAtFrom(double distance, Coordinate givenCoordinate)
    {
        ArcGeometry.IPoint outletPoint = new ArcGeometry.PointClass();
        ArcGeometry.ICurve curve = _polyline;
        curve.QueryPoint(ArcGeometry.esriSegmentExtension.esriExtendTangentAtFrom, distance, false, outletPoint);
        return new Point(outletPoint);
    }

but it gives from point or others.

Any idea? I am using .net in arcgis 10.0


texts: How to choose cell assignment method in GDAL/rasterio rasterize functions?

texts:
I have a big (2GB) landcover shapefile I want to convert to a GeoTIFF file. For this purpose, I am using the gdal.Rasterize() function with the following options:
rasterize_options = gdal.RasterizeOptions(
    outputType=gdal.GDT_Byte,
    xRes=1000,  # 1km
    yRes=1000,  # 1km
    noData=255,
    allTouched=True,
    attribute='CODE'  # Shapefile field with each landcover class
)
gdal.Rasterize(dst, src, options=rasterize_options)

Although this works, it is not clear what cell assignment method the function is using to assign the output value to each pixel. Here is an example of the input shapefile and the output raster (with each color representing a different class):

I want to specify a cell assignment method where the output value in a given pixel corresponds to the class of the polygons that occupy most of the area in the pixel. Thus, the output raster shown above would end up being mostly orange and probably just a few pixels, if any, would be green.
This seems to be possible in other tools such as ArcGIS's Polygon to Raster, where one can specify the cell_assignment parameter:

The method to determine how the cell will be assigned a value when
more than one feature falls within a cell.

CELL_CENTER —The polygon that overlaps the center of the cell yields the attribute to assign to the cell.
MAXIMUM_AREA —The single feature with the largest area within the cell yields the attribute to assign to the cell.
MAXIMUM_COMBINED_AREA —If there is more than one feature in a cell with the same value, the areas of these features will be combined. The combined feature with the largest area within the cell will determine the value to assign to the cell.

However, after reading the documentation, GDAL does not seem to have the option to specify the cell assignment method. I also took a look at the rasterio.features.rasterize() function but my best guess is that this function is just a wrapper of GDAL's rasterize function and does not offer the option to specify the method either.
As I don't have access to ArcGIS and would rather use open source tools for this, the only workaround I have thought of is rasterizing the shapefile to a smaller cell size (30m maybe as the shapefile was produced using Landsat imagery) and then resampling the output using a mode algorithm. However, this does not seem like a good option given the extent of the original shapefile. What are my options here?


texts: Creating color ramp based on quantitative range

texts:
I have a set of polygons with correlation values. I would like to color each polygon based how strong their correlation value is where 0=lightest and 1=darkest.
This means if I have a set of three polygons whose correlation values are: 0.4, 0.41, 0.42, then their colors should be relatively similar. However, when I set the layer symbology to a graduated color ramp, ArcMap splits these values into their own class where 0.4=light, 0.41=medium, 0.42=dark.
In other words, how would I create a color gradient for values based on a range rather than distinct classes?


texts: Outgoing connection to an IRC channel from GeoNetwork

texts:
I'm using a GeoNetwork 3.10 service installed on an Ubuntu 18 server, and noticed that this server often experiences what seems to be Denial-Of-Service attacks, overwhelming Apache and Tomcat.
I then ran the following command to check the outgoing connections made by the server :
netstat -putw

I noticed that GeoNetwork apparently connects an external IRC channel :
tcp6       0      0 geocatalogue.afric:http irc.efnet.nl:ircd       SYN_RECV    -

Is there any explanation for this ? It seems that it is a backdoor being exploited. Is there any possibility to block this suspicious script and to identify the part of the GeoNetwork code doing this ?
EDIT : I'm not 100% sure that the outgoing connection is performed by GeoNetwork (netstat doesn't list nor java nor the pid of the involved process) but the problem appeared after having upgraded GeoNetwork from 3.0 to 3.10....


texts:
I'm trying to compile liblas since I can't find any Mac executables out there(?). I have OSX 10.8 and XCode 4.6.3.

Since I'm not a programmer I get really confused when I follow the instructions on liplas.org . In those instructions older version of XCode is being used so with the new XCode I can't find anything..

I get stuck here

And if I just skip all the steps and try compile in XCode I just get this message:

  make -f
  /Users/oskarlin/Downloads/libLAS-1.7.0b2/xcode/CMakeScripts/ReRunCMake.make
  make1: * No rule to make target
  /Users/oskarlin/Downloads/libLAS-1.7.0b2/xcode/CMakeFiles/2.8.10/CMakeCCompiler.cmake',
  needed byCMakeFiles/cmake.check_cache'.  Stop. make: *
  [/Users/oskarlin/Downloads/libLAS-1.7.0b2/xcode/CMakeFiles/ZERO_CHECK]
  Error 2 Command /bin/sh failed with exit code 2

Has anyone here complied it for Mac and could perhaps give me an updated instruction?


texts:
Is there a geometry extension that can be enabled in PostGIS, in order to save Shapefile data in it?

I imported Shapefile  data in to PostGIS using shp2pgsql plugin, having the notion that PostGIS should automatically add geometry attribute to the data because it is a polygon data, but when I exported it to QGIS, the data appeared as just table with no geometry and I cannot work with it.

Please how can I resolve this issue?


texts: Getting ValueError while creating VRT using python

texts:
I am trying to create VRT file using python(GDAL). I am getting "ValueError: array larger than output file, or offset off edge".

How to solve this error?

Code written for the same is as follows.

import numpy as np
import gdal
from gdalconst import *
import osr
ds = gdal.Open('/home/vigna/Downloads/2740115314.GMODO-SVM12_npp_d20170824_t0721423_e0727209_b30176_c20170824140610010359_nobc_ops.h5')

subds = ds.GetSubDatasets()
#print(subds)
bt=gdal.Open('HDF5:"/home/vigna/Downloads/2740115314.GMODO-SVM12_npp_d20170824_t0721423_e0727209_b30176_c20170824140610010359_nobc_ops.h5"://All_Data/VIIRS-M12-SDR_All/BrightnessTemperature')
bt_array=bt.ReadAsArray()
#print(bt_array)
long=gdal.Open('HDF5:"/home/vigna/Downloads/2740115314.GMODO-SVM12_npp_d20170824_t0721423_e0727209_b30176_c20170824140610010359_nobc_ops.h5"://All_Data/VIIRS-MOD-GEO_All/Longitude')
long_array=long.ReadAsArray()
lat=gdal.Open('HDF5:"/home/vigna/Downloads/2740115314.GMODO-SVM12_npp_d20170824_t0721423_e0727209_b30176_c20170824140610010359_nobc_ops.h5"://All_Data/VIIRS-MOD-GEO_All/Latitude')
lat_array=lat.ReadAsArray()
def raster(ds):
    xsize = ds.RasterXSize
    ysize = ds.RasterYSize
    GeoT = ds.GetGeoTransform()
    Projection = osr.SpatialReference()
    Projection.ImportFromWkt(ds.GetProjectionRef())
    DataType = 'UInt16'
    X_dataset=long_array
    Y_dataset=lat_array
    return xsize, ysize, GeoT, Projection,DataType,X_dataset,Y_dataset
def Create_vrt(Name, Array, driver,xsize, ysize, GeoT, Projection,X_dataset,Y_dataset):
#for DataType:
    DataType = gdal.GDT_Float32
    NewFileName = 'bt3'+'.vrt'
    DataSet = driver.Create( NewFileName, xsize,ysize, 1, DataType )
        # the '1' is for band 1.
    DataSet.SetGeoTransform(GeoT)
    DataSet.SetProjection( Projection.ExportToWkt() )
    # Write the array
    DataSet.GetRasterBand(1).WriteArray( Array )
    return NewFileName
DataSet = ds
Band = DataSet.GetRasterBand(1)
# Open as an array.
Array = bt_array
NewArray = Array
print(NewArray)
xsize, ysize, GeoT, Projection, DataType ,X_dataset,Y_dataset= raster(ds)
driver = gdal.GetDriverByName('VRT')
NewFileName = Create_vrt('bt3', NewArray, driver,
                        xsize, ysize, GeoT, Projection,    X_dataset,Y_dataset)
#Create_vrt()
print("ok")


texts: arcpy.MakeFeatureLayer_management function not creating feature layer in ArcGIS Pro

texts:
I'm attempting to run a simple tool script in ArcGIS Pro 2.4 that creates a feature layer from an existing layer. The script creates the feature layer just fine if running the script in PyCharm, as well as the Python interpreter within ArcGIS (adds visible feature layer under Maps).

Why would the layer not be created and be visible in ArcGIS pro if running the script as a tool, compared to running it externally or in the interpreter?

Anyone else having this issue with 2.4 when running scripts as a tool?

No errors.

p = arcpy.mp.ArcGISProject("CURRENT")
a_map = p.activeMap
lyrs = a_map.listLayers()
# Checking layers within activeMap, they exist
for l in lyrs:
 arcpy.AddMessage("{0}: activeMap layer: ".format(l))

# Feature layer not created when used as a tool, but works in the Python interpreter & external PyCharm run
arcpy.MakeFeatureLayer_management(a_map.listLayers('countries')[0], 'test')


texts: Authoring map with advanced gradient symbol for ArcGIS web service?

texts:
I need to prepare a web service which will display the gradient of polygon points values. It has to look just like in this picture (please ignore the arrows):

Each polygon is constructed always by 9 points. The logical order of points on the polyline forming the polygon is constant. Values of points will differ. The angles and lengths of lines may differ. Basically, it will always form a trapeze.

First I played with splitting the 'large' polygon into 4 small ones using its central point. In this way, I could adjust the angles of gradients defining small polygons so that they formed this large polygonal gradient. There were however two issues - connections of the small polygons were looking rather ugly and most important - as the values of points may differ, it only suited the default scenario presented in the picture.


texts: Using Select by Location and Select by Attributes in single ModelBuilder model?

texts:
I'm working on a project involving wind turbines. I'm using ArcMap 10.2 and ModelBuilder to count the number of wind turbines within a 5km buffer zone of a given project (let's call it Project A), and then remove the turbines associated with Project A so I just have the "foreign" turbines within the buffer (e.g., those not associated with the project which establishes the median center of the buffer).

The image below shows a simple example involving just two projects, but I've got 58 different projects I need to analyze, some of which overlap with multiple wind farms which is why I am attempting to use ModelBuilder to automate the process.

The script does what I need when I step through the components individually (see figure below), but when I run the last process, it just returns an empty result.

I'm pasting the manual steps I'm trying to replicate in ModelBuilder below. I have a feeling I'm overlooking something obvious, but despite a lot of tweaking I can't figure it out.

Select by Location:

Selection method: Select features from
Target layer: ALL WIND TURBINES
Source layer: PROJECT A (projects saved as separate layers)
Spatial selection method "are within a distance of the source feature" 5,000 meters

Select by Attributes

Method: Remove from current selection
Where site_name = ‘PROJECT A’
Count using attribute table

Add value of "foreign" turbines to project shapefile

Start Editor and add value by hand to Fr_trb field


texts: Bing Aerial (OpenLayers) won't show below 1:2257

texts:
I'm having this problem with OpenLayers (plugin v. 1.3.5 - QGIS v. 2.6.1 Brighton) where it does not display the map when zooming in past 1:2257. If I zoom to 1:2556 (or anything below it), I won't have any background image from Bing Aerial.

This is a problem specially when I want to draw lines/polygons with more detail, or view information (points, mainly) when they are very close together.

Any idea why this happens, and any solution?


texts: postgis_tiger_geocoder nationscript.sh does not set the location of PSQL correctly

texts:
This generated nationscript includes the following line: PSQL=${PGBIN}/p (...)
I have tried 'umpteen' variations of replacements, including the logical one: the physical path to the psql file.

I am running the script as the postgres user who has permissions on the path to the psql file. This is an installation of enterpriseDB postgres (standard) on a big-honking mac running Yosemite. So the postgres user owns the /library/postgresql/9.5/bin folder.

And I have executed sudo su postgres in the folder where the nationscript.sh resides, and tried to execute it without that as well.

In any case the first use of the ${PSQL} macro is on line 15 through 22 or so. and that's the heart of the script and nothing will make these lines do anything but fail.

So it's probably some bash syntax arcanity that I can't deduce. But this is all 'supposed' to be generated script out of the postgis-tiger-geocoder extension.

TMPDIR="/volumes/metro/_data/gisdata/temp/"

UNZIPTOOL=unzip

WGETTOOL="/usr/bin/wget"

export PGBIN=/library/postgresql/9.5/bin

export PGPORT=5432

export PGHOST=localhost

export PGUSER=postgres

export PGPASSWORD=uggabuggababytellmewhyowhyowhy

export PGDATABASE=tiger

PSQL=${PGBIN}/p (...)"

cd /volumes/metro/_data/gisdata

wget ftp://ftp2.census.gov/geo/tiger/TIGER2015/STATE/ --no-parent --relative --

recursive --level=1 --accept=zip --mirror --reject=html

cd /volumes/metro/_data/gisdata/ftp2.census.gov/geo/tiger/TIGER2015/STATE

rm -f ${TMPDIR}/*.*

${PSQL} -c "DROP SCHEMA IF EXISTS tiger_staging CASCADE;"

${PSQL} -c "CREATE SCHEMA tiger_staging;"

for z in tl_*state.zip ; do $UNZIPTOOL -o -d $TMPDIR $z; done

for z in */tl_*state.zip ; do $UNZIPTOOL -o -d $TMPDIR $z; done

cd $TMPDIR;

${PSQL} -c "CREATE TABLE tiger_data.state_all(CONSTRAINT pk_state_all PRIMARY KEY (statefp),CONSTRAINT uidx_state_all_stusps  UNIQUE (stusps), CONSTRAINT uidx_state_all_gid UNIQUE (gid) ) INHERITS(tiger.state); "
${SHP2PGSQL} -D -c -s 4269 -g the_geom   -W "latin1" tl_2015_us_state.dbf tiger_staging.state | ${PSQL}
${PSQL} -c "SELECT loader_load_staged_data(lower('state'), lower('state_all')); "
    ${PSQL} -c "CREATE INDEX tiger_data_state_all_the_geom_gist ON tiger_data.state_all USING gist(the_geom);"
    ${PSQL} -c "VACUUM ANALYZE tiger_data.state_all"
cd /volumes/metro/_data/gisdata
wget ftp://ftp2.census.gov/geo/tiger/TIGER2015/COUNTY/ --no-parent --relative --recursive --level=1 --accept=zip --mirror --reject=html
cd /volumes/metro/_data/gisdata/ftp2.census.gov/geo/tiger/TIGER2015/COUNTY
rm -f ${TMPDIR}/*.*
${PSQL} -c "DROP SCHEMA IF EXISTS tiger_staging CASCADE;"
${PSQL} -c "CREATE SCHEMA tiger_staging;"
for z in tl_*county.zip ; do $UNZIPTOOL -o -d $TMPDIR $z; done
for z in */tl_*county.zip ; do $UNZIPTOOL -o -d $TMPDIR $z; done
cd $TMPDIR;

${PSQL} -c "CREATE TABLE tiger_data.county_all(CONSTRAINT pk_tiger_data_county_all PRIMARY KEY (cntyidfp),CONSTRAINT uidx_tiger_data_county_all_gid UNIQUE (gid)  ) INHERITS(tiger.county); "
${SHP2PGSQL} -D -c -s 4269 -g the_geom   -W "latin1" tl_2015_us_county.dbf tiger_staging.county | ${PSQL}
${PSQL} -c "ALTER TABLE tiger_staging.county RENAME geoid TO cntyidfp;  SELECT loader_load_staged_data(lower('county'), lower('county_all'));"
    ${PSQL} -c "CREATE INDEX tiger_data_county_the_geom_gist ON tiger_data.county_all USING gist(the_geom);"
    ${PSQL} -c "CREATE UNIQUE INDEX uidx_tiger_data_county_all_statefp_countyfp ON tiger_data.county_all USING btree(statefp,countyfp);"
    ${PSQL} -c "CREATE TABLE tiger_data.county_all_lookup ( CONSTRAINT pk_county_all_lookup PRIMARY KEY (st_code, co_code)) INHERITS (tiger.county_lookup);"
    ${PSQL} -c "VACUUM ANALYZE tiger_data.county_all;"
    ${PSQL} -c "INSERT INTO tiger_data.county_all_lookup(st_code, state, co_code, name) SELECT CAST(s.statefp as integer), s.abbrev, CAST(c.countyfp as integer), c.name FROM tiger_data.county_all As c INNER JOIN state_lookup As s ON s.statefp = c.statefp;"
    ${PSQL} -c "VACUUM ANALYZE tiger_data.county_all_lookup;"


texts:
I'm trying to connect my network to a points layer with a threshold of 4500 meters. I have about 11000 or so entries in my polyline and about 24 in my points layer.

I have already used v.distance, hence the threshold. However, I receive an error (image attached) any time I try the following:

v.net input=network points=cdw1c output=network_v1 operation=connect threshold=4500

Here is the error:

Error content is:

  Assertion failed: !p->taken[i], file split_q.e, line 69

Can anyone help me to understand what it means?

After a long consulation with a generous user from GRASS community, I still haven't came to a solution. I've tried on two computers with Win 7 as operating system, the following GRASS versions: 6.4.2 (through QGIS and without the plugin), 6.4.3 RCs, 7.0.0

For the first two, I have received the same error; for version 7, I didn't receive the error but the process have never finished. On the other hand, my helper managed to complete it in seconds using my files.

Can the problem have something to do with Windows 7, some background or system process, or so? Should I install any extension or plugin in order to run network analysis programs on Win 7?


texts:
I've just started with arcgis and I've got a number of different "maps" as individual layers. I've grouped them so it's easier to turn them on and off.
They are choropleth maps, so show counts of things.

Is it possible to assign the same grouping and colour schemes across the layers in one go?

I was thinking : select group layer, apply symbology, categories, join count, and then magic all the layers in the group have the same colour scheme.

Is that possible?

How?

OK Edit:

I've used the information below and learned the model builder and the apply symbology tool will apply any symbology from one layer to a group of layers (which I had to add individually, rather than adding the group layer).This is good. And works quickly!

While it applies the symbology, what I think I want is to set the same grouping across each layer. At the moment I'm counting number of houses per street

0 -   5 Red
6 -  10 Green
11 - 20 Yellow
21 - 30 Orange
31 - 40 Purple
41 - 50 Blue
51+ - Black

How can I apply a standard grouping and colour scheme across the multiple layers?

The maps currently show categories, should it show quantities?
Not every layer will have > 10 houses per street, should I build a template layer that uses all the grouping classes and apply that across the other layers (and ultimately the other map books)?


texts: Why is KML always selected as a standard service when publishing a ArcGIS Map Service?

texts: On publishing an ArcGIS 10.x Map Service, the capabilities list always has KML selected, and I always have to deselect it. Why is this always selected, and can I configure ArcGIS Server to keep it unselected when publishing services?


texts: Would a QGIS qt5 build use qt5's OpenGL rendering capabilities?

texts:
Related to this question, I am considering building QGIS with qt5 (caveats acknowledged) in order to improve rendering speed with qt5's OpenGL support. Before I attempt this I need to know whether QGIS (2.8/2.10) built this way would actually use qt5's OpenGL capabilities for improved rendering.

Can anyone confirm this and perhaps provide a guide as to the degree of rendering performance increase one may expect over a standard qt4 build (obviously context dependant, but order of magnitude indication would be nice)


texts: Plotting points using a list of lat & long coordinates

texts:
I have been given a list of latitude and longitude coordinates in a .csv file and asked to plot them on a map.  There are multiple other columns of data about each point included in the excel file.

I'm using QGIS to try to project them, but the procedures I have found both on this website, and on others have not given me the answer I need.  I have plenty of experience geocoding addresses on a list in an excel sheet, but have not had to plot coordinates before.  I've imported the list as a Delimited Text Layer, then saved it as a shapefile.  However the points do not project.

Is there a step I'm missing, or a plugin I need to acquire to get this done?


texts: Error in performing transect in ambur R?

texts:
I'm trying to filter a transect, window size is 5 and filter individual baselines is set to 0. I select the transect shapefile, then the inner baseline shapefile and it gives me this error.

Error in SpatialLinesDataFrame(shape.prep2, trandata) : row.names of data and
Lines IDs do not match.

I did everything correct up to this point, add fields to shapefiles, construct transects, stuck on filter transects. Anyone know how to fix this error or knows what it means?

The syntax for SpatialLinesDataFrame is

function (sl, data, match.ID = TRUE)
{
if (is.character(match.ID)) {
    row.names(data) = data[, match.ID[1]]
    match.ID = TRUE
}
if (match.ID) {
    Sl_IDs <- sapply(slot(sl, "lines"), function(x) slot(x,
        "ID"))
    data_IDs <- row.names(data)
    mtch <- match(Sl_IDs, data_IDs)
    if (any(is.na(mtch)))
        stop("row.names of data and Lines IDs do not match")
    if (length(unique(mtch)) != length(Sl_IDs))
        stop("row.names of data and Lines IDs do not match")
    data <- data[mtch, , drop = FALSE]
}
if (nrow(data) != length(sl@lines))
    stop("length of data.frame does not match number of Lines elements")
new("SpatialLinesDataFrame", sl, data = data)
}
<environment: namespace:sp>

not too familiar with code sorry. The error is happening because one of those if statements. Looking at the documentation on the R website I see that shape.prep2 is the object in SpatialLines class and trandata is an object in the data.frame class. "object of class data.frame; the number of rows in data should equal the number of Lines elements in sl." It must be something wrong with one of the shapefiles.


texts: In QGIS, add to existing dem tiles and expand corresponding contour coverage within a project

texts: In my first QGIS project I have dozens of digitized layers that have spread beyond the contour and dem tile area that I started with. Please recommend how to safely expand or replace the dem and contour layers with larger ones in the project.


texts:
I'm a total beginner with PostGIS, and trying to create a database with the geonames dataset.

Their table has latitude and longtitude coordinates:

  latitude          : latitude in decimal degrees (wgs84)
  longitude         : longitude in decimal degrees (wgs84)

Since the data is global and in WGS84, I've used  the  type:
geography(point,4326), and inserted it using ST_MakePoint(lat, long).

However, I get a lot of warnings:
NOTICE:  Coordinate values were coerced into range [-180 -90, 180 90] for GEOGRAPHY

Did I use the wrong type?


texts: Summarizing string values in QGIS 3.0 for a Point layer?

texts:
Is there a way to summarize string attributes using QGIS 3.0 without having to export to excel?

I have a point layer as an output of a spatial join. I would like to summarize the information so that it looks more like the "desired output" below using QGIS 3.0. (I recreated what I wanted out of QGIS for demonstrative purposes of what I am trying to do I'm QGIS).

I have tried multipart to singlepart (assume since points this does not work).

Tried the "list unique values tool" and either I don't know how to work the tool, or that is DEFINITELY not the answer.


texts:
I'm looking for a zigzag line symbol in QGIS. Is there perhaps an easy way to do this that I'm missing?
I have tried creating a marker line using a simple triangle marker (^) and adjusting the size of the marker and the marker placement interval until the traingles touched each other and appeared to make a nice zigzag line. This works for straight lines but around curves there are gaps between the triangles because the triangles aren't actually connected. Is there perhaps a way to join the markers together? Or another way to go about this? I would be very grateful for any suggestions!
(using QGIS 2.4.0)


texts: Calculating first quartile of a raster by zonal statistics?

texts:
I want to calculate the 1st quartile (25th percentile) of a raster (beneath polygons) using zonal statistics but it does not this option.

Zonal statistics has: MIN, MAX, MEAN, MAJORITY, RANGE, and does not have 1st quartile. How can I do that?


texts: Geoserver WMS Time Support - Absolute Interval

texts:
as the title suggests, I am currently working with a wms provided by geoserver 2.8.x, with the focus being put on a time-attribute. The official documentation is really lacking details on how to configure a layer with a date-field in geoserver, especially the "dimensions"-tab.

So the .shp i uploaded onto geoserver comes with two date-fields (begin and end), each containing, you guessed it, a date where a special action has begun and ended. I want my wms to only query features whose begin- and end-field are within a given period, e.g. by submitting &TIME=2014-03-22/2014-04-01 on getmap i want it to return features that lie completely within these boundaries.
However, geoserver returns any value where the given time period is somewhat included.

Does anyone know a way to work around that problem? Or maybe how to configure the layer itself in geoserver?


texts: I have a 6mb KML file I am using with the Leaflet API/Ominvore plugin. It is just a centerline of a long distance trail. What would be the best way to make it smaller? Reduce complexity etc.


texts: How to create desktop shortcut for QGIS?

texts: I installed QGIS version 2.8.1-Wien on Win 8.1 but there is no start or desktop icon to run QGIS. I can only execute using OSGeo4Shell which give me a command com line and I type in qgis to run. How can I create a desktop shortcut


texts: Importing CSV into QGIS with Python, Specific Parameters?

texts:
I would like to know where I can find a list of parameters that can be specified when importing a CSV into QGIS using python.  The documentation I've found only shows a few options.

QGIS Cookbook

I would like to specify the number of lines to skip, and the fields of the X and Y when there isn't a header.  How do I find the complete list of options?


texts:
I'm no expert with GIS development but I have created a tool that is used to interface with selected data in the Attribute Table that is currently open. Everything has been working fine since I wrote the tool but now all of a sudden my tool cannot find the CurrentLayer. My general question is "Is there a configuration setting that would cause this behaviour?"

Using ArcMap 10 on Windows 7. SQL backend, VS 2010 with C#.

In the code snippet below _map has zero Layers when I know they are there!

    private IMap _map;
    public MyEditorExtension()
    {
        _extensionStateService = Services.ExtensionStateService;
        _map = ArcMap.Document.FocusMap;
        _attributeWindow = (AddIn.FromID<AttributesWindow.AddinImpl>(ThisAddIn.IDs.AttributesWindow));
    }

The code below should iterate all Layers in the map and find the correct layer based on the featureClassName string passed in, but because the map has no Layers (code above) it can't return the Layer that I need to edit. I'm 99.99% sure that this is a configuration problem because it worked and then it stopped. I can't even get my production Compiled AddIn to work either, but my colleague has no problem, however I'm the only one developing the Tool.

    public static IFeatureLayer GetFeatureLayerByFeatureClassName(IMap map, string featureClassName)
    {
        IFeatureLayer result = null;
        for (var i = 0; i < map.LayerCount; i++)
        {
            var layer = map.Layer[i];
            var featureLayer = layer as IFeatureLayer;
            if (featureLayer == null) continue;

            var layerFc = featureLayer.FeatureClass;
            var layerDs = (IDataset)layerFc;

            if (layerDs.Name.Trim() == featureClassName.Trim())
            {
                //print the layer name in TOC..
                result = featureLayer;
                break;
            }
        }
        return result;
    }


texts: Reducing a box while keeping its center?

texts:
So I have to draw a second box inside the original box but reduced by a factor. I have the first coordinates of the box and I have to calculate the new coordinates that would draw the new box.

def getCenter(coords):
    coord1 = coords[0]
    coord2 = coords[1]

    x0 = coord1[0]
    y0 = coord1[1]

    x1 = coord2[0]
    y1 = coord2[1]

    middleX = (x0 + x1) / 2
    middleY = (y0 + y1) / 2

    return (middleX,  middleY)

def drawBox(t, coords):
    coord1 = coords[0]
    coord2 = coords[1]

    x0 = coord1[0]
    y0 = coord1[1]

    x1 = coord2[0]
    y1 = coord2[1]

    t.up()
    t.width(1)
    t.goto(coord1)
    t.setheading(0)  # East
    t.down()
    t.goto(x1, y0)
    t.right(90)
    t.goto(coord2)
    t.right(90)
    t.goto(x0, y1)
    t.right(90)
    t.goto(coord1)

def reduceSize(coords):
    middle = getCenter(coords)
    middleX = middle[0]
    middleY = middle[1]
    print(middle)

    coord1 = coords[0]
    coord2 = coords[1]

    x0 = coord1[0]
    y0 = coord1[1]

    x1 = coord2[0]
    y1 = coord2[1]

    newX0 = x0 - (middleX * FACTOR)
    newY0 = y0 - (middleY * FACTOR)
    newX1 = x1 + (middleX * FACTOR)
    newY1 = y1 + (middleY * FACTOR)

    return ((newX0, newY0), (newX1, newY1))

FACTOR = 0.10

So far I have this but it doesn't work for the 4 quadrants of a plot.

Edit:
With the currend implementation it works with the coords((-150, 150), (0, 0)) and a factor of 0.90 but if I give it coords in other quadrant the center of the second box changes


texts: Cannot find all tile services hosted in my AGOL via REST

texts:
I have already posted this problem to geonet, but I thought it might be good to post it here, as people on GIS SE might have already seen such a problem (although I could not find any similar question).

Basically, I have more than 100 tile layers hosted on my ArcGIS Online, and I need to list them through the REST URL. However, some of my tiles are not listed there. It seems like the interface is limiting the number of elements to 100 (regardless the publication date and alphabetic order).

Moreover, I can confirm:

all the services I am not finding in the REST are working properly.
the server hosting the tiles is the same (or at least it should, as it is in the item's detail of AGOL)

Just to calrify: if from the AGOL interface, I enter in one of these tile layers' details and then click on visualize (at the bottom right part of the detail page), it opens up the details of the service via the REST interface. If from this same page I come back to the home of REST, where all services should be listed, the same service is not present there.


texts: Defining ROI from shapefile in pyQGIS for use with google earth engine plug-in

texts:
I've just started experimenting with the QGIS Google Earth Engine plugin, but I am new to Python and struggling with what seems a relatively trivial task.
I would like to mimic the geometry drawing feature of the browser version of Earth Engine, and I think such a code snippet would be useful to many.
My first step was to create and draw new polygon in QGIS, save as a shapefile, and then try with the following code:
fn = "C:/path/to/shapefile.shp"
layer = QgsVectorLayer(fn, '', 'ogr')
feats = layer.getFeatures()
for feat in feats:
    geom = feat.geometry()
    polygondat= geom.asMultiPolygon()

roi = ee.Geometry.Polygon(polygondat)
#... now use roi in some fun gee code

This prints the following:
[[[<QgsPointXY: POINT(Xnum_1 Ynum_1)>,... <QgsPointXY: POINT(Xnum_n Ynum_n)>>]]]

... and unsurprisingly throws an error as input to the earth engine geometry definition.  I've tried approaches with 'fromPointXY' and 'fromPolylineXY' with no luck, although I could just be using bad syntax in my ignorance.


texts:
I am using OpenLayers 2.13.1 to display Zip-areas. I have a PostGIS database with multi-polygons, that I show on a map. The user can select/deselected them.

Some polygons have an area (I think they are called donut holes?) cut out, like the example below. When the large polygon is selected, it shows the hole in the middle even when the smaller polygon that is underneath it, is selected (as seen on image 2.). It looks like its a bug, because the polygon is transparent, but still I can't see the polygon underneath.

Also the small polygon in the middle can't be selected/deselected at all.

The first problem is the biggest, the other can be solved with a simple workaround.

Image1: Area 7722 and 7721 are selected, but 7721 looks like it isn't

Image2: Area 7722 is deselected. Now you can see that area 7721 is

UPDATE:

I created an basic JSFiddle example here. You cannot select the polygons, but the center of the colored polygon has another one over / under it that should be gray as well.

Hope someone can point me in the right direction...


texts: Buffers based on table relationship using ArcGIS Desktop?

texts:
I have what seems like a trivial problem, but I can't come up with a solution.

So I have two spatial tables (points) -

Buildings ~1000 records (addresses; with a common field identifying the builder)
Builders ~100 records (addresses, can be linked (join or relate) to 1st table based on builder ID).

My goal is to select all the buildings that are outside of 1 km zone (buffer) from their builder.

Is there a tool-based solution to that? I'm sure that it can be done using python, but I'm not familiar enough to come up with the code.
I'm using ArcGIS Desktop 10.3.


texts: NDVI reclassify with mean +1SD as threshold in Google Earth Engine

texts:
I'm trying to do ndvi reclassify in google earth engine. I know exactly how to do it in Arcgis but I'm very new with GEE and I got stuck with javascript. What I want to do is to use mean + standard deviation of ndvi as a threshold. Pixel value between -1 to mean +SD reclassify to 0 and pixel between mean +SD to 1 reclassify to 1. The two classes can be displayed into different colour as a same layer, but ideally I would like to have the "0" area with no colour so I can check how well the "1" area overlaps with vegetation.
Here's the code I have worked out so far:

var Blueskin =ee.Image('COPERNICUS/S2/20160721T223713_20160722T000334_T59GMK');
Map.centerObject(Blueskin, 10);
var Color = {bands:['B8','B4','B3'], max: 3000}
var Mask = ee.FeatureCollection(table);
Mask = Mask.geometry();
Map.centerObject(Mask);
var Blueskin = Blueskin.clip(Mask);
var ndvi = Blueskin.normalizedDifference(['B8', 'B4']);
var palette = ['blue', 'white', 'green'];
var ndviParams = {min: -1, max: 1, palette: palette};
Map.addLayer(ndvi, ndviParams, 'NDVI image');

var reducers = ee.Reducer.mean().combine({
  reducer2: ee.Reducer.stdDev(),
  sharedInputs: true
});

// Use the combined reducer to get the mean and SD of the image.
var stats = ndvi.reduceRegion({
  reducer: reducers,
  bestEffort: true,
});

// Display the dictionary of band means and SDs.
print(stats);

The aoi is the estuary in the image. I can get the value of mean and SD so I can actually manually set the threshold. But I don't know the code to set the threshold, reclassify pixels with the value lie between the same interval to a same integer (0 or 1) and give each area a colour.


texts:
I have a high resolution image, which I've broken down into tiles and named the tiles according to x,y co-ordinates and zoom levels. The image is not projected in any geospatial/gis (WGS84 ellipsoid, spherical mercator) ellipsoid.

If it is possible to use OpenLayers (version 2.x) to display the image in the browser, then how should I go about doing that?

Note that I also want to be able to overlay the image with other data, and OL allows me to do this easily on a map but not entirely sure if I can do this with an image.


texts: I created a very good symbol layout for a current map.  I will be creating similar maps just different regions and want to keep the symbols and markers that way they are.  The new data will have same attributes.  I know with ArcMap I could create just a layer file and the symbols and line colors would be saved.  In QGIS how do I do this.  I have looked on the board and saw you can save just one to a .qml file.  I want to save everything to the same file.  Is that possible even.  I tried grouping all my symbols together, but could not get it to save to a .qml.


texts: calculate difference between two interpolated grids using python

texts:
I'm trying to calculate difference between two interpolated grids. The grids are of different shapes, so cannot use np.subtract to calculate the difference. Any help on how to calculate the difference?

#create two interpolated grids from two set of xyz points and calculate difference
ixx = np.array(self.ix)
iyy = np.array(self.iy)
izz = np.array(self.iz)

iminx = min(ixx.flat)
imaxx = max(ixx.flat)
iminy = min(iyy.flat)
imaxy = max(iyy.flat)

# post surface - get extents
pxx = np.array(self.px)
pyy = np.array(self.py)
pzz = np.array(self.pz)

pminx = min(pxx.flat)
pmaxx = max(pxx.flat)
pminy = min(pyy.flat)
pmaxy = max(pyy.flat)

ixi = np.arange(iminx, imaxx, self.cellsize)
iyi = np.arange(iminy, imaxy, self.cellsize)
iX,iY = np.meshgrid(ixi,iyi)
iextent = (min(ixi), max(ixi), min(iyi), max(iyi))

pxi = np.arange(pminx, pmaxx, self.cellsize)
pyi = np.arange(pminy, pmaxy, self.cellsize)
pX,pY = np.meshgrid(pxi,pyi)
pextent = (min(pxi), max(pxi), min(pyi), max(pyi))

igrid = griddata((ixx, iyy), izz, (iX, iY), method='linear')
pgrid = griddata((pxx, pyy), pzz, (pX, pY), method='linear')

# would like to calculate difference, but this won't work, due to different shapes
changegrid = np.subtract(pgrid, igrid)

igrid shape is (181, 130)
pgrid shape is (184, 126)

Another important information is these are spatial grids with overlapping spatial extents. If we just broadcast/extend array shapes, would that be enough? I do not want to subtract the grids at wrong x/y locations. So basically like in GIS/Arcmap how we calculate raster difference at each cell at a particular xy cell location.

I'm not sure why this question is under hold and not considered as GIS. The grids are spatial grids with extent falling under GIS coordinate system.


texts: No Routes found for Locations soft-restricted network element

texts:
I am using OD Cost Matrix to find the least cost path for a number of locations and destinations.

There are few Origins that have this warning:
"Warning: Location "Location123" in "Origins" is on a soft-restricted network element."

And no Routes are found for these Origins.

I looked up in my map so make sure that the location is covered by network roads and it seems find on the map:

I tried to Google what does "soft-restricted network element" mean but could not really find anything online.

How do I solve this problem?


texts: Error with save path. Does not exist or is not supported

texts:
The final save raster file errors without message with Does not exist or is not support. But it does exist. If you can find anything wrong with this script I would be very greatfull. I had a vast array of errors with this bit of code. But much of it is already working in other parts of the script.

path ="C:/gisdata/india/imagery/"
parentdir ="off"

print "===== using Not isfile and listing folder names only ======="
for fname in os.listdir(path + parentdir):
  if not os.path.isfile(os.path.join(path + parentdir, fname)):
#   pathfolders.append(os.path.join(path, fname))
    dirname.append(fname)

print dirname

#Loop through image directories and list the tif raster files within each and
#process to reflectance values
for dir in dirname:
    print "Setting workspace"
    #Set workspace to current working directory
    env.workspace = path + parentdir + "/" + dir
    print "Workspace set!"

print "Creating local variables for input bands to Dark Pixel Subtraction"
#Local variable = Raster("LE71460402010085ASN00_B1.TIF")
Rt_B1dp = Raster(path + parentdir + "/" + dir + "/" + "Rt_b1.TIF")
Rt_B2dp = Raster(path + parentdir + "/" + dir + "/" + "Rt_b2.TIF")
Rt_B3dp = Raster(path + parentdir + "/" + dir + "/" + "Rt_b3.TIF")
Rt_B4dp = Raster(path + parentdir + "/" + dir + "/" + "Rt_b4.TIF")

sample = "C:/gisdata/india/imagery/sample"
bands = ("Rt_B1.TIF", "Rt_B2.TIF", "Rt_B3.TIF", "Rt_B4.TIF")
dpRasters = (Rt_B1dp, Rt_B2dp, Rt_B3dp, Rt_B4dp)
fieldNames = ("Rt_b1", "Rt_b2", "Rt_b3", "Rt_b4")

for dpRaster in dpRasters:
    for band in bands:
        for fieldName in fieldNames:
            #Sampling raster values using select points from Delhi_dp.shp
            print "Sampling raster values using selected points"
            arcpy.sa.Sample( band,"C:/gisdata/india/imagery/Delhi_dp.shp","sample","NEAREST")

            queryString = '"{0}" IS NOT NULL AND "{0}" > 0'.format(fieldName)
            #Use SearchCursor to query non NULL values for creating a new variable
            with arcpy.da.SearchCursor("sample", (fieldName),queryString) as cursor:
                for row in cursor:
                    dpValue = row[0]

                    print "Starting dpRaster subtraction"
                    if arcpy.Exists(dpRaster):
                        arcpy.Delete_management(dpRaster)
                    print "Subtracting " + str(dpValue) + " from " + str(dpRaster)
                    raster_dp = dpRaster - dpValue
                    print "Saving " + str(dpRaster) + " as image file"
                    raster_dp.save(path + parentdir + "/" + dir + "/" + dpRaster + ".TIF")
                    print str(dpRaster) + " has been created"


texts: Exporting layers from MXD to geodatabase maintaining same TOC structure?

texts:
I have imported several shapefiles from different sources and have set up an MXD project with groups and subgroups of layer. I want to export the layers to a geodatabase while maintaining the same structure as the TOC, sending each group of layers to a different dataset.

How could I?


texts:
I'm trying to use file output parameter in arcobjects geoprocessing tool which will be used as rest service on ArcGis server 10.1

IGPParameterEdit3 outParam = new GPParameterClass();
outParam.DataType = (IGPDataType)new DEFileTypeClass();
outParam.Value = (IGPValue)new DEFileClass();
outParam.Direction = esriGPParameterDirection.esriGPParameterDirectionOutput;
outParam.ParameterType = esriGPParameterType.esriGPParameterTypeDerived;
outParam.DisplayName = "ResultFile";
outParam.Name = "outResultFile";

What i do know is how to access file if that is input parameter:

IGPParameter inFileParam = (IGPParameter)paramvalues.Element[inFileParamID];
IDEFile deFile = inFileParam.Value as IDEFile;
IDataElement dataElement = deFile as IDataElement;
string filePath = dataElement.CatalogPath;

But what I'm trying to figure out is how to set output file. At the moment some things i have tried is:

IGPParameter outFileParam = (IGPParameter)paramvalues.Element[outFileParamID];
IDEFile deFile = outFileParam.Value as IDEFile;
IPersistStream pStream = deFile as IPersistStream;

ESRI.ArcGIS.esriSystem.FileStream fileStream = new ESRI.ArcGIS.esriSystem.FileStream();
fileStream.LoadFromFile(someFilePath);
IStream iStream = fileStream as IStream;
pStream.Load(iStream);

But last line throws COMException "Error HRESULT E_FAIL has been returned from a call to a COM component."


texts: Making part of lines layer transparent based on circle

texts: I have lines from a centerpoint going outwards at 10 degree intervals. I want to make the inner portion within the black circle of those lines transparent (25%). How can I do this?


texts: Smaller raster file natively supported by GeoServer

texts:
I have a 800MB GeoTIFF file generated from (ASTER GDEM v2 tiffs) depicting relief to be used as part of a basemap for an application I am working on. I intend GeoServer to serve this raster but it's way too big for my disk-quota allocation. When converted to JPEG200 it's about 15MB but can't use it since I can't install GDAL support on the host machine .

My question is: how else can I get a file < 100MB supported by GeoServer whilst retaining the relief depiction of my original GeoTiff.


texts: Google Earth Engine: Labeling polygons according to their attributes

texts: I went through different websites and GEE Guides to find a way to label polygons. But, I could not find a solution for it. Is there anyone who knows how to label polygons based on their attribute data (as in ArcGIS or QGIS)?


texts:
I found the gtfs2pgrouting project from pgRoutings' GitHub page. But it's three years old and doesn't work with the (current) full GTFS specification as it looks.

First of all I would like to know if there is a simple way to route GTFS data (preferably inside PostgreSQL). I think OpenTripPlanner is too full blown for me. And it has problems with huge/multiple graphs as it looks.

Otherwise I would write my own solution, connected to my already imported OSM/pgRouting instances.

So part two: Are there any fully functional import tools to get GTFS data into PostgreSQL 9.3+ that employ PostGIS 2+?


texts:
We changed a password on an Oracle schema, but apparently we have some rogue SDE clients - probably old ArcGIS map services on a dev machine somewhere - attempting to connect with the old password. SDE is causing Oracle to lock out the schema because of the rapidfire attempts to connect with the bad password.

Is there any way to detect where the bad login attempts are coming from? I have turned on the ArcSDE intercept log, but I don't see anything related to connect attempts there.

We're using SDE 9.3.1.