texts: Recommended way to use a setting file to initialize OpenLayers

texts:
I am new to OpenLayers. I am trying to figure out what is the best way to have a configuration file for the map I am creating.

For example, I am creating a map with the following style :

    var defaultStyle = new OpenLayers.Style(
                        {
                            pointRadius: 10,
                            fillColor: "yellow",
                            strokeColor:"#0500bd",
                            strokeWidth:2,
                            fillOpacity:0.4,
                            labelXOffset: "${xOffset}",
                            labelYOffset: "${yOffset}",
                            labelAlign: "${align}"
                        });

I would like to have the pointRadius and fillColor in a config file so I don't have to modify my source code but just the config file.

Shall I create a JS file; config.js with the following content:

      var configMap = {
         pointRadius: 10,
         fillColor: "yellow"
      }

And then, in my other JS that create the style, does:

    var defaultStyle = new OpenLayers.Style(
    {
        pointRadius: configMap.pointRadius,
            fillColor: configMap.fillColor,
        strokeColor:"#0500bd",
        strokeWidth:2,
        fillOpacity:0.4,
                    labelXOffset: "${xOffset}",
            labelYOffset: "${yOffset}",
        labelAlign: "${align}"
    });

What is recommended?


texts: Loading lots of points on leaflet through time slider.

texts:
I'm trying to add more than 13 thousand point using Time Slider plugin, but as you can imagine the loading process of all these points puts a heavy burden on the map. Even the slider start to get sluggish. Is there any way to increase the application's performance? I played  bit with the Marker Cluster plugin and it manage to handle all the points very well, but this was without the time slider. I tried to combine both but could not do it. Is it possible to combine both plugins? Or is there a better combination involving other plugins?

best,


texts: OpenLayers only keeps the first drawn feature

texts:
I would like to save all my drawings on a vector layer stored in a PostgreSQL database and shared by GeoServer for a web mapping application.

I wonder why this piece of OpenLayers/JS code (as found here : http://codepen.io/barbalex/pen/fBpyb) doesn't work. Only the first drawn feature is kept.

Why?

Is it a bug with OpenLayers?


texts: Batch editing of circle diameters using attribute table in ArcGIS Desktop?

texts: I need to edit the diameter of a bunch of circles to reflect their correct diameters, and I wondered if there was a way of drawing from the attribute table to make all the circles equal their correct diameters rather than redrawing all the circles with the correct diameters by hand?


texts: Multiplying two iterating rasters in ModelBuilder or Python?

texts:
I have several Landsat scenes that I wish to remove clouds from. I can do this by multiplying the Landsat layers by their associated cloud layers in the Raster Calculator. How can I set up iterations for both raster, followed by multiplication in Raster Calculator, in ModelBuilder?

The Landsat layer and cloud layer has the same parameters but are stored in separate folders. Ideally, I seek a ModelBuilder solution, but I am also open to python.

I have tried using a submodel, which lets me iterate through the two sets of raster folders, but the model fails at the Raster Calculator. The equation Raster 1 x Raster 2 works outside ModelBuilder. Does anyone know where I have gone wrong?


texts: Access geoprocessing service via REST, get JSON back

texts:
I have an ArcGIS Server 10.0 geoprocessing task with a REST endpoint something like this:
http://server/arcgis/rest/services/MyService/GPServer/MyTask

The task takes a parameter. I was hoping I could pass it like so:
http://server/arcgis/rest/services/MyService/GPServer/MyTask?param=Value&f=pjson

... and get the result back as JSON. But it doesn't seem to respond to the parameter - all I get is static metadata about the task. There is the option of submitting a job, but my JavaScript will have no idea when it finishes, so I would have to keep checking the job till it's done and then get the results via individual URLs:

http://server/arcgis/rest/services/MyService/GPServer/MyTask/jobs/a72h72hd6d6g26gdh7/results/someParam?f=pjson

Is there really no way to just get the results as JSON? And pass the parameter in the first place?


texts: Calculating file name with field calculator

texts:
Is it possible to calculate the file name in an attribute table?

My shapefile has a name like "cities.shp". In this shapefile I want to have a attribute "cities").


texts: How to show QGIS maps in android application?

texts:
I have created map using QGIS software and trying to display data-set records on that map. Also developed android application where user want to see that map.
But I am stuck in how to show this QGIS map on android application.

I am working on crime profiling and showing on map to user phone.


texts: How to Create SCS (now NRCS) Curve Number Grid using GIS with European map (ESDAC dataset)?

texts: After finishing the watershed delineation, and after creating the HMS schematic in Arcmap, i was trying to include the land use map and  soil map in the process, by using the CNgrid preparation by using Analysis tools-Overlay-Union command. But in order to do that the attribute table for both layers must contain specific column headings. Unfortunately I do have only the names of the each feature in the soil layer, so I have to populate the other features on my own. In the soil layer there must be 4 columns:"PctA”, “PctB”, “PctC”, and “PctD”.  These columns contain the percentage of each hydrologic soil group (A, B, C, and D) in the soil layer polygons, which I populated with our geologist. But also there are two other columns there are very important in order to do the join of the features of the soil layer, these columns are MUSYM and MUKEY, where  MUSYM - is a map unit symbol that should only be used when working with the data from 1 survey. I do not have these data and i dont know to populate them on my own, but without them i can not join the layers. i searched at dataset ESDAC but i dont know which of them to download and what method to use.


texts: Specify data type when writing a GeoTIFF

texts:
I'm using org.geotools.gce.geotiff.GeoTiffFormat to write a GeoTIFF. I would like to use Float32 instead the default Float64. I've tried:

    OutputStream output = ...
    Hints hints = new Hints(Hints.SAMPLE_DIMENSION_TYPE, SampleDimensionType.REAL_32BITS);
    GridCoverageWriter targetWriter = new GeoTiffFormat().getWriter(output, hints);

But it still writes as Float64 (as verified with gdalinfo). Any ideas?


texts: Excluding layers from GeoExt 2 Legend Panel

texts:
I'm trying to exclude some legend in GeoExt.panel.Legend.

I've created a vector layer named 'City Extent' via a function but I want the legend to be removed in the list of legends in my legend panel.

The 'filter' option doesn't seem to work

Here's my vector layer for City Extent

city_ext = new OpenLayers.Layer.Vector("City Extent", {
                strategies: [new OpenLayers.Strategy.Fixed()],
                protocol: new OpenLayers.Protocol.WFS({
                    url: geoserverURL,
                    version: "1.1.0",
                    featureType: "gadm2",
                    featureNS: "http://www.adburbaninfo.org",
                    srsName: "EPSG:900913"
                }),
                filter: new OpenLayers.Filter.Logical({
                    type: OpenLayers.Filter.Logical.AND,
                    filters: [
                    new OpenLayers.Filter.Comparison({
                        type: OpenLayers.Filter.Comparison.EQUAL_TO,
                        property: "id_2", // city
                        value: params_2
                    }),
                    new OpenLayers.Filter.Comparison({
                        type: OpenLayers.Filter.Comparison.EQUAL_TO,
                        property: "id_1", // province
                        value: params_1
                    }),
                    new OpenLayers.Filter.Comparison({
                        type: OpenLayers.Filter.Comparison.EQUAL_TO,
                        property: "id_0", // country
                        value: params_0
                    })
                ]
                }),
                eventListeners: {
                    'loadend': function(evt) {
                        map.zoomToExtent(city_ext.getDataExtent());
                    }
                }
            });

Here's my GeoExt Legend Panel

var legendPanel = Ext.create('GeoExt.panel.Legend', {
            defaults: {
                labelCls: 'mylabel',
                style: 'padding:5px'
            },
            bodyStyle: 'padding:5px',
            width: 250,
            autoScroll: true,
            layerStore: mapPanel.layers,
            // dynamic: false,
            loader: {
                filter: function(record) {
                    return record.get("layer").name.indexOf("City Extent") == -1
                }
            }
        });

Could someone teach me the proper way to exclude/remove the legend?


texts:
Which software is better to do atmospheric and radiometric corrections in Landsat images?

After the corrections, I hope them to be used in supervised classification according to land use. So is it ok to correct the image by the MTL file of the image in ENVI? Or do I need to correct the 3 RGB bands seperately for this task?


texts:
Let's say I have a layer made up of pre-selected features from another layer:

selects = arcpy.SelectLayerByAttribute_management (S_FIRM_Pan, "NEW_SELECTION", "DFIRM_ID = '{}'".format(fipC))

If I get a count of the number of these selected features, I get says, 199.

getCount = arcpy.GetCount_management(selects)
print(str(getCount))
>>>199

Checking in arcMap, I can see that this is the correct number. I then try to use arcpy.da.SearchCursor to get all of the values for a certain field for selects:

raw_firm_pan = [row[0] for row in arcpy.da.SearchCursor(selects, "FIRM_PAN")]

If I use print(len(raw_firm_pan)) it only gives me 198 items. When I print it out and compare it to the actual rows for selects in arcMap, it appears that the last row is being cut off! How can I fix this?


texts: How to keep layer style/symbology when saving as KML?

texts: When saving as a kml, is it possible to keep the style of the QGIS layer?  I have tried a couple plugins that have been mentioned on here: qgis2kml ('format not yet supported') and qgis2google2 (crashed QGIS upon activating). Is there anything else?


texts: Understanding Error network dataset has no edge elements It may be unbuilt?

texts:
I'm trying to alter an Integrated Transport Network (ITN) to remove all of the roads that are within a flood zone (a separate shapefile).

The ITN Network Dataset I'm working with was created by Clipping the original feature classes from within the downloaded ITN .gdb and then rebuilding the Network Dataset so that it is restricted to the county boundary I'm working within. The rebuilt Network Dataset works fine for routing.

However, I need to remove roads that are within a flood zone so I can compare travel times. When I use the Erase tool to delete the parts of the individual feature classes of the clipped ITN network that are within the flood zone and then rebuild the Network Dataset from the new feature classes I get the following error when I try to run Closest Facility or Service Area routing:

Error: The network dataset has no edge elements. It may be unbuilt.

I made the Network Dataset exactly the same way that I made the clipped dataset, just using the new feature classes that had the flood zone roads erased. The elevation, turn, directions and so on settings are all the same and the Network itself looks fine (the correct bits have been erased).

Have I done something wrong? Is there a better way to remove the roads that are within the flood zone than using the Erase tool?

I followed @Mapperz's suggestion below to build the network dataset again and I end up with the same output (see below image) and get a long list of build warnings/errors:

SourceName: ITN_RoadLinks_CLIP_FLOOD_1, ObjectID: 1, Invalid hierarchical value.

SourceName: ITN_NetworkDataset_Junctions_CLIP_FLOOD_1, ObjectID: 36332, Standalone user-defined junction is detected.

SourceName: ITN_TurnRestrictions_CLIP_FLOOD_1, ObjectID: 1751, Cannot find edge element corresponding to turn identifier  1.

The original Network Dataset and the Clipped Network Dataset seem fine, with both reporting to be Built. But then when I use the Erase tool to create new feature classes and then create a new Network Dataset, it never builds and looks weird:


texts: Downloading Points of Interest from Google Maps?

texts:
Is there any legal way to search and download google maps POI information, based on user location and search radius (50 km) + search criteria (eg restaurant) ?

Output would ideally be KML or CSV.


texts:
I performed radiometric correction analysis on Landsat 8 image ( format: integer unsigned 16 bit) in ERDAS IMAGIN 2014.

However, the output image would be created if I could select format: float 32 bit.

It was not created in unsigned 16 bit.

My question is, was this float 32 bit output correct or not?


texts:
I have looked at the first answer and the answer to this question and I'm on QGIS 2.8.2 (Ubuntu 14.04) but the Web menu doesn't appear even if OpenLayers plugin is shown to be in the Installed plugins.

Where can get the OpenLayers plugin to show off ?


texts: Meaning of numPixels parameter for unsupervised clustering in Google Earth Engine

texts:
I did an unsupervised classification with google earth engine. numPixels parameter change the classification result hugely. I am very curious about how this parameter work and what does it mean and what is the optimal values for best results. My clustering code below here.
var training=clipdata3.sample({
  region:geometry,
  numPixels:101,
  scale:10
});
var grupo= ee.Clusterer.wekaKMeans(3).train(training); //K-Means Method with 3 classes
var classified= clipdata3.cluster(grupo);


texts: Help with technique to create contours from ASTER DEM in Phnom Penh, Cambodia

texts:
I have the ASTER version 2 DEM of Cambodia. I want to use it to create 1m contours of Phnom Penh, Cambodia. If anyone has this already please let me know. Now, onto the problems:

I am using ArcGIS so do not have advanced raster processing tools at my disposal such as ER Mapper
The ASTER DEM seems very noisy in Phnom Penh, using the contour tool in ArcGIS simply creates lots of contours around what looks like a field of knolls
Using the focal statistics tool, even with filter sizes up to 30m, does not appear to fix the problem.
I am comparing the output to 10 and 20m contour lines digitised from a topo map. These lines sweep across the city like contour lines should. They are in no way mimicked by the contour lines I have so far been able to generate from the ASTER DEM.

Does anyone have any advice for generating contours from the ASTER raster?


texts: Selecting only parallel lines along center line?

texts:
I have read through how to find close and parallel lines but I am still looking for a simpler solution. Found several ArcGIS forums post on the topic but no answers there either.

Does a tool exist in ET GeoWizards, ArcToolbox, or QGIS that will select only parallel lines (with a given tolerance)?


texts: Browser crash when using marker cluster in leaflet

texts:
When using the following portion of code, the clusters and some markers are shown but if I zoom or click anywhere the browser crashes.

var cam = L.ExtraMarkers.icon({
             icon: 'fa-video-camera',
             markerColor: 'red',
             shape: 'square',
             prefix: 'fa'
          });

var cameras = omnivore.kml('json/cameras.kml')
            .on('ready', function(layer) {
               this.eachLayer(function(marker) {
                 marker.setIcon(cam);
                 var markers = L.markerClusterGroup();
                 var popupCam = marker.toGeoJSON().properties.name;
                 marker.bindPopup(popupCam);
                 markers.addLayer(cameras);
                 map.addLayer(markers);
               })
             });

Without using the clustering, all markers (~300) are shown pretty well so I guess it should not be performance related.


texts: How to estimate time spent within a location using GPS/timestamp data in PostGIS?

texts:
I am new to PostGIS, and I am trying to estimate time spent in a location based off of interpolating time between GPS points (I am using PostgreSQL 12.2 and PostGIS 3.0).

I have GPS points stored with timestamps (With data columns: Latitude, Longitude, Date, and Time) of a path taken. The time intervals in between each point are not the same and they vary. I also have a polygon of interest through which the GPS path cuts (multiple times entering and leaving)- see image below for a conceptual example of the data with which I am working, with blue dots being the points and the red line being the hypothetical path entering and exiting the black polygon.

I'm interested in estimating the time spent in the polygon based off of the time represented by the segments inside the polygon (when a segment is partially in/partially out, I am aiming to calculate the time represented by the fraction that is in- I am accepting the assumption that a fraction of a segment corresponds to an equal fraction of time in between the points).

I know that ST_Intersection would simply return the section of the overall line within the polygon, but I don't believe I'd be able to just take the proportion of the length inside compared to the total length, because different segments represent different lengths of time (Please correct me if I would actually be able to use ST_Intersection in this case). It is straightforward when both the beginning and end points in a segment are within the polygon but less so when one is inside and the other is outside. My ideal output would be a sum of the lengths of time represented by segment/segment fractions within the polygon.

I have looked at:

Converting GPS points with timestamp to polylines in PostgreSQL?
How to interpolate GPS Positions in PostGIS
Storing GPS tracks with timestamps in PostGIS
Store GPS trace with timestamps, then interpolate location by time

but most of these seem to look at estimating location based off of time (e.g. predicting where location will be in between existing data points), whereas I sort of want to do the opposite.

Would it better to do calculations for individual pairs of points given the uneven time between each pair (i.e. compare fractions of segments via ST_LineLocatePoint and/or ST_LineSubstring) or would it be better to do as one linestring/trajectory?

Ideally, I am hoping to do this all within the database/PostgreSQL/PostGIS.


texts: How do I calculate the coordinates of the median point of a polyline?

texts:
I have a shapefile with some polylines. For each one of these polylines I'd like to calculate the median of the polyline (i.e. the point that's equidistant from each end) and then, calculate its coordinates.

I don't know if this is feasible.
The only tools I've found regarding this matter are Mean Center and Median Center, however this is not what I'm looking for.

Any help would be much appreciated.


texts:
I have a DEM in the form of a raster that I would like to update in certain locations. Ideally, I'd like to be able to create a separate polygon or point shapefile with the values, in the correct spatial location, to update the original DEM with. Is there a way to do this?

I have a spatial analyst extension and ET Geowizards. I do not have 3D Analyst.

I need to use ArcGIS for this.


texts: Extracting raster value to points in ArcGIS

texts: I extracted my slope map to the points. I had to add the slope values to the points through extracted values to point but in the result's table of attribute I can't see any slope value. Can you help me with that?


texts:
I'm using the Count points in polygon tool to estimate points and make some time-series in QGIS 2.18.4. I edit shapefile of points (just delete rest) to make a certain date range.
As example:
1) 1995-2015

But in order to see the difference I must copy (Right click on layer - Styles - Copy) the styles from 1) to

2) 1975-1995.

Attribute table number of points has the same name. What can cause such a result? I tried to use different shapefiles from GADM and Natural Earth


texts: How to calculate degrees per pixel in a given image? Or .nth of degree?

texts:
I am trying to make very rough estimates of Lat/Long coordinates in some random parking lot images. If I visit any such random parking lot on a Google map, I can click some driving area and Google gives me a location. For example.

Here is a parking lot image from Florida. (Hosted on wikimedia commons.)

The image is 1280 pixel width and I can find some pretty close coordinates for the corners of the building on the google map. Here is the corresponding map view served up at Google.

The left side of the Big Lots is very close to 30.180806, -82.670137. And the right side of the store is very close to 30.180752, -82.669594

So there is a LAT difference of 0.000054 and a LONG difference of 0.000543.

Let's say I measure the pixels of the Big Lots CarPark in the image and find it uses 800 pixels. With that, I can approximate the width of each pixel.

0.000543 / 800 = 0.000000679 per pixel.

My First Question (not so important) is: what is that 0.000000679 measurement called? (One 0.000000679 th of a degree per pixel?)

Anyway, to the real question.

If I have some other image, for which I cannot reliably find the location on map or GIS, how can I make a similar estimate, given close guess of the the locations LAT and LONG?

Example:

Here is another parking lot image for which I cannot know the exact location. I can make a guess that the store front is something like 100 feet wide, by counting the parking spaces. I might even guess the image is 120 feet wide, just for the sake of making calculations.

Then let's say I assign an arbitrary LAT / LONG to the pixel at the left side of the image, at ground level with the building foundation.

LAT  34.065154
LONG -118.468884

Note: I am also arbitrarily assuming the camera is facing directly north.

I can find on the internet that at this latitude, there are 363922 feet in a degree of longitude.

So, how can I calculate the 0.0000XXX per pixel?

I am given:

rough estimate of ground distance, 120 feet
estimated 363922 feet per degree latitude
left most pixel location:  34.065154, -118.468884
image pixel width 1280

I should be able to estimate the LAT / LONG on the other side of the image right?

Here is my uneducated guess:

120 feet / 1280 pixels = 0.09375 ft/pixel
0.09375 ft/pixel / 363922 ft/deg/lat = 0.000000258 nth degree per pixel?

If so, do I calculate the right side of the image like this?

-118.468884 + (0.000000258 * 1280) =
-118.468884 + 0.000329741 =
-118.468554259

LAT  34.065154
LONG -118.468554

Is the same working for the latitude if I estimate the distance to the camera by counting the cars for instance?


texts: Creation of a polynomial function from old and new coordinate system pairs

texts:
I have data points in the old Roussilhe Stereographic coordinate system of Suriname (u,v), I have the same data points also in WGS84 UTM21N coordinates (x,y). I want to create a polynomial approximation conversion function f(u,v)=(x,y) for use in excel for people who do not have access to geodetic software.
I have access to the tools of proj.org to generate as many points as needed. How do I create a polynomial function from these point pairs (u,v,x,y) ? 5 meter accuracy would be sufficient for this purpose.


texts: QGIS: why does simplesvg export plugin generate poor quality output?

texts: I am exporting a QGIS project to SVG via the simplesvg plugin. However the output is of really poor quality: all coordinates appear to have snapped to a rather coarse grid (while original data, as shown in QGIS main window, have not). Why is this happening?


texts: Batch Create Circle Sections based on center point, radius and azimuth range

texts: I have over 1000 points with radius and min/max azimuth. I need to create a circle section based on these three columns for each point (buffer each point by radius but only keep the section between min and max azimuth). Anyone has a good way to do this?


texts:
I try with QGIS to insert new features in a multi-table Oracle Spatial view that is updatable by adding an "instead of" trigger. The view is "key-preserved" and I don't have any problem when I insert a new feature directly in SQL.

But when I try to do this in QGIS (QGIS 2.18.24 as QGIS 3.2) I get this error :

Oracle error while adding features: Oracle error: Could not retrieve feature id -2
SQL: Error: SELECT "OBJECTID" FROM "MY_SCHEMA"."MY_VIEW" WHERE ROWID=:f

Despite this error, the feature is inserted in my database anyway. The problem is, if I save layer edits another time, my feature would be recorded 2 times.

Is this a limitation of the Oracle Spatial QGIS connector or I can do something to have a workaround solution?


texts: Count values from another layer based on value of current feature

texts:
I have two layers. One is a shapefile of U.S. states where each state is labeled with the attribute STATE_NAME, and the other is a CSV data layer named rid with many points that are labeled with a State attribute. In this CSV data layer, there are many points for each state.

I want to create a column in the shapefile that contains the count of points per state.

Using the field calculator in the shapefile, I've attempted

aggregate(layer:= 'rid',
          aggregate:= 'count',
          expression:='State',
          filter:=attribute($currentfeature, 'STATE_NAME'))

but this results in 0 for each state. If I attempt without the filter, I get the total number of points regardless of state, resulting in the same number for each state.

Is there a filter that will allow me to count the number of points per state?


texts: Non finite transformation detected - lidR spTransform

texts:
I would like to re project the following laz file using lidR and spTransform from sp.
The EPSG code is not defined but it seems to me that the CRS string is complete:
tmp_laz <- readLAS("./USGS_LPC_CA_FEMA_R9_Russian_2017_10SDJ993053_LAS_2018.laz")
epsg(tmp_laz)
#>[1] 0
tmp_laz@proj4string@projargs
#> [1] "+proj=lcc +lat_0=37.6666666666667 +lon_0=-122 +lat_1=39.8333333333333 +lat_2=38.3333333333333 +x_0=2000000.0001016 +y_0=500000.0001016 +ellps=GRS80 +units=us-ft +vunits=us-ft +no_defs"

I tried to project to the global metric datum EPSG 3857 and got the following error non finite transformation detected:
tmp_laz_2 <- spTransform(tmp_laz, sp::CRS("+init=epsg:3857") )
#> non finite transformation detected:
#>   X       Y
#> [1,] 6275252 2094773 Inf Inf
#> [2,] 6275249 2094771 Inf Inf
#> [3,] 6275246 2094769 Inf Inf
#> [4,] 6275243 2094768 Inf Inf
#> [5,] 6275251 2094768 Inf Inf
#> [6,] 6275251 2094775 Inf Inf
#> [7,] 6275248 2094774 Inf Inf
#> [8,] 6275245 2094772 Inf Inf
#> [9,] 6275242 2094770 Inf Inf
#> [10,] 6275239 2094768 Inf Inf
#> ... goes on till reaching print limit
#> Erreur : C stack usage  100062272 is too close to the limit
#> De plus : Warning messages:
#>   1: In showSRID(uprojargs, format = "PROJ", multiline = "NO") :
#>   Discarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs
#> 2: In showSRID(uprojargs, format = "PROJ", multiline = "NO") :
#>   Discarded datum WGS_1984 in CRS definition
#> 3: In sp::proj4string(x) : CRS object has comment, which is lost in output
#> 4: In sp::spTransform(spts, CRSobj) :
#>   12346945 projected point(s) not finite

I tried another (local) projection, EPSG 6417 and got the same error:
tmp_laz_6417 <- spTransform(tmp_laz, sp::CRS("+init=epsg:6417") )

#> non finite transformation detected:
#>   X       Y
#> [1,] 6275252 2094773 Inf Inf
#> [2,] 6275249 2094771 Inf Inf
#> [3,] 6275246 2094769 Inf Inf
#> [4,] 6275243 2094768 Inf Inf
#> [5,] 6275251 2094768 Inf Inf
#> [6,] 6275251 2094775 Inf Inf
#> [7,] 6275248 2094774 Inf Inf
#> [8,] 6275245 2094772 Inf Inf
#> [9,] 6275242 2094770 Inf Inf
#> [10,] 6275239 2094768 Inf Inf
#> ... goes on till reaching print limit

#> Erreur : C stack usage  100062272 is too close to the limit
#> De plus : Warning messages:
#>   1: In showSRID(uprojargs, format = "PROJ", multiline = "NO") :
#>   Discarded datum NAD83_National_Spatial_Reference_System_2011 in CRS definition
#> 2: In sp::proj4string(x) : CRS object has comment, which is lost in output
#> 3: In sp::spTransform(spts, CRSobj) :
#>   12346945 projected point(s) not finite


texts: Obtaining average value around polygon in raster for different cross sections?

texts:
I have elevation data in a raster, river (blue line) and the flooded area (black polygon) as shown in the image below:

I want to find the water depth in the flooded area (inside the black polygon). Outside of the polygon, water depth is obviously zero. Inside the polygon, the land elevation is lower than outside (At least it should be)

What I tried to do and did NOT work was this: I assigned nodata value to the points inside the polygon. Then I used Eucilidean Allocation to find the closest points just outside the polygon and assign it to nodata points. Then, I Subtracted the original land elevation raster from this new raster. However, I get big negative values because after all, it was not a good method at all.

Now, I am thinking: I need to follow the water stream line. At each point, draw a perpendicular line to the river and let it intersect minimum two times with the polygon. Find the average value of cells inside a circle (with those intersections as center). Find the mean of these averages. apply this value too all points in that cross section (perpendicular line to the river). Now, is there any way on earth I can do this? Hopefully I don't need a for loop for this? Do you guys have any other suggestions, perhaps easier ways?

I am using ArcGIS 10.1 with Spatial Analyst and 3D Analyst extensions.


texts: Change projection of list of rasters in R

texts:
I am handed down some RData that contain a list of of raster stacks. The file structure is pretty challenging to manipulate. After doing load(), the data looks like this:
[[1]]
class      : RasterStack
filename   : raster1
dimensions : 82, 101, 8282, 1  (nrow, ncol, ncell, nlayers)
resolution : 2.008502, 1.991304  (x, y)
extent     : 488240.2, 488443, 4950773, 4950937  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=18 +datum=NAD83 +units=m +no_defs
names      :    Y2020
min values : 24.14956
max values : 24.14958

[[2]]
class      : RasterStack
filename   : raster2
dimensions : 435, 157, 68295, 1  (nrow, ncol, ncell, nlayers)
resolution : 1.997977, 1.999499  (x, y)
extent     : 488481.3, 488795, 4950998, 4951868  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=18 +datum=NAD83 +units=m +no_defs
names      :    Y2020
min values : 20.94342
max values : 22.16785

...

[[n]]
class      : RasterStack
filename   : raster_n
dimensions : 225, 261, 58725, 1  (nrow, ncol, ncell, nlayers)
resolution : 1.998339, 1.998643  (x, y)
extent     : 484329.9, 484851.5, 4950577, 4951027  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=18 +datum=NAD83 +units=m +no_defs
names      :    Y2020
min values : 20.62384
max values : 20.62384

I would like to reproject the list of rasters using spTransform() through dplyr. Here is my code:
utm18 <- "+proj=utm +zone=18 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"
new_crs_myrasterlist <- sapply(myrasterlist, function(x) spTransform(myrasterlist[[x]], CRSobj=utm18))

I also tried,
new_crs_myrasterlist <- sapply(myrasterlist, function(x) spTransform(myrasterlist[[x]], CRS(utm18)))

Both code gave me this error:
error in evaluating the argument 'x' in selecting a method for function 'spTransform': invalid subscript type 'S4'

What am I missing? Is there a better alternative?


texts: LRS plugin dropdown selections not working

texts:
I have just loaded the LRS plugin but when I try to access the layer with measures drop down selections on the measure tab they are not working.  Calibration drop downs work and so to do the layer directly below layer with measures

I have reloaded the plug in, rebooted my computer, not sure what else to try.

calibration tab works fine but the locate, events and measures top option all have the same issue.


texts: Why won't my csv file imported as a vector layer have attribute table headers?

texts:
I have a .csv file that has unique identifiers that I'd like to join with some legislative district shapefiles in QGIS. I've done this before, and the joins all worked fine -- and technically it still works if I memorize the order of the columns, though when I upload the .csv using "Add Vector Layer" like normal -- my .csv gets imported properly, except my headers are all in as the first row and the headers are all "field_1," "field_2," etc.

Converting the file to .dbf in OpenOffice and I attempted converting the line breaks by using this in Terminal:

tr '\15' '\n' < infile.csv > outfile.csv

Still no luck. Any other ideas to import my .csv and retaining all my column headers?


texts: error creating Geoserver layer from Multipolygon PostGIS table

texts:
I imported the table from a multipolygon shapefile using PostGIS Shapefile and DBF Loader 2.2.    I tried to use the 'generate simple geomtries ...' option in the tool, but got back an error : 'We have a Multipolygon with 3 parts, can't use -S switch! Shapefile import failed.', so I imported the shapefile as a multipolygon.

It seemed to import ok, and a preview of the table look ok as well.  But when I try to create a published layer in Geoserver from this table, nothing renders in my Javascript app, and a preview of the layer shows that the multipolygons are missing.

What can I do here?

below is the schema from the PostGIS table:


texts: Overlapping Polygon becoming hole when converting geojson to shapefile

texts:
I'm trying to convert a geojson MultiPolygon which contains two overlapping polygons to a shapefile, however I noticed that the resulting shape has a smaller area. Upon inspection it's because the overlapping polygon is becoming a hole in the shapefile. Do I need to change a setting or do something different to preserve the overlapping polygon? I'm using GDAL 1.10.1 and the following code:

def convert():
    data = {
        "type": "FeatureCollection",
        "features": [
            {
                "geometry": {
                    "type": "MultiPolygon", "coordinates": [
                    [[[-98.7, 49.6], [-98.7, 49.7], [-98.8, 49.7], [-98.8, 49.6], [-98.7, 49.6]]],
                    [[[-98.74, 49.64], [-98.74, 49.66], [-98.76, 49.66], [-98.76, 49.64], [-98.74, 49.64]]]
                    ]
                },
                "type": "Feature",
                "properties": {}
            }
        ]
    }

    geometry = []
    if (data['type'] == 'FeatureCollection'):
        for feature in data['features']:
            geom = ogr.CreateGeometryFromJson(json.dumps(feature['geometry']))
            print 'Geom Area', geom.GetArea()
            geometry.append({
                'geometry': geom})

    driver = ogr.GetDriverByName("ESRI Shapefile")
    srs = osr.SpatialReference()
    res = srs.ImportFromEPSG(4326)
    print 'Import Res', res
    data_source = driver.CreateDataSource('out.shp')

    layer = data_source.CreateLayer('shape-layer', srs, ogr.wkbMultiPolygon)
    featureDefn = layer.GetLayerDefn()
    layerDefn = layer.GetLayerDefn()
    for data in geometry:
        feature = ogr.Feature(featureDefn)
        feature.SetGeometry(data['geometry'])
        layer.CreateFeature(feature)
    data_source.Destroy()

Viewing the geojson and resulting shapefile in QGIS, this is what I see:

Initial Geojson

Resulting Shapefile


texts: Using the C++ GDAL API how do I set the GDAL_DATA path?

texts: If I want to change the directory of the GDAL_DATA path or point to a different version of a csv file for EPSG numbers how do I do that using the C++ API?


texts:
I have a dataset of manholes and pipes; I want to snap the lines start and end nodes to the manholes if they are within a specified distance. Been trying to find away to do this using PostGIS or QGIS because I don't have access to the Arc Desktop tools to do it. I found a post of a similar situation but I cannot get the function to work, I don't know what information is suppose to go where in the expression.

    SELECT f.gid as gid,
        ST_Snap(f.Geometry, g.Geometry, 2) as geom
    FROM
        pipe as f,
        (SELECT ST_Collect(Geometry) as Geometry
         FROM mh) as g

I thought it was suppose to work like this;

            SELECT f.id as gid,
                ST_Snap(f.geom, g.geom, 2) as geom
            FROM
                smains as f,
                (SELECT ST_Collect(geom) as Geometry
                 FROM samh) as g

I get column g.geom does not exist, but the column samh.geom does and to my understanding g = samh. When I change it to samh.geom I get missing FROM-clause entry for table "samh" which sounds to me that the variable samh is not defined, but it isn't a variable it is a table.


texts:
Whenever I do a presentation about anything related to GIS, I consider the level of knowledge of GIS of the audience. The worst times are definitely when I need to explain GIS to a group that has almost no knowledge about it. I said virtually, because by know everyone has interacted with some type of electronic map in one way or another (elephant in the room is Google Maps), so I usually start from there and explain everything that is necessary to create Google Maps. I have found that works quite well.

I know there are several formal definitions of GIS.

The questions are:

Which formal definition of GIS do you use?
How do you explain it to somebody that doesnt know much about the field in elevator pitch style? (30 seconds)
How do you explain it when you do have more time? (10 min)


texts: I have a georeferenced image and I need to use it as a basemap map, instead of OpenStreetMap. is there a way to do it with simple leaflet?


texts: Editing capabilities document in GeoServer

texts:
I've got several layers in GeoServer, one of which I'm trying to use as a WMS in desktop programs (Manifold 8). I've set it up correctly, but whenever I bring it into Manifold, it projects the layer into EPSG:4326, which unfortunately can't work with my map - I need it in EPSG:3857.
The layer in question is wrapped inside a more overarching layer definition, which contains all the layers - and coordinate systems that allow them all to be used flexibly. However, I want to put this layer on its own in the capabilities document, and only allow it to use EPSG:3857.
How, in GeoServer, would I edit, or override, a capabilities document?


texts: Using SRTM Global DEM for Slope calculation?

texts:
I have downloaded SRTM GDEM (~90 km resolution).

I am using ArcGIS 10.

I have tried to use spatial analyst to compute for slope.

However, I can not compute for the slope.

The output values has only two ranges 0 and 0.1-90.

I am not really sure what is the problem?


texts: Transforming line into polygon with Geometry Generator

texts:
I would like to turn my line into polygon in the symbology.
I found some hint here: https://anitagraser.com/2017/04/08/a-guide-to-geometry-generator-symbol-layers/ but I am not sure of the expression, which should be included.
I tried something like this, but the fill is not visible at all. I've got a line only.

How can I make fill for these enclosed lines?


texts: Seeking Mississippi Township Range Section shapefiles?

texts:
Does anyone know where I can get a full complete T/R/S for Amita County in MS?  Geocomm shp files looks incomplete.  The counties in Mississippi are not helpful.

I ended up finding the TRS from this website http://www.maris.state.ms.us/HTM/DownloadData/Statewide.html


texts: My company I started working for is looking to move from ArcInfo and MapPoint to an online GIS for easier data access/availability.  Obviously the GIS capabilities needed are rather basic, but security is a must.  What sort of options are out there?  In order to ensure security of data, would we need a hosted platform?  Any help would be appreciated!


texts: Calculating minimum/maximum width for overlay of two polygon strips in ArcGIS or some Opensource GIS

texts:
Is there some solution for calculating minimum and maximum width for overlay of two polygon strips? When I clip these strips the result is another strip with variable width about 1 km and 20 km long. How can I get minimum and maximum width for the clipped strip in shorter axis? I try script from Calculating Polygon Width in ArcGIS Desktop? but it doesn´t suite with this issue.
Is there some solution for ArcGIS 10.1? example: clipped strip with red line for minimum width, blue for maximum
Or is it some solution solving the problem straightway without clipping inputs?


texts: How to save buffer created with ArcObjects to a database?

texts:
I have created a buffer as follows

IMap map = axMapControl1.ActiveView.FocusMap;
IFeatureSelection fSelection = map.get_Layer(0) as IFeatureSelection;
ISelectionSet selSet = fSelection.SelectionSet;
ICursor cursor;
selSet.Search(null, true, out cursor);

IFeature feat = ((IFeatureCursor)cursor).NextFeature();
ITopologicalOperator topoOpr = feat.Shape as ITopologicalOperator;

IPolygon buffPolygon = topoOpr.Buffer(100) as IPolygon;

IActiveView activeView=axMapControl1.ActiveView;
IScreenDisplay screendispaly = activeView.ScreenDisplay;
screendispaly.StartDrawing(screendispaly.hDC, (System.Int32)esriScreenCache.esriAllScreenCaches);

ISimpleFillSymbol simFillSymbol = new SimpleFillSymbolClass() ;
simFillSymbol.Style = esriSimpleFillStyle.esriSFSSolid;
IRgbColor color = new RgbColorClass();
color.Red = 255;
simFillSymbol.Color = color;

ISymbol symbol = simFillSymbol as ISymbol;
screendispaly.SetSymbol(symbol);
screendispaly.DrawPolygon(buffPolygon);

Now I can see a buffer around a feature and I want to save it as a feature class into a database.

How can I do it?


texts: Generating tiles from multiple rasters - gdal2tiles.py

texts:
I have 3 rasters, each is in the same SRS but with differing extents

I want to combine them into a single tile set using gdal2tiles.py

eg:-

gdal2tiles.py -z '2-4' -e -p raster -r average --s_srs EPSG:27700 raster1.tif /tiles
gdal2tiles.py -z '5-8' -e -p raster -r average --s_srs EPSG:27700 raster2.tif /tiles
gdal2tiles.py -z '9-12' -e -p raster -r average --s_srs EPSG:27700 raster3.tif /tiles

When I try this the tiles are not centred correctly when zooming from 4 to 5 and then 8 to 9. I suspect this is due to inconsistent extents of the rasters. Or am I missing something else? Should I be wrapping them in a virtual raster


texts: Reading Null values in field using ArcPy SearchCursor?

texts:
I'm using a Python script with ArcGIS 10 to count the number of unique string values in a feature class. My script works if there are values in the field but if there are any NUll values in the field it crashes. Any Ideas?

My code is below:

import arcpy
import operator

inTable = arcpy.GetParameterAsText(0)
stringField = arcpy.GetParameterAsText(1)

stringList = [row.getValue(stringField) for row in arcpy.SearchCursor
                                       (inTable, "", "",stringField)]


texts: QGIS: How to append one Geopackage layer into another Geopackage layer?

texts: At 3.8.3, I need to to append one Geopackage layer to another Geopackage layer, both in the same Geopackage.  Let’s say my Geopackage has two layers, A and B.  How do I append B into A?


texts: Draw several lines to represent the distance to a certain point for several points


texts:
I have a few MULTIPOLYGON geometries that are very complex to work with, and I need to find POINT geometries from a ~40M rows that are inside these MULTIPOLYGON geometries.
Right now I'm working with ST_Within but things are going on extremely slowly.

I don't need all of the POINT geometries, just a few of them, but they all have to be inside of the MULTIPOLYGON (I need a few points per geometry).

I thought that if there was something reverse to an envelope, that is to find the largest box that is within a geometry.

I couldn't find it in the docs or Google or here.

Is it possible to create such a box? If so, how?

I do not want to use simplify since it is:

Not accurate enough.
Still slow.


texts:
I am working on a project that uses SQL server for geometry datatype and ArcGIS Server and asp.net website and entity framework for data access and ArcGIS JavaScript API for map viewer.

I want to split polyline object in equal distance for example if I have a polyline with 60 km I want to split it in 10 polyline with 6 km distance and send polylines on ArcGIS JavaScript map api.

I can use T-Sql and C# for splitting polyline.
How can I do it?


texts: PostGIS package problem on Ubuntu 12.04.2 LTS

texts:
I've started getting an error on Ubuntu when trying to upgrade packages:
"Not all updates can be installed
Run a partial upgrade, to install as many upgrades as possible"

When I run a partial upgrade, I get this:
"Could not calculate the upgrade
An unresolvable problem occurred while calculating the upgrade:
The package 'postgresql-9.1-postgis' is marked for removal but it is in the removal blacklist.
This can be caused by:
 * Upgrading to a pre-release version of Ubuntu
 * Running the current pre-release version of Ubuntu
 * Unofficial software packages not provided by Ubuntu

This is most likely a transient problem, please try again later."

I have PostGres 9.1.9 and PostGIS 1.5.3.

Sounds similar to this problem: https://gis.stackexchange.com/questions/67143/postgresql-9-1-postgis-2-0-ubuntu-12-04-installation-problem

Thanks

Mini


texts:
Since installing 2.14 (and manually fixing errors with MSYS / Grass 7 carryover from previous installation of 2.12), when I try to search for a command in processing toolbox it appears to be re-indexing the available tools and causes a significant delay (say 40 secs +) before starting to display matching results.

I also have advisories confirmign that I can add more algorithms to the toolbox (referring to Grass 6) which I don't need, given I have Grass 7 activated.

How can I improve performance in accessing processing toolbox algorithms?


texts: "Save As" truncating string fields to 255 characters

texts:
I'm trying to convert a postgis table to an ESRI FileGeodatabase with a nightly build of QGIS, using Save As which utilises OGR.

One of the fields is a text field with > 255 characters. I have noticed that the field was truncated to 255 characters and then QGIS did not export the records which had more than 255 characters in the field. I then found this part of the QGIS code which seems to truncate String fields to 255 characters.

I understand that this has always been a limitation for Shapefiles, but it is not for file geodatabases, so it seems bizarre to truncate all text even when the provider supports it. Or is there some other limitation I should be aware of?Maybe should go on QGIS-dev but thought I would get folks' view here first.


texts:
When creating a GeoPackage vector layer in QGIS 3.10.x, it provides the following geometry types: Point, Line, Polygon, Multipoint, MultiLine, MultiPolygon, CompoundCurve, CurvePolygon, MultiCurve, MultiSurface.

And when creating a shapefile-layer only the geometry points: Point, Line, Polygon, Multipoint are available.

How are the geometry types MultiLine, MultiPolygon, CompoundCurve, CurvePolygon, MultiCurve and MultiSurface represented when exporting a GeoPackage layer to a Shapefile?


texts:
I've calculated relative heights in ArcGIS, now I'd like to calculate relative heights (range between maximum and minimum cell value in window 3x3) in SAGA. It's important that I don't want standardized values (0-1).

There is a module called "Relative Heights and Slope Positions", it creates few outputs, but they don't correspond directly to results calculated in ArcGIS (Focal Statistics)

These are sample results from ArcGIS:

Does anyone know how to convert the results from SAGA into values like I've calculated in ArcGIS with mentioned module or other SAGA module?


texts: Extracting Multi Values to Points and Dissolve afterwards with ModelBuilder?

texts:
I want to extract Multi Values from 2 Rasters to Points and dissolve these points afterwards with an ArcGIS ModelBuilder model. In the Dissolve function I want to make some statistics of the extracted Values (Min. and Max.)

My Problen: In the statistic fields of Dissolve tool I can’t select the new generated fields of the "Extract Multi Values to Points"-Tool . The new generated fields can also not be selected for dissolving.
Questions:
1) Why are the new generated fields not offered in the dissolve tool?
2) How can I make my model run?


texts:
Is there a way to use IN filters (like the IN operator from SQL) in OpenLayers 2?

For example, let's say, I want a filter like this:

city = 'Paris' OR city = 'London' OR city = 'Lisbon'

In OpenLayers, I would have to combine logical and comparison filters and it would be a big filter for only 3 cities. It would be something like this:

var filter = new OpenLayers.Filter.Logical({
    filters: [
        new OpenLayers.Filter.Comparison({
            property: 'city',
            type: OpenLayers.Filter.Comparison.EQUAL_TO,
            value: 'Paris'
        }),
        new OpenLayers.Filter.Comparison({
            property: 'city',
            type: OpenLayers.Filter.Comparison.EQUAL_TO,
            value: 'London'
        }),
        new OpenLayers.Filter.Comparison({
            property: 'city',
            type: OpenLayers.Filter.Comparison.EQUAL_TO,
            value: 'Lisbon'
        })
    ],
    type: OpenLayers.Filter.Logical.OR
});

I was wondering if there is a way to do it like this:

var filter = new OpenLayers.Filter.Comparison({
    property: 'city',
    type: OpenLayers.Filter.Comparison.IN,
    values: ['Paris', 'London', 'Lisbon']
});

I know that it does not work with OpenLayers, but how can I make it work?


texts:
I want to obtain a square from a circle using simplification function in sf.

This is my circle:

With this code, I obtain a rhombus

plot(st_simplify(st_geometry(sfdf),dTolerance = 1000))

How can I obtain a square? in QGIS I manage to obtain it using SnapToGrid option in simplification process.


texts: Extracting latitude/longitude from javascript and passing to html form

texts:
I have the following problem

I use a google map inside an html form where the user can enter an address (autocomplete) and the marker goes to the place. The marker is draggable.

What I want is to extract latitude and longitude and to pass the values to two different hidden fields In order to end up to two text fields (decimal)

The JavaScript

// This example adds a search box to a map, using the Google Place Autocomplete
// feature. People can enter geographical searches. The search box will return a
// pick list containing a mix of places and predicted search terms.
// This example requires the Places library. Include the libraries=places
// parameter when you first load the API. For example:
//

function initAutocomplete() {

var map = new google.maps.Map(document.getElementById('map'), {
         center: {lat: 37.9371157, lng: 22.9303688},
         zoom: 13,
         mapTypeId: google.maps.MapTypeId.ROADMAP
});

        // Create the search box and link it to the UI element.
        var input = document.getElementById('pac-input');
        var searchBox = new google.maps.places.SearchBox(input);
        map.controls[google.maps.ControlPosition.TOP_LEFT].push(input);

        // Bias the SearchBox results towards current map's viewport.
        map.addListener('bounds_changed', function() {
          searchBox.setBounds(map.getBounds());
        });

        var markers = [];
        // Listen for the event fired when the user selects a prediction and retrieve
        // more details for that place.
        searchBox.addListener('places_changed', function() {
          var places = searchBox.getPlaces();

          if (places.length == 0) {
            return;
          }

          // Clear out the old markers.
          markers.forEach(function(marker) {
            marker.setMap(null);
          });
          markers = [];

         // For each place, get the icon, name and location.
          var bounds = new google.maps.LatLngBounds();
          places.forEach(function(place) {
            var icon = {
              url: place.icon,
              size: new google.maps.Size(71, 71),
              origin: new google.maps.Point(0, 0),
              anchor: new google.maps.Point(17, 34),
              scaledSize: new google.maps.Size(25, 25)
            };

            // Create a marker for each place.
            markers.push(new google.maps.Marker({
              draggable: true,
              map: map,
              icon: icon,
              title: place.name,
              position: place.geometry.location
            }));

            if (place.geometry.viewport) {
              // Only geocodes have viewport.
              bounds.union(place.geometry.viewport);
            } else {
              bounds.extend(place.geometry.location);
            }
          });
          map.fitBounds(bounds);
        });
      }

This code works well. Now i need to extract latitude and longitude values and to pass them in a hidden field inside an html form

input type="hidden" name="maps_latitude" id="maps_latitude" value=""
input type="hidden" name="maps_longitude" id="maps_longitude" value=""

maps_latitude and maps_longitude are two text fields in a cms (decimal) in order to get the values

Can anyone help how I need to modify the JavaScript and the hidden fields?


texts: gdal_polygonize returns with a strange error

texts:
i'm trying to execute gdal_polygoinize.py command.
I have a main python method that takes in input a raster file and a output shapefile, and then calls a method named polygonize.py that builds and try to execute the gdal_polygonize command.

The code is:

format = "ESRI Shapefile"
command_poligonize = "gdal_polygonize.bat " + \
                        + input_raster + ' '\
                         '-f ' + '"' + format + '"'+ \
                         shape_output
status=subprocess.call(command_poligonize, shell=True)

i tryed also with the os.system() method:

os.system(command_poligonize)

Executing the command:

print('DEBUG: comand: ' + command_poligonize)

i obtain the correct gdal_polygonize invocation, like:

gdal_polygonize.py input.tif -f "ESRI Shapefile" output.shp

The error description is the following:

  File "C:\OSGeo4W64\bin\gdal_polygonize.py", line 209, in <module>
     callback = prog_func )
 File "C:\Python35\lib\site-packages\osgeo\gdal.py", line 1636, in Polygonize
     return _gdal.Polygonize(*args, **kwargs)
 RuntimeError: Object given is not a Python function

Any suggestion?

I tryed to run the command by itself, like:

gdal_polygonize.py input.tif -f "ESRI Shapefile" output.shp

and it return with the following error:

  File "C:\OSGeo4W64\bin\gdal_polygonize.py", line 209, in <module>
     callback = prog_func )
 File "C:\Python35\lib\site-packages\osgeo\gdal.py", line 1636, in Polygonize
     return _gdal.Polygonize(*args, **kwargs)
 RuntimeError: Object given is not a Python function

I'm using Python35 and Gdal 2.0.2 released 2016/01/26.


texts: Detecting If A Coordinate Is Near Polygon Edges In Google Maps API

texts:
Is there a way to add a test to see if a Lat/Long coordinate is near the edge of a polygon?

I am building very complicated polygons that have rivers and roads as borders. Because of this, it is very difficult to make the polygons 100% accurate.

A test like this could add appropriate warnings in these scenarios.


texts: SP1 for ArcGIS 10.1 has extension for 64-bit Python for background geoprocessing -- implications for other Python use

texts:
If the 64-bit Background Geoprocessing for ArcGIS Desktop 10.1sp1 is installed how does that affect Python configuration for general use?

Esri has released SP1 for the ArcGIS 10.1 product, and for the Desktop component has made available a 64-bit Python 2.7.2 instance that when installed replaces Esri 32-bit Python 2.7 calls with 64-bit calls for background geoprocessing (processing sent out side of an active ArcGIS Desktop GUI session).

Installation is optional--i.e. without it ArcGIS Desktop continues with only 32-bit user memory space available.  With it, a \bin64 entry is made to the Esri program folder, and registry keys are established for 64 bit processing at "HKEY_LOCAL_MACHINE\SOFTWARE\ESRI\ArcGIS 10.1 SP1 for Desktop Background Geoprocessing (64-bit)" and background geoprocessing is then done in 64-bit memory space.

Both Esri's 32-bit and now its 64-bit Python are full 2.7.2 installations.

During background geoprocessing Esri's registry, path and configuration files take care of calling the 64-bit Python executables. While use of 32-bit Python is retained in the ArcGIS Desktop session where needed.

So, what is left to each user is decide whether to integrate the installed 64-bit Python 2.7.2 into their computing needs outside of ArcGIS. And how best to do that when both 32 and 64 bit versions of Python are resident on the system?

Asked another way--Esri takes care of integration of 64-bit Python within its ArcGIS Desktop, so what ways will work to configure it for the rest of the use cases beyond ArcGIS?


texts: Add row and column name to square grid feature in ArcGIS

texts:
I have a square grid feature class of 250 x 250 tiles.  This will serve as a tile index for a a project.  They currently have no identifying attributes other than a unique ID.  I need to add a field to rename them by their row and column, e.g. "R001_C001" or something.

Is there a way to do this easily?  I'm okay with python.


texts: IndexError: list index out of range - item report ArcGIS API for Python

texts:
I'm quite new to Python and ArcGIS API for Python.
I'm trying to create a report of all the items owned by one of my ArcGIS online organization user.
For each item owned by the user I want to record some item info in a Pandas dataframe, then export a CSV as a report.
# Import the modules
import os
import pandas as pd
from arcgis.gis import GIS

# Login to arcgis online
gis = GIS(url='https://myorganization.com',
         username= 'myUsername')

# Define the user
user = 'UserToCheck'

# Search for contents owned by a specific user
all_contents = gis.content.search(query="owner:" + user,
                                    max_items=3000
                                   )

# Check the number of contents
len(all_contents)

# Define the columns title to insert in the dataframe we want to export
title = 'TestTitle',
description = 'testDescription'
item_id = 'TestID'
item_type = 'TestItemType'
date_creation_unix = int("1415400560000")
date_modified_unix = int("1415969577000")
url = 'TestUrl'
usage = int("142348")

# Define the columns
columns = {'content title': title,
          'description': description,
          'id' : item_id,
          'item_type' : item_type,
          'date_creation_UNIX' : date_creation_unix,
          'date_modified_UNIX' : date_modified_unix,
          'url' : url,
          'usage_last_year': usage
          }

columns

# Create the dataframe using the columns defined in the previous step
df = pd.DataFrame(data=columns)

# print the dataframe
print(df)

x = 0
for content in all_contents:
    df_usage = content.usage('1Y') # create a report dataframe about last year (it include 2 columns, "Day" and "Usage". I'm only interested about the "Usage" and i need to sum this value along all the year)
    sum_usage = sum(df_usage["Usage"]) # store only the sum of "Usage" column

    to_append = (content.title,
                content.description,
                content.id,
                content.type,
                content.created,
                content.modified,
                content.url,
                sum_usage
                )

    # append "to_append" as a row to the dataframe
    df_lenght = len(df)
    df.loc[df_lenght] = to_append
    x = x + 1
    print (x, " raw added to the dataframe")

it seems working well, but sometimes I have the following error:
IndexError                                Traceback (most recent call last)
In  [24]:
Line 5:     df_usage = content.usage('1Y') # create a report dataframe about last year (it include a row per day)

File C:\Users\Me\AppData\Local\Programs\ArcGIS\Pro\bin\Python\envs\arcgispro-py3\lib\site-packages\arcgis\gis\__init__.py, in usage:
Line 9060:  if not self.layers[0].container:

IndexError: list index out of range
---------------------------------------------------------------------------

I noticed that the first time the error refers to a certain row (120), and the first 119 rows were well recorded in the dataframe. The second time, using the same script, happened with another row (somewhere around 200)
Why do I have this error and how could I fix it?


texts: Avoiding small sliver areas and Reducing pages in ArcGIS Pro Map Series?

texts:
Automating some map production.

Current high level approach:

Have user providing an Area of Interest (AOI) as a polygon extent
If AOI > Threshold Area size Then generate a grid map
Use the origin of the AOI, create number of grids based upon PageUnit size for an A4 landscape template
Then use this Grid to generate a Map Series.

Problem

There can be some pages that will only show a small 'sliver' of the intersecting AOI.  As crudely illustrated by the red rectangles in the below image:

What is a better approach to automate this kind of map production.  Make the assumption that the incoming AOI can be any size, any shape.  Scale and layout orientation are flexible.

The final multi-page map book needs to be easily understood by end-user, hence why we assumed a simple grid map would be sufficient.

I anticipate others will have encountered this issue before and found good workarounds


texts: In PostGIS, does casting an `ST_Envelope` to a geography take the earth's curvature into account?

texts:
Using PostGIS, I have a table of data with geographic points defined as last_known_coordinates GEOGRAPHY(POINT,4326). I want to find all of these within a geographic bounding box with given top left and bottom right geographic coordinates, and be sure that the calculations are geographic, not geometric, e.g. that the results account for the curvature of the earth.
The query I'm using is roughly:
SELECT * from table WHERE last_known_coordinates && ST_MakeEnvelope(?, ?, ?, ?, 4326)::geography

As you can tell I'm a bit new to this, and I know what I want, but not necessarily how to do it correctly. Other answers hint that this is correct, but I'd like to know for sure: Does the ST_Envelope produced by ST_MakeEnvelope, when cast to a geography, take the curvature of the earth into account?


texts:
I love the idea of ArcGIS Pro, especially for non-GIS analysts to get one-off maps done easily and quickly.

The big question is: Does ArcGIS Pro join to excel worksheets like ArcMap does?


texts: How do you merge overlapping .asc files in QGIS with topography and bathymetry?

texts:
This is my first time trying to manipulate GIS data.

I have two files. The first is higher resolution (SRTM) topography data for a region of islands (the Caribbean), I don't know what the surrounding water values are set to.

The second file has lower resolution topography and bathymetry data over a similar area (ETOPO1).

What I want to do is: combine the two files so that I have the higher resolution topography data (from the first file) surrounded by the lower resolution bathymetry data (from the second file).

Is this possible? If so, how? I have spent ages trying various combinations of merge's, clip's and masks to no avail.

Maybe some further complication are that I only need the higher resolution data for one island of interest, and that the first file has a much larger spatial extent than the second (most of which I don't need due to the single island requirement).

Hopefully that all makes sense, and someone knows how to do what I'm asking.


texts: How can I generate Kernel density raster/stars with the same dimensions

texts:
I would like to generate different rasters from kernel density functions with the same resolution, so that I can compare then.
Take this example from here:
library(MASS)
library(raster)
set.seed(123); xy = cbind(rnorm(1000), rnorm(1000))
k = kde2d(xy[,1],xy[,2],h=.4,n=100)
r1 = raster(k)

xy = cbind(rnorm(1000), rnorm(1000))
k = kde2d(xy[,1],xy[,2],h=.4,n=100)
r2 = raster(k)

par(mfrow = c(1,2))
plot(r1); plot(r2)

r1 - r2
#Error in compareRaster(e1, e2, extent = FALSE, rowcol = FALSE, crs = TRUE,  :
#  different origin

Obs: I know I can resample the rasters, but I would like to define a standard raster/stars object
Obs2: I can compare as stars objects. Is it correct?
s1 <- st_as_stars(r1)
s2 <- st_as_stars(r2)
s3 <- s1 - s2

plot(s3, axes = T, breaks = "equal")

Thanks


texts: How to extract clear cut areas from dried trees in forests

texts: I have 40 cm resolution aerial image with visible and near-infrared bands. Purpose is to extract dried trees in forest. I used multiresolution segmentation algorithm to extract individual trees and used NDVI to classify dried and healthy trees, but there are clear cut areas in the study area, and they give the same NDVI value as dried trees. I even checked RGB and NIR values of segments and they also give the same value as dried trees. Any suggestions?


texts: GeoDjango: Order by distance for model that has a ForeignKey relation to another model that has a PointField

texts: Model A has a ForeignKey relation to Model B. Model B has a PointField. Is there a way to order Model A objects by distance from a POINT?


texts: Reprojecting raster using gdalwarp "EPSG_Code: 6204"

texts:
I want to reproject a raster using the projection of other raster using gdalwarp, however, I am confronting some problems because the raster does not line up when loading them into a GIS Software.

This is my code:

from osgeo import gdal,osr
out_ras_500m_buff = r'MK_020m_nat_buffer500m_v2.tif'
boundaries_in = r'MK_020m_nat_buffer100m_v2.tif'
ds=gdal.Open(boundaries_in)
prj=ds.GetProjection()
print(prj)
gt = ds.GetGeoTransform()
cols = ds.RasterXSize
rows = ds.RasterYSize

srs=osr.SpatialReference(wkt=prj)
if srs.IsProjected:
    print(srs.GetAttrValue('projcs'))
print(srs.GetAttrValue('geogcs'))

ext = GetExtent(gt, cols, rows)
print(ext)
xmin = ext[0][0]
ymin = ext[2][1]
xmax = ext[2][0]
ymax = ext[0][1]
print(f"XMin:{xmin}")
print(f"YMin:{ymin}")
print(f"XMax:{xmax}")
print(f"YMax:{ymax}")

EPSG_CODE = "6204"
warped_raster = out_ras_500m_buff.replace(".tif","_" + EPSG_CODE+ ".tif")
warped_raster = warped_raster.replace("02", "03")

prj_file = OutDir + "/" + os.path.basename(out_shp).replace(".shp", ".wkt")

with open(prj_file, "w") as f:
    f.write(prj)
#subprocess.call("gdalsrsinfo -o wkt {raster_boundaries_nat} > target.wkt")

#gdal.Warp(warped_raster, out_ras_500m_buff, dstSRS= f"EPSG:{EPSG_CODE}")

#subprocess.call(f'gdalwarp -t_srs {prj} {out_ras_500m_buff} {warped_raster} -te {xmin} {ymin} {xmax} {ymax}')
command =  f'gdalwarp -t_srs {prj_file} {out_ras_500m_buff} {warped_raster} -te {xmin} {ymin} {xmax} {ymax} -tr 10 10 -r near -co "COMPRESS=LZW"'
subprocess.call(command)

I am using a GetExtent function which get the extent of the raster I want to use to reproject the other one:

def GetExtent(gt,cols,rows):
    ''' Return list of corner coordinates from a geotransform

        @type gt:   C{tuple/list}
        @param gt: geotransform
        @type cols:   C{int}
        @param cols: number of columns in the dataset
        @type rows:   C{int}
        @param rows: number of rows in the dataset
        @rtype:    C{[float,...,float]}
        @return:   coordinates of each corner
    '''
    ext=[]
    xarr=[0,cols]
    yarr=[0,rows]

    for px in xarr:
        for py in yarr:
            x=gt[0]+(px*gt[1])+(py*gt[2])
            y=gt[3]+(px*gt[4])+(py*gt[5])
            ext.append([x,y])
        yarr.reverse()
    return ext

Although the result contains the CRS that I want, the projected raster does not have the following information that is found in a .prj file for the target spatial reference:

PROJCS["Macedonian_State_Coordinate_System",
    GEOGCS["GCS_MGI_1901",
        DATUM["D_MGI_1901",
            SPHEROID["Bessel_1841",6377397.155,299.1528128],
            TOWGS84[-521.7476,-229.4892,-509.9207,-4.02878,-4.48836,15.52067,-9.7803]],
        PRIMEM["Greenwich",0.0],
        UNIT["Degree",0.0174532925199433]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["False_Easting",500000.0],
    PARAMETER["False_Northing",0.0],
    PARAMETER["Central_Meridian",21.0],
    PARAMETER["Scale_Factor",0.9999],
    PARAMETER["Latitude_Of_Origin",0.0],
    UNIT["Meter",1.0],
    AUTHORITY["EPSG",6204]]

My question is:

TOWGS84[-521.7476,-229.4892,-509.9207,-4.02878,-4.48836,15.52067,-9.7803]]

Given that information, how can I use the .prj file to reproject a raster with the correct parameters (shift to WGS84)?

It seems that only giving the EPSG CODE in the subprocess.call does not make the work.

I found that using pyproj, and setting a string like the following, one can accomplish that. However, I cannot understand how to plug the shift to wg84 parameters when reprojection a raster.

isn2004=pyproj.Proj("+proj=lcc +lat_1=64.25 +lat_2=65.75 +lat_0=65 +lon_0=-19 +x_0=1700000 +y_0=300000 +no_defs +a=6378137 +rf=298.257222101 +to_meter=1")


texts: In a Network is Valve a Line or Point in GIS Standard?

texts: I am designing a GIS Database related to Gas Network. In my opinion Valve equipment is a link equipment which can be closed or opened, So it should be a Line. But lots of clients gather this data as Point layer, which cause lots of problem when analysis the network. Is there any standard that clarify this issue?


texts: Geoserver does not return JSON when creating a feature

texts:
I'm trying to convince Geoserver that it returns JSON after creating a feature in WFS:

POST http://localhost/ows?outputFormat=application/json&exceptions=application/json

BODY:

<wfs:Transaction xmlns:wfs="http://www.opengis.net/wfs" service="WFS" version="1.1.0" 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">
    <wfs:Insert>
        <feature:GEOIP_NETWORKS xmlns:feature="http://edit.tmb.cat">
            <feature:GEOMETRY>
                <gml:Point xmlns:gml="http://www.opengis.net/gml" srsName="EPSG:4326">
                <gml:pos>2.167296 41.381919</gml:pos></gml:Point>
            </feature:GEOMETRY>
            <feature:NETWORK>AAA</feature:NETWORK>
            <feature:DESC_NETWORK>asdasd</feature:DESC_NETWORK>
       </feature:GEOIP_NETWORKS>
    </wfs:Insert>
</wfs:Transaction>

However it seems to ignore outputFormat parameter and keeps returning XML:

<?xml version="1.0" encoding="UTF-8"?>
<wfs:TransactionResponse xmlns:xs="http://www.w3.org/2001/XMLSchema" xmlns:wfs="http://www.opengis.net/wfs" xmlns:gml="http://www.opengis.net/gml" xmlns:ogc="http://www.opengis.net/ogc" xmlns:ows="http://www.opengis.net/ows" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" version="1.1.0" xsi:schemaLocation="http://www.opengis.net/wfs http://tmbapi.tmb.cat:80/geoserver/schemas/wfs/1.1.0/wfs.xsd">
    <wfs:TransactionSummary>
        <wfs:totalInserted>1</wfs:totalInserted>
        <wfs:totalUpdated>0</wfs:totalUpdated>
        <wfs:totalDeleted>0</wfs:totalDeleted>
    </wfs:TransactionSummary>
    <wfs:TransactionResults/>
    <wfs:InsertResults>
        <wfs:Feature>
            <ogc:FeatureId fid="GEOIP_NETWORKS.19"/>
        </wfs:Feature>
    </wfs:InsertResults>
</wfs:TransactionResponse>

The exceptions work fine though, does return JSON.

Is it just not implemented yet?

I'm using Geoserver 2.11.1.


texts: Installing r.denoise on GRASS GIS (Windows) doesn't work

texts:
I try to install the r.denoise add on for GRASS GIS on my Windows machine. If I start the installation without copying the html manual by hand, I get the errnr 2 that this manual can't be found. After copying it to the right place, it doesn't find its binary. Therefore, I created a .bat file to start the .py file.

Right now, it doesn't find the binary, but if I copy/paste the content of the .bat into GRASS console, Python doesn't find my pyproj library (although installed for my home python installation).

Is there any workaround to get this one working?


texts: Generated ArcGIS tile cache ignores DEM hillshade effect?

texts:
I'm attempting to generate a tile cache in ArcGIS 10.3 where one of the input layers is a DEM raster. The resulting tile cache displays all layers as expected, except for the raster layer. The output tile images show the raster as if it had all the correct settings but simply seems to ignore the hillshade effect that has been applied. The images below illustrate how the raster is originally symbolised (left) vs how it is showing in the tile cache (right).

Min-max stretch - Expected vs Actual

No stretch - Expected vs Actual

The resulting tiles do seem to take into consideration at least some of the stretched raster settings (e.g. the gradient, stretch type, etc), just not the hillshade.

Should this be happening?


texts: How to export layers with rule based symbology with qgis2web

texts:
I am working as volunteer on a simple map for an NGO, I have some point layers with rule-based symbology. When I load these layers in qgis2web the symbology disappears, both in Leaflet and OpenLayers. I think I have the same problem posted by Sophia R. and solved (Can qgis2web display rule-based layer style?).

The matter is that I'm new in QGIS and I don't know how to "preselect your layer with "X" = 'a' and export this selection as a new layer" as suggested by Riccardo.

How can I do this?


texts: I got two different number of items in a merged or dissolved layer

texts: I'm using QGIS and I have the locations toxic waste sites in the U.S. on my map. I also have the weather warnings as a polygon on my map. I want the number and list of toxic waste sites in those polygons. The first time, I merged the weather warnings into one file, and clipped the toxic waste sites in the area of weather warnings, with the toxic waste sites as the input layer and the merged weather warning polygon file as the overlay. The second time, I dissolved the weather warning polygon, and did "count points in polygon." I'm getting different numbers. Is there a reason why this is happening?


texts: Changing description of layer using ArcPy?

texts:
I would like to change the file description with a python code on ArcMap (10.6). I could not find a way to do it or resource explaining it.

Do you have any idea?


texts: Getting polygon shapefile node coordinates and point order

texts:
I am trying to get the coordinates of each node of a polygon in QGIS. What I basically did is use the following tools:

Vector -> Geometry tools -> Extract nodes

Vector -> Geometry tools -> Export Geometry column

However only around the first 300 points had successfully obtained coordinate values in the new xcoord and ycoord fields. The other three thousand plus points had null values. What can I do to get the other points' coordinates? Also is there an easy way to get the point order of each node?


texts: Color style in the RGB field without “#” or how to copy part of the text from one field to another in the layer attributes?

texts:
How to insert the color of the style without the key "#" in the field 'CLRRGB' (Reefmaster does not understand "#" but QGIS does not understand when without "#"). Forced to leave two attributes. With the help of the script, I fill in the layer style in 'QGIS_COLOR'. I need another script so that in the field 'CLRRGB' was the same color only without the "#".

Perhaps craftsmen can improve this script:

prefix = "'"
layer = iface.activeLayer()
attr = layer.renderer().classAttribute()
attrColor = 'QGIS_COLOR' # Название поля для хранения цветов
fieldIndex = layer.dataProvider().fieldNameIndex(attrColor)
attrFeatMap = {}

for cat in layer.renderer().categories():
  expr = "\""+attr+"\"="+prefix+unicode(cat.value())+prefix
  for f in layer.getFeatures(QgsFeatureRequest(QgsExpression(expr))):
    attrMap = { fieldIndex : cat.symbol().color().name()}
    attrFeatMap[ f.id() ] = attrMap

layer.dataProvider().changeAttributeValues( attrFeatMap )


texts: OpenLayers Multiple Maps Problem (maps leave blank)

texts:
I try a round with OpenLayers and was very impressed how mighty it is.

But I have a problem (maybe a easy one, because I am not so trained in js).

After adding additional maps all maps leave blank, and I don't know why. Maybe I need a simple hind

    <script src="http://openlayers.org/api/OpenLayers.js" ></script>
    <script type="text/javascript">

    var supplyMap, demandMap, voluminaMap;
    var supplyMapLayer, demandMapLayer, voluminaMapLayer;
    var supplyMapMarkers, demandMapMarkers, voluminaMapMarkers;

    function init(){
        var wms_all = new OpenLayers.Layer.WMS("OpenLayers WMS", "http://vmap0.tiles.osgeo.org/wms/vmap0", {layers: 'basic'});

        //set supplyMap
        supplyMap = new OpenLayers.Map('supplyMap', {});
        supplyMapLayer = wms_all.clone();
        supplyMapMarkers = new OpenLayers.Layer.Markers( "SupplyMarkers" );

        supplyMap.addLayers([supplyMapLayer, supplyMapMarkers]);
        supplyMap.setCenter(new OpenLayers.LonLat(10.5,51.25),6.5);

        //set demandMap
        demandMap = new OpenLayers.Map('demandMap', {});
        demandMapLayer = wms_all.clone();
        demandMapMarkers = new OpenLayers.Layer.Markers( "DemandMarkers" );

        demandMap.addLayers([supplyMapLayer, supplyMapMarkers]);
        demandMap.setCenter(new OpenLayers.LonLat(10.5,51.25),6.5);

        //set voluminaMap
        voluminaMap = new OpenLayers.Map('voluminaMap', {});
        voluminaMapLayer = wms_all.clone();
        voluminaMapMarkers = new OpenLayers.Layer.Markers( "VoluminaMarkers" );

        voluminaMap.addLayers([voluminaMapLayer, voluminaMapMarkers]);
        voluminaMap.setCenter(new OpenLayers.LonLat(10.5,51.25),6.5);

       // setMarker(52.519945, 13.404713);
       // setMarker(48.141176, 11.568800);

    }

  function setMarker(lat, lon){

        var markerVal = new OpenLayers.LonLat(lon, lat);
        var size      = new OpenLayers.Size(21,25);
        var offset    = new OpenLayers.Pixel(-(size.w/2), -size.h);
        var icon      = new OpenLayers.Icon('http://www.openlayers.org/dev/img/marker.png',size,offset);

        markers.addMarker(new OpenLayers.Marker(markerVal,icon.clone()));

  }

  function clearMarker(){
       markers.clearMarkers();
  }

  </script>
    <style>
      @media screen
      {
        .supplyMap{width: 850px; height: 725px; border: 2px solid black;}
        .demandMap{width: 850px; height: 725px; border: 2px solid black;}
        .voluminaMap{width: 850px; height: 725px; border: 2px solid black;}
      }
    </style>


texts:
I'm missing Processing toolbox in my QGis 3.4 (Mac-version). I know the feature can be activated when going to the Plugins and clicking it on. But it's not there. I've tried to solve this thing for one week and I'm desperate. Processing toolbox is the key feature in the whole program and I cannot get further in my practice.

I don't know any script and I'm pretty new with QGIS.


texts: Alter Esri geoprocessing environment settings with C#

texts:
I maintain an extension to ArcMap (currently at 10.2.2) for coworkers' internal use.  Several people lately have been tripped up by un-set spatial grid settings when trying to export data into a file geodatabase from our non-SDE SQLServer database. My extension isn't involved in the export itself--rather that happens using standard Esri tools. Consequently I need to set the "Output Spatial Grid" size settings once universally. I'd like to include a check at startup that will look at these settings as they are stored on disk and change 0,0,0 to some useful values.  No user involvement is needed or wanted.  I've perused the ArcObjects doc, scoured my registry, and Googled for an entire afternoon and ended up having to ice my Googler . . . no luck.  Can anyone provide information on either:

Where on disk the "Output Spatial Grid" settings are stored, and/or
ArcObjects references for how to change them

I've tried the following with no luck.

    var gp = new ESRI.ArcGIS.Geoprocessing.GeoProcessor();
    gp.SetEnvironmentValue("spatialGrid1", 50);
    gp.SetEnvironmentValue("spatialGrid2", 200);
    gp.SetEnvironmentValue("spatialGrid3", 800);

But even if that works, I believe it would only alter the settings for the current session.  What I'd prefer is to write the values to disk as if the user had made the settings changes manually in ArcCatalog.


texts: Converting classification image to aggregated polygon in shapefile using gdal/python?

texts:
My goal is convert a 2D numpy array to shapefile. The matrix results from a classification, so the numbers have a meaning.

I solved my problem, but the solution is somehow convoluted and complicated.

Does anyone has a better implementation?

The actual code is too long lines (you cand find it here implemented as a function array -> shapefile) for a relatively simple action.

One of the big problem is gdal_polygonize giving separate polygons for disjoint regions.

gdal_calc.py* seems a pretty good solution, but it isn't in the osgeo module.

Here an example of dummy data from (Raster Layers — Python GDAL/OGR Cookbook 1.0 documentation)*

plot_matrix = np.array([[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
                    [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
                    [ 0, 1, 1, 1, 1, 0, 2, 2, 2, 2, 0, 1, 1, 1, 0, 2, 0, 0, 0],
                    [ 0, 1, 0, 0, 0, 0, 0, 2, 0, 2, 0, 1, 0, 1, 0, 2, 0, 0, 0],
                    [ 0, 1, 0, 1, 1, 0, 0, 2, 0, 2, 0, 1, 1, 1, 0, 2, 0, 0, 0],
                    [ 0, 1, 0, 0, 1, 0, 0, 2, 0, 2, 0, 1, 0, 1, 0, 2, 0, 0, 0],
                    [ 0, 1, 1, 1, 1, 0, 2, 2, 2, 2, 0, 1, 0, 1, 0, 2, 2, 2, 0],
                    [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
                    [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
                    [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]])

The data represent 2 disjoint region, with 2 different values of the classification (1 and 2) plus a background of "No Data" values, represented by 0s.

My actual solutions are :

Workflow #1 (Write some files around)

write the 2D numpy array to GeoTiff

use (gdal_polygonize.py)* to convert to "vector polygons for all connected regions of pixels in the raster sharing a common pixel value" (see documentation). This produce 5 different polygons, for the 5 regions in my data.
open the produced shapefile and aggregate the polygon sharing the same value in multi polygon collection, in my case there will be 2 MultiPolygon (0 represent no data in the original data)

Workflow #2 (everything in memory)

The same as above, but all driver in 1 is 'MEM' , so no temporary .tiff is created.
Also done  in memory (ogr driver 'Memory' )
Same as above.

The result is a shapefile with 2 MultiPolygon, each with an unique value of the Field "class"

I need at least 10 reputation to post more than 2 links, so the references are all here:

Raster Layers — Python GDAL/OGR Cookbook 1.0 documentation : pcjericks.github.io/py
gdalogr-cookbook/raster_layers.html#create-raster-from-array
gdal_polygonize.py: www.gdal.org/gdal_polygonize.html
gdal_calc.py : www.gdal.org/gdal_calc.html


texts: Summing thousands of rasters in GRASS GIS (after reach r.series limit)

texts:
I am using GRASS GIS 7.6 and want to sum 30.000 small rasters from a total viewshed calculation. I tried the following approach: 1-> r.list to get a list of the names of the rasters. 2 -> r.series (sum) to sum all the rasters, but there is a limit on the number of rasters that I can sum in each take...

Can anyone provide help with this?

Maybe the solution is some python coding which loops the sum.


texts: Google Earth Engine Image (Error) User memory limit exceeded. How can I fix it?

texts:
This is the code that has an error. I don´t know how to fix it.
I am trying to get the albedo values of a range of time: 2001-2018 and just of 4 months of each year, September to December.
var col = ee.ImageCollection('MODIS/006/MCD43A3')
                  .filter(ee.Filter.date('2001-01-01', '2018-12-31'))
                  .filterBounds(geometry).select('Albedo_WSA_shortwave');

// Functions to stack colletion series into image bands

var stackCollection = function(collection) {
  var first = ee.Image(collection.first()).select([]);
  var appendBands = function(image, previous)
  {return ee.Image(previous).addBands(image);};
  return ee.Image(collection.iterate(appendBands, first));};

// Function to generate mask for each Scene.

var getSceneMask = function(scene) {
  var dateString = ee.Date(scene.get('system:time_start')).format('yyyy-MM-dd');
  var mask = scene.select('Albedo_WSA_shortwave').gt(-9999).updateMask(ee.Image(1));
  return mask.rename(dateString);};

//Stacked mask for collection

var Maskcol = stackCollection(col.map(getSceneMask));

// Generate a datelist for collection

var datelist = ee.List(Maskcol.bandNames()).map(function(date){return ee.Date(date).format('yyyy-MM-dd')});

// Stack bands and rename by date

var WSA_stacked = stackCollection(col.select('Albedo_WSA_shortwave')).multiply(0.001).rename(datelist).clip(geometry);

Map.addLayer(WSA_stacked.select(20), {min:0.0, max:400, gamma:2});

print(WSA_stacked);
print(datelist);

// Exportar datos //

Export.image.toDrive({
  image: WSA_stacked.clip(geometry),
  description: 'MODIS_WSA_stacked',
  region: geometry,
  folder: 'ejemplo',
  scale: 500,
  crs: 'EPSG:4326'});


texts: Changing the path to the untwine executable in QGIS? - Point cloud data will not display in QGIS 3.18.1

texts:
Describe the bug
QGIS will load my LAZ file, but it will only show the file extents. It won't show points in 2D or 3D map. I have installed and uninstalled QGIS several times and I have used the rebooted OSGEO4W installer as well as the Standalone installers from OSGeo4W testing packages (MSI).
How to Reproduce
Drag and drop an LAZ file from the browser onto the map. The bounding box will be displayed immediately. Here is the LAZ file I used - nb_2017_2418000_7497000
The ept folder was created (ept_nb_2017_2418000_7497000) but it was empty.
The log messages for Point clouds included:
2021-04-10T13:20:22 INFO Using executable C:/Program Files/QGIS31~1.1/apps/qgis/./untwine.exe

However, the actual path on my laptop is "C:\Program Files\QGIS 3.18.1\apps\qgis\untwine.exe". I noticed this problem with the untwine path in the C:\OSGEO4W install and the C:\Program Files\QGIS 3.18.1 install.
Can I change the path to the untwine executable in QGIS?
QGIS and OS versions
QGIS version
3.18.1-Zürich
QGIS code revision
202f1bf
Compiled against Qt
5.15.2
Running against Qt
5.15.2
Compiled against GDAL/OGR
3.2.2
Running against GDAL/OGR
3.2.2
Compiled against GEOS
3.9.1-CAPI-1.14.2
Running against GEOS
3.9.1-CAPI-1.14.2
Compiled against SQLite
3.35.2
Running against SQLite
3.35.2
Compiled against PDAL
2.2.0
Running against PDAL
2.2.0 (git-version: 24dd45)
PostgreSQL Client Version
13.0
SpatiaLite Version
5.0.1
QWT Version
6.1.3
QScintilla2 Version
2.11.5
Compiled against PROJ
8.0.0
Running against PROJ
Rel. 8.0.0, March 1st, 2021
OS Version
Windows 10 Version 1709
Active python plugins
db_manager;
MetaSearch;
processing
Windows 10 Enterprise


texts: Shapefile lines disappear, but data remains

texts: I'm working in QGIS 3.4 on a line type shapefile layer with +/-5000 records.  I made an edit to the start/end nodes of a couple of lines to move them to a point node from another layer; then when I saved my line layer, half of the lines from that layer disappeared from my canvas.  The records still exist in the attribute table, I can select them based on their values, but it appears as if the geometry has gone null - even with a record selected, I can not zoom or flash features.


texts: Is it possible to give offset to symbol of point feature  using sld in openlayers or geoserver?


texts: AutoIdentifying spatial reference from shapefile fails on second attempt

texts:
I read a shapefile with GDAL 2.4.3 in C# and try to determine the used spatial reference. This works fine for the first time reading the file. When I close the file and read it again, the determination fails. Here is my code for reproduction:

GdalConfiguration.ConfigureGdal();
GdalConfiguration.ConfigureOgr();
if (!GdalConfiguration.Usable)
  throw new Exception("GdalConfiguration failed");

var filename = @"someshapefile.shp";

for (int i = 1; i <= 2; i++)
{
  using (var driver = OSGeo.OGR.Ogr.GetDriverByName("ESRI Shapefile"))
  using (var dataSource = driver.Open(filename, 0))
  using (var layer = dataSource.GetLayerByIndex(0))
  {
    layer.ResetReading();
    using (var feat = layer.GetNextFeature())
    using (var geom = feat.GetGeometryRef())
    {
      var spatRef = geom.GetSpatialReference();
      spatRef.AutoIdentifyEPSG();
      spatRef.ExportToWkt(out var wkt);
      Trace.WriteLine($"{i}. WKT = {wkt}");
      if (!int.TryParse(spatRef.GetAuthorityCode(null), out int srid) || srid <= 0)
        throw new Exception("unknown coordinate system."); // <= will be thrown when i = 2; srid is 0 then
    }
  }
}

The output differs for both attempts:

  WKT = PROJCS["ETRS89 / UTM zone 32N (zE-N)",GEOGCS["ETRS89",DATUM["European_Terrestrial_Reference_System_1989",SPHEROID["GRS
  1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6258"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4258"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",9],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",32500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","4647"]]

which is correctly identified as 4647.

  WKT = PROJCS["ETRS_1989_UTM_Zone_N32",GEOGCS["GCS_ETRS_1989",DATUM["European_Terrestrial_Reference_System_1989",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",32500000.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",9.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]

which results in srid == 0.

Any ideas?


texts: browser crashes when loading legends at startup

texts:
After much perusing for solutions, I've decided to post my problem as a question. I'm trying to load the legends that are built into the ArcGIS Service Layers I'm using. I'm loading about 7 Dynamic Service Layers and 7 Feature Layers. Each one of the Dynamic Service Layers has it's own legend whereas, each of the Feature Layers is sharing the same very long legend that shows all the color references for all the different layers.

According to the documentation on legends, it should be pretty simple. My code snippet is small but I really believe I'm missing something here. If people whine for more, I can include more later.

var legend = new Legend({
   map: map
}, "key");
legend.startup();

The legend contents load to a div called "key". From what I can tell before the browser freezes and crashes, is that some legend contents are loading to the div. This is where it gets hung up. It seems like there are too many legend contents in different places or that it cannot locate one or some of the legend contents on the server. Is there something that I'm missing when it comes to loading legends?

This is the error I get...
"Cannot read property '_layers' of undefined"

EDIT After trying to isolate some of the layers to see which ones are causing an issue. I still get the error when I'm trying to load more than 4 legends from different layers. Is there a way to make loading the legends more efficient? How do you only load one legend at a time?

EDIT Here's some more information describing more details of my problem. I've changed the name of my div from "key" to "loadlegend"

I'm trying to get legends loaded from a feature service that contains multiple legends within the same layer. The browser loads forever and sometimes gets some of the legends loaded and then crashes at a certain point. It looks like all the legends are contained in the same layer on the feature service. You can see my development website with all the firebug notices here.

I've also tested other feature services with multiple legends and not had this issue. I'm not sure if it has to do with the shear quantity of legends mine contains or what. I think there are at least 10 legends in the single feature layer.

You can also see where I've got this little test of the sample server on JSFiddle, although the results are a little funny looking. If you copy and paste the code into ArcGIS Sandbox, it should work though. Sidenote, why does the legend look different when I go to the reference link?


texts: How to use reclassify tool in ModelBuilder?

texts:
I have a model which buffers line features by a variable and then rasterizes them. I then need to reclassify before I am finished processing.

My issue is that I cannot access the values to reclassify in the model because they have not been created yet...  I have tired to use the add field button in the tool and generalize the reclassification I need (very simple - nodata is 0 and all values 0 - 100 are 1 so I get 0 or 1 in my raster...


texts: Formula for charting a course around the globe, given a location and degree heading?

texts:
I'm looking for a formula that takes a long/lat coordinate and a degree heading, (e.g. "I'm standing in San Francisco and I'm facing 45º from north"), and which can draw a line across a map projecting your path along a globe from that point, along that heading, and back around to you.

So, for example...

If you were standing at the equator facing due west, it would just be a horizontal line along the equator.
If you were facing due north, it would be a vertical line from your position, and another vertical line at the longitude exactly halfway around the world.
If you are standing at a random point and facing a random heading, it would create one of those "sin-wave" paths along the map, like we often see on maps of space stations or shuttles, starting and ending at your current location.

Does a formula like that exist?  It seems like it would have to, since it would be one of the most basic formulas for calculating orbital paths and such.  But I haven't been able to find anything like it...


texts: How do you set the MS_MAPFILE environment variable in IIS?

texts:
I'm trying to set up MapServer as a WMS on a windows server running IIS. All works well until I try to cascade the service through GeoServer and the getcapabilities URL is flagged as invalid as it includes a full path to the map file in it and : is not a valid character.

In the short term I can work around this by using a relative path (and / instead of \) but I'm not happy having external URLs depend on my local file layout, so I looked at hiding the map= section of the URL.

This is usually done by setting the MS_MAPFILE variable, but I can't work out how to make IIS see it. I even tried setting it in the system environment but it seems to make no difference.

In the good old days when I last tried this I would have wrapped the .exe in a .bat file and set the variable there to make sure, but IIS won't let me run anything but a .exe (or .dll).

So is there any way to set an environment variable in IIS so that MapServer can see it.


texts:
The area of interest is 1450 km2.

My options are:

50cm newly collected aerial natural color imagery for $28,300

OR

50cm recent archival satellite imagery (includes near IR band) for $20,600

The primary use will be land appraisal, geology(rock outcrop identification), and general GIS work such as digitizing of hydrology and transportation features, orientation maps, etc.

The company hired to fly the area and process the imagery is having various technical difficulties so our product has not yet been delivered.

Recent satellite imagery is available for the entire area. I am aware of the pan-sharpening multispectral imagery process and that aerial image (when available) will be show more detail on the ground but what other pros and cons are there?

Which type of imagery is better for the uses I listed above?