texts: How to push newly added or updated domain values to ArcGIS rest service?

texts:
I'm currently using ArcGIS Server 10.4 for my REST service and ArcGIS Online. I have created a domain list and assigned it to a field. All of it works fine. I'm currently using Oracle for my geodatabase. When i added the the rest service to a webmap, i can see the domains as a select list.

However, from time to time, additional values must be added to the domain or need to be updated. I tried updating it in ArcMap and it reflects nicely. In editing mode you can see it straightaway in the select list. This is not the same for the REST service. After adding in the value, it is still not available when i try to check it via webmap.

I then proceed to restart the REST service and now it starts to show the new entries. Is this the standard way? in order for it to reflect the new entries, i would constantly need to restart the service?

Is there a way to refresh the connection to the geodatabase maybe without the need of fully restarting a rest service. I can see that if you add/or edit data in ArcMap and save editing, it can straightaway reflect. But this is not the case for domain.

Is there a way around restarting the service?


texts: I want to write a game that will draw real geographical map. What software systems can often make requests for free? I need to stay at drawing a map.


texts: Difficulty rendering a multipolygon from a GeoJSON using JQuery

texts:
I'm pulling polygon data (neighborhood boundaries in San Francisco) from an API on https://data.sfgov.org/Geographic-Locations-and-Boundaries/SF-Find-Neighborhoods/pty2-tcw4/data. While I made sure on my console that I'm indeed pulling the data, I am having trouble visualizing the layer as displayed on the original page.

I kept the JS fiddle simple. I've tried different D3 methods but I keep running into dead ends, usually with a feature being undefined.


texts: Making mass/bulk update by attribute value in ArcGIS Desktop?

texts:
I have feature class with 70,000+ records (file geodatabase).
I want to make in bulk update on an attribute where the value is NULL.
Select by attribute shows that I have 5,000 records like that. Doing it one by one is going to be a pain.

In sql pseudo code it should be something like that:

 UPDATE <MyAttributeName> = 'NewValue'
 WHERE <MyAttributeName> IS NULL

How do I do that with ArcMap 10.1?


texts: What happened to the changelog on QGIS 3 plugin info page?

texts:
In QGIS 2.X the changelog from the metadata.txt was/is displayed on the plugins info page (and besides, the info whether a plugin is trusted or not):

In QGIS 3 I miss both informations:

Users are justifiably interested in such changelogs, and often do not want to browse to some metadata.txt.

Do I miss a setting or is this intented? Or a bug?


texts: Get UTM zone with GeoPandas from a shapefile

texts:
I want to know to which UTM zone belong a shapefile.

With GeoPandas I can read a shapefile and know its EPSG code:

import geopandas as gpd
dataframe = gpd.read_file('shapefile.shp')
print(dataframe.crs)

    {'init': 'epsg:32616'}

Is there anyway to get the UTM zone with GeoPandas or to transform the EPSG code to UTM zone? for example epsg:32616 would be the UTM zone 16 N.


texts: Displaying two separate Excel files on same map?

texts:
I have started to learn ArcMap.

I have two separate Excel files which are:
Border coordinates,
CSD coordinates.

For each Excel files, I have information for longtitude and latitude.
My aim is to see both border coordinates and CSD coordinates on same map in ArcMap. To see clearly which points are border, I want to make different colour for border coordinates than CSD coordinates.


texts: Reloading OpenLayers layer from GeoServer when underlying data changes?

texts:
I have a webpage that displays a GeoTIFF served via GeoServer as an OpenLayers3 ImageWMS layer.
What I need is to (via a button) reload the OL3 layer after the underlying GeoTIFF changes.

I searched around, and tried

layer.dispatchChangeEvent();

on the layer in question. But the map/image isn't updated.

I have tried

1) removing the layer, constructing a new layer from the GeoTIFF and add it back using JavaScript.

2) disabling caching for the raster layer in GeoServer.

3) using layer.getSource().dispatchChangeEvent(); as suggested by @ahocevar's answer.

None solved the problem.

The only thing that works is to reload the whole webpage, in which case the map shows the new image but all states in the page is lost.

How can I reload the layer without having to reload the whole page?


texts: Landsat 8 cloud detection without band 9

texts: What is the best possible way to create a cloud mask of a Landsat 8 scene without using band 9? I think the best way will be to determine the cloud cover by thresholding in the Panchromatic band and verify the output with the Thermal band (clouds are cold...) - but I am open for any other ways.


texts: Set property (date) of images in Imagecollection from value in date band in Earth Engine

texts:
I have a Imagecollection with 2 bands per image. One is the 'date' band. The pixel values from this band are equal in every image. The value of one of the pixels i want to use to add as date to the properties of the image. In order to use filter.Date(). i dont understand the error because '20141018' is a string and 'yMd' is the right format to transform or isn't it? http://joda-time.sourceforge.net/apidocs/org/joda/time/format/DateTimeFormat.html

var new_collection = vvAscending
    .map(function(image) {
      return ee.Feature(null, {'date': ee.Date.parse('yMd',image.reduceRegion(ee.Reducer.first(),geometry,10).get('date'))})
    })
    .distinct('date')
    .aggregate_array('date')
print("new_collection ",new_collection )

The Error i get is:

ComputedObject (Error)
Collection.reduceColumns: Error in map(ID=0):
Date.parse: Date/time '20141018' cannot be parsed with format 'yMd'.

Line 107:113
https://code.earthengine.google.com/6afaa18508be5b108581056077efe6ea


texts: I've downloaded the latest QGIS version and I cannot add basemap. I installed OpenLayers plugin and while adding any map from list QGIS crashes. I tried to enable Google Mercator CRS first, yet it's not in the list of available CRSs. Does anybody know how to fix this and how can start adding basemaps? Or is it possible to add National Geographic map as a basemap?


texts: How to import multiple shapefiles without any dbf to a PostGIS database?

texts:
I have a polyline shapefile (.shp and .shx) without .dbf. I would like to import the shp into PostGIS.

In order to do that, I created the corresponding .dbf using the r code provided by mdsumner in this post  Adding Attribute data to shapefile?

The .dbf consists of just one row, while the .shp corresponds to a polyline.

Now, if I load the shp into QGIS, the attribute table lists the features ids of the polyline in read-only mode.
I know that I could save the file loaded in QGIS then I would be able to import it into PostGIS.

However I have thousands of similar shapefiles and I would like to import all of them into PostGIS. I suppose that this is not possible because the .dbf is read-only because I got this message: "dbf file (.dbf) can not be opened".

I read in another post, that it is the shx file that connects the shp to the dbf but this did not help me.

Apologies if the question is not relevant or it is not the correct way to proceed.


texts:
I wrongly deleted Data Management Tools from ArcToolbox in ArcGIS 9.3. Now I can't use some features such as "Dissolve", "Eliminate", etc.

How could I make it as default?


texts:
When I try to make 3D scene from DEM with Qgis2threejs it creates an empty HTML page. I also try to make 3D thematic map from shapefile with the same result.

I have also check my browser and it supports WebGL. I have use this plugin last year on the same computer last year, but since I updated the plugin I can't get it to work. Any suggestion would be great.


texts: Returning raster values from multiple raster tables via Lat/Lon in PostGIS?

texts:
I am new to PostGIS and am unsure if I'm thinking of this correctly, but I'd like to return values from multiple raster tables in a PostGIS database in one go.

The purpose is to query legend values for a webmap based on the lat/long values of a cursor click. So, if there are multiple raster tables in the database e.g., raster1 and raster2, I would like to shoot lat and lon to the DB and return the values from both rasters in some sort of dictionary.

How would I go about doing so with an SQL command? Would I need to do some sort of spatial join? Multiple select statements?


texts: Visualizing multitemporal rasters scenes, side-by-side in QGIS?

texts:
I have to classify land cover/change every 2 months. Every 2 months about 5000 random selected points are classified, using L7/L8 images. For this, images from 5 different dates are show, side-by-side. Its enough if a can see a 'minitile' of the region around the selected point (1:30.000 scale).

I can do this using ENVI software, but I want to know if I can do this in QGIS.
There is a plugin for this or I'll have to figure one out?


texts: Adding text to end of existing value in QGIS

texts:
I have a column in a large attribute table. The values in this column are numeric value.

Example entries in this case would be:

1, 2.0512, 8.5744, 8.944

How can I add the word 'm' to the end of each entry so it'd look like:

1m, 2.0512m, 8.5744m, 8.944m

I don't want to use the technique of sorting them, then selecting them, then using field calculator to define them. I just want to know how to add text at the end of existing value, no matter what the existing value is.


texts: What is the unit of the area when geometry attributes are calculated in QGIS?

texts:
My CRS is WGS 84. I intersected two files containing polygons. Each file has an attribute table stating the area of each of the polygons (in m2). After intersecting them I navigate to: Vector >> Geometry Tools >> Add Geometry Attributes. The area of the new polygons are calculated but they are orders of magnitude different to the original polygons. For instance, when comparing a new (smaller) polygon to the original (larger) polygon (which are only slightly different), the new area is calculated as 28.08 and the original is 243,943,129,649 m2

There is no way that these can be the same units. How can i find out the units of the new polygons?


texts: Dissolved area of polygons using the gArea function of rgeos?

texts:
I am taking occurrence point data (in the form of individual Lat/Long points), and adding a 1km buffer around them using gBuffer() in R package rgeos.

I'm then using the function gArea() to calculate the total area of the generated shapefiles. However, the result is the sum of each polygon's area, such that overlapping polygons are counted multiple times.

What I am hoping to get is only the total area of their footprint (essentially equivalent to a dissolved shapefile in ArcGIS).

Because I need to do this process for over 120 species, I'm looking for an automated approach in R rather than manually dissolving the shapefiles in ArcGIS.

# Generate 1km point buffers. Projected CRS for function: USA Contiguous albers equal area
geodata_Albers <- spTransform(dat_sp, crs.albers)

# Create 1 km buffer and add the total area to the AOOs dataframe.
sp_buffer_1km <- gBuffer(geodata_Albers, width = 1, byid= TRUE)

AOOs[i,]$buffer_1km <- gArea(sp_buffer_1km, byid = FALSE)


texts:
I am using node to capture gps data from my devices after every 30 seconds, i want to store the data into a table but when executing the query statement it throws an error

Failed.........error: bind message supplies 3 parameters, but prepared statement "" requires 1, i am using prepared statements. My geometry field is a PointField, database Postgresql with postgis. any one having a clue.

My query code:
`

client.query("INSERT INTO gpsdevices_coordinate (gpsdevice, geometry) VALUES ($1, ST_GeomFromText(\'POINT($2 $3)\', 4326))", [this.getUID(), data.longitude, data.latitude], (err, res) => {
console.log("Failed........." + err);
 //Do something where done
        })

`


texts:
Here I have two lines in QGIS. There is a part when they overlap. I want to merge them so it will be one line in the middle.

I have tried the solutions already mentioned previously but no success.
What I tried:

answer to Merging close lines using v.clean
How to merge lines with slightly different endpoints in QGIS?
Merging lines (with and with no overlaps) in QGIS
Merging overlapping lines into one line. QGIS
Join / merge lines that touch and intersect in QGIS
Merging adjacent lines in QGIS
How to connect two line ends in QGIS without resorting to other software?

This is what I want to get:

How do I merge them?


texts: Proper way to check if a Factory Code is Projected or Geographic

texts:
I would like to create a SpatialReference with ArcObjects, but I won't know whether it is Projected or Geographic, I'll only have the SRID.  Is it fine practice to use 'ISpatialReferenceFactory.CreateGeographicCoordinateSystem' with a 'try', see if it fails, and if so, use 'CreateProjectedCoordinateSystem', and if that fails, give a message that it won't work?  I have a feeling that I shouldn't use Exceptions in this way, but I'm not aware of an alternative.

I am aware of these answers, but they don't quite answer the question:

Arcobjects: Is my spatial reference Geographic or Projected?

How to programmatically determine whether a spatial reference factory code is valid?

Edit: This is after the fact, but I have become aware of this answer


texts: How do I prevent my QGIS Plugin from opening twice?

texts:
I have made a new QGIS Plugin, almost everything is working fine, but I can click on my Plugin-Button several times and it starts several times. So I see my Plugin 4 or 5 times. How can I prevent that?

Almost no Plugin out there shows this behavior, but I can't figure out what I do have to change.


texts: On hover show feature information of features inside buffer ol5

texts:
I have created a buffer when user enter radius. Points are selected inside buffer but I dont know how to show feature information of selected features.

Here is my code:

 var select = new ol.interaction.Select({
    style: pointStyle
 });

 map.addInteraction(select);
 var selectedFeatures = select.getFeatures();

 var feature ;
 map.getViewport().addEventListener("dblclick", function(e) {
// circleLayer.getSource().clear();
 // select.getFeatures().clear();
 if(feature!=null){
        map.removeLayer(featureRequest);
        vectorSource.removeFeature(feature);
        marker.removeFeature(featureMarker);
        select.getFeatures().clear();
 }

    var coordinate = map.getEventCoordinate(e);
    console.log("coordinate = "+coordinate);
    //add marker
         var latLong = ol.proj.transform(coordinate, 'EPSG:3857', 'EPSG:4326');
        featureMarker = new ol.Feature(new ol.geom.Point(coordinate));
        marker.addFeature(featureMarker);
    //add marker
     feature = new ol.Feature(new ol.geom.Circle(coordinate, dist));

//GET CIRCLE EXTENT
    var ext = feature.getGeometry().getExtent();

    //SELECT FEATURES BEFORE ADDING CIRCLE TO PREVENT THE CIRCLE FROM BEING SELECTED
   vectorSource.forEachFeatureIntersectingExtent(ext, function(features) {
        selectedFeatures.push(features);
    });

    //ADD CIRCLE
    vectorSource.addFeature(feature);

    /*** jsts ***/
                var geom_select = feature.getGeometry();
                console.log("geom_select = "+geom_select);
                var pt = reader_wkt.read('POINT (' + geom_select.getCenter().join(' ') + ')')
                var jsts_geom_select = pt.buffer(geom_select.getRadius());
                console.log("jsts_geom_select = "+jsts_geom_select);

    /*** jsts ***/
featureRequest = new ol.layer.Vector({
 title: 'Filtered point Map',
  visible: true,
  source: new ol.source.Vector({
    loader: function(extent) {
      $.ajax('http://localhost:8080/geoserver/wfs', {
        type: 'POST',
        data: {
          service: 'WFS',
          version: '1.1.0',
          request: 'GetFeature',
          typename: 'export:Ground_Control_Point',
         outputFormat: 'application/json',
          srsname: 'EPSG:3857',
          bbox: extent.join(',') + ',EPSG:3857'
        }
      }).done(function(response) {
            console.log("response "+response);
             features = new ol.format.GeoJSON().readFeatures(response);

             // extent_swe=vectorSource.getExtent();
            features = features.filter(function(el) {
                  if (jsts_geom_select.contains(parser.read(el.getGeometry()))) {
                    return true;
                  }
                });
                console.log("features "+features);
                 select.getFeatures().extend(features);
                // vectorSource.addFeatures(features);
      });
    },
    strategy: ol.loadingstrategy.bbox,
    projection: 'EPSG:3857'
  }),
    style: new ol.style.Style({
          fill: new ol.style.Fill({
            color: 'rgba(255, 255, 255, 0.2)'
          }),
          stroke: new ol.style.Stroke({
            color: '#ffcc33',
            width: 2
          }),
          image: new ol.style.Circle({
            radius: 7,
            fill: new ol.style.Fill({
              color: 'aqua'
            })
          })
        })
});

 map.addLayer(featureRequest);
 });

How to get information on hover of selected points?


texts: Calculate distance in point collection algorithm

texts:
I've got collection of points which are sorted by x,y and distance between n and n+1 (array of objects [{x:0,y:0},{x:1,y:1}....] I would like to calculate distance between all points (point 1 distance to all points, point 2 distance to all points and so on) and add weight to point depends of points in distance parameter (for ex. 500m). Can someone help me with algorithm, because below algorithm time is O(n^2) so maybe there is faster way than that..

_groupPoints : function(points) {
    var points = points.sort(function(a, b) {
        return Math.sqrt(Math.pow((a.y - b.y), 2) + Math.pow((a.x - b.x), 2));
    }), tmp = [], idx = 0, length = points.length, i = 0, j = 1, res = 0;

    tmp.push([points[0]]);
    idx = 0;
    for (i; i < length; i++) {
        j = i + 1;
        for (j; j < length; j++) {
            res = Math.sqrt(Math.pow((points[i].x - points[j].x), 2) + Math.pow((points[i].x - points[j].x), 2));
            if (Math.floor((((res * Math.PI) / 180) * 40076000) / 400) <= 500) {
                tmp[idx].push(points[j]);
            }
        }
        idx++;
        tmp[idx] = [points[j]];
    }
}


texts: Manipulating ASCII (.asc) files in QGIS?

texts:
I have an .asc file and I want to import in QGIS.

The standard way to open it is from "add Raster layer" ("import Raster Layer") but is there another way?

I need to convert the .asc file to .tif or to another raster format with QGIS.


texts: Merging polygons in shapefile based on similar attributes

texts: I have parcel data for several counties that contain landowner information.


texts: pyqgis 3.0 : how to project raster layer, avoid pop-up prompts interface

texts:
rasterlayer = QgsRasterLayer(self.infile, layername)
raster_crs = rasterlayer.dataProvider().crs()

if raster_crs is None:
     self.rasterlayer.setCrs(QgsCoordinateReferenceSystem(4326,
 QgsCoordinateReferenceSystem.EpsgCrsId))

how to avoid the pop-up prompt interface for projection, the above code is useless, what should i do?


texts:
Here's an extract of a feature from the web that I converted to geojson using qgis. However, the coordinates are not standard are at all and when trying to actually preview it on the web it fails, however, qgis has no problem reading it. My question is how do I transform this set of coordinates using QGis into "standard" coordinates. Data comes from here http://donnees.ville.montreal.qc.ca/dataset/affectation-du-sol/resource/c38ea3f0-fa70-48c9-9421-92113c581b69

{
"type": "FeatureCollection",
"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:EPSG::900913" } },

"features": [
 { "type": "Feature", "properties": { "categorie": "rural" }, "geometry": { "type":     "MultiPolygon", "coordinates": [ [ [ [ -8232754.718677375465631, 5691987.542297831736505, 0.0 ], [ -8232826.554151621647179, 5691863.061614270322025, 0.0 ]

Here is the prj file

PROJCS["Ontario_MTM_Zone_8_east_of_75_degrees_W_NAD_83_datum",GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["Geodetic_Reference_System_of_1980",6378137,298.2572221008916]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-73.5],PARAMETER["scale_factor",0.9999],PARAMETER["false_easting",304800],PARAMETER["false_northing",0],UNIT["Meter",1]]


texts: Spectral profile (z-profile) viewer for ArcGIS

texts:
I have multilayer raster with 6 bands. Is there a way to select one specific pixel and to display the z profile (the values of each band) for this pixel as a chart or table? I am looking for a tool similar to the spectral profile viewer of ENVI: link

As a workaround I tried to cast my raster pixels into points but there are way too many of them for a smooth display.


texts:
I'm attempting to crop a raster with gdal warp using python bindings, the code works fine as I use it to crop a raster approximately 3k times, of which a couple fail with the following error and produce no outputs:

gdal warp cutline invalid

So following some online suggestions I added the following line of code:

gdal.SetConfigOption("GDALWARP_IGNORE_BAD_CUTLINE", "YES")

which then changes the error I get to the following:

gdal Points of LinearRing do not form a closed linestring

So from this I can see that the polygons that I generated from a linestring then buffered using ogr aren't correct. I've then gone and extracted the geometry as a WkT so I can check them, I've included the full list at this URL: https://pastebin.com/d60L3nni .

Now I'm still stuck because the first and last point are at exactly the same coordinates:

-1406502.73228006 6411785.18363632

and

-1406502.73228006 6411785.18363632

Sure because the first and last points are on the same coordinates then it should be a valid polygon?
I've tried passing this throught ArcMap Repair Geometry and it also passes the geometry.IsValid() check in ogr.


texts: Python code to count business days in ArcGIS field calculator gives ERROR 999999?

texts:
I'm trying to field calculate using a python script the number of business days between two given dates.  Right now the dates are hard coded, but once I get the script working I will pass dates from two other fields to it.  The script works fine in IDLE, giving me the correct number of business days, but when I put it in the field calculator it errors out:

  ERROR 999999: Error executing function

  The field is not nullable.

The field is set as "Short" with Precision at 5.

Pre-Logic Script Code:

import datetime as dt
import numpy as np

def workdays():
    start = dt.date(2018, 1, 31)
    end = dt.date(2018, 2, 12)
    days = np.busday_count(start, end)
    return days

Field=

workdays()


texts: Labelling polygon with multiple values from Attribute table

texts:
What is a good way to put a list of names in one polygon feature?

For example, I'm creating a County Facilities layer. There are cases when you have multiple departments in the same building.

I don't want separate polygons for these departments, but be able to label them from the same polygon.


texts:
I am trying to add layers from this ArcGIS MapServer to my Leaflet webmap:
http://igems.doi.gov/arcgis/rest/services/igems_haz/MapServer

This code added the weather warnings layer (12), but used the default blue symbol for all of the feature types.

var URL2 = 'http://igems.doi.gov/arcgis/rest/services/igems_haz/MapServer/12/'

var warnings = L.esri.featureLayer({
    url: URL2,
    useCors: false
    });

warnings.addTo(weather);
weather.addTo(map);

Looking at this service description page, it seems that there is are preset symbols configured on the server:
http://igems.doi.gov/arcgis/rest/services/igems_haz/MapServer/12

From example for an Ash Warning:

Symbol:
Style: esriSFSSolid
Color: [169, 169, 169, 255]
Outline:
  Style: esriSLSSolid
  Color: [0, 0, 0, 0]
Width: 0

Is there an option for L.esri.featureLayer to use the preset symobology?


texts: Adding a table to a web map made in Folium

texts:
I have made a nice Folium map however, the amount of data is small, but increasing. Right now it is grabbing a wildfire perimeter GeoJSON from NIFC and then getting points from my organization and doing a spatial join for the points within the fire perimeter. Then I filter the joined points from the original list and displaying one as a normal points layer and the other as a MarkerCluster because there are thousands of points.
It would be great to display a table of the spatially joined points and polygons in two separate tables on the side or bottom of the map. It doesn't necessarily need to be interactive (although that would be cool), I think a static HTML table would work fine for now. My googling has lead me to the to_html from the Pandas library that I'm already using with Folium. I haven't found any examples or ways to inject/pass that HTML into the HTML file for the map when Folium generates it.
Should I try and manipulate the HTML file separately after Folium generates it?


texts: qgis 2.8.1 arcgis filegdb incomplete list in win-version

texts:
i am running qgis 2.8.1 (OSGeo4W 64bit-version) on mac (10.10.3) and win (7)

... and was just recently confronted with different behavior of both versions when trying to open (ESRIs) file.gdb.
the following img shows the FeatureClass listing in ArcCatalog (see the two polygon FCs at the end)

this img shows the "add vector layer / directory / ESRI FileGDB" result on mac:

and tho following one the corresponding result on Win7

on win7 one table and the point featureclass are displayed properly
.

Why are the polygon feature classes missing?


texts: R spatial: add new attribute to S4 class object?

texts:
I am trying to check if my shapefile polygons have open edge or not. To have an open edge, the polygon either has no neighbors, or is >5 m taller than its neighbours. So, I want to add this new field open_edge TRUE or FALSE to my attribute table.

My idea was to loop through each polygon row by row (similar to arcpy.da.searchCursor in ArcGIS), and depending of the number of adjacent polygon fill in attribute value, or move to next row.

This is how I can fill in all values in a field:

polys@data$open_edge <- TRUE

But within for loop i.e. subsetting by one row at time polys[i, ]@data$open_edge <- TRUE I got an error:

Error in as.vector(data) :
  no method for coercing this S4 class to a vector

How can I insert new value to the attribute table row by row?

My code:

polys= readOGR(getwd(),
                layer = "polys")

# Crete new vector to store open_edge value
polys@data$open_edge <- NA

for (i in seq_along(polys)) {
  one = polys[i, ]
  left = polys[-i,]

  # Identify the number of neighbors for each polygon
  nbrs <- left[which(gTouches(sp::geometry(one),
                              sp::geometry(left),
                              byid = TRUE)),]

  if (nrow(nbrs) == 0) {
    polys[i, ]@data$open_edge <- TRUE  ## how to fill this???
    print("in if statement")
  } else {
    print ("something else")
    # compare here the tree heights of stands
    # ...
  }
}

I expect to get back filled attribute table in new filed:

polys@data$open_edge
TRUE TRUE FALSE TRUE TRUE FALSE...


texts:
I am attempting to use the Esri resource-proxy to proxy all WMS requests to other servers. In doing so I am attempting to get around a CORs issue that does not require a user browser plugin. I can get the proxy working for any site that does not require user authentication with a PKI certificate, but when I attempt to access a site that requires a cert, the request fails and throws a "PKIX path building failed" error. It appears as though the user's certificate is not being passed from the proxy to the end site.
Is there a way to pass the user's PKI credential through the proxy to the desired destination?
I currently have a Tomcat instance configured to use PKI for user authentication (this part is working) on my local machine, but am not having success passing that credential to the end site.
I have noticed the proxy config can be set to use different authentication methods (access token, username/password, etc..) but all of these seem to be preconfigured and not dynamic per user.
I am not too familiar with PKI and fairly new to web development.


texts: QGIS Download File Algorithm Ignore SSL Errors

texts:
I am requesting QGIS to download a file using Download File Algorithm in a pyqgis script. The URL is as follows https://www.mrlc.gov/geoserver/mrlc_display/NLCD_2016_Impervious_L48/ows?version=1.3.0&service=WMS&layers=NLCD_2016_Impervious_L48&styles&crs=EPSG:3857&format=image/geotiff&request=GetMap&width=421&height=261&BBOX=-11698584.467247624,4925812.478408471,-11685967.801238675,4933631.516008289&
Note
This URL is working on off, but the issue can still be tested on https://wrong.host.badssl.com/

The QGIS is throwing up this error and once I save SSL exception by doing Save&Ignore option. This window will not pop up but the algorithm will still fail.

I checked trough the browser and there is a certificate issue on the website but I can ignore it and proceed forward.

How can I request qgis file downloader to ignore this error? Here is my code:
    alg_params = {
    "URL": request_URL,
    "OUTPUT": QgsProcessing.TEMPORARY_OUTPUT,
}
outputs["DownloadNlcdImp"] = processing.run(
    "native:filedownloader",
    alg_params,
    context=context,
    feedback=feedback,
    is_child_algorithm=True,
)


texts: Data-defined styles in Mapbox vector tiles

texts:
Is it possible to style Mapbox vector tiles (Mapbox GL JS) using a feature's attributes? Assuming each feature has an attribute color which is a hex color string, I tried

"paint": {
    "fill-color": function(d) { return d.color }
}

but received an error saying that a string was expected, not a function. Looking through some examples it seems that features of similar styles needed to be grouped into their own layers, but that seems a bit unwieldy in the event of a single layer with hundreds of possible styles.

What is the best way to style polygons based on a feature's attributes?


texts: The table was not found error - ArcGIS 9.3.1, Python

texts:
I have a Python script that does the following on multiple tables

maketableView - not shown in code
Addjoin
Copyrows
MakeTableView

it was working OK, until all of a sudden it stopped.

Here's part of my code:

gp.AddJoin_management("HDATA_L_LICENCE_view", "HIRE_COMPANY", "HDATA_L_HIRE_COMPANY_view", "HIRE_COMPANY", "KEEP_COMMON")
gp.copyrows_management("HDATA_L_LICENCE_view", myTempWorkspace + os.sep + "HDATA_L_LICENCE_0")
gp.maketableview_management(myTempWorkspace + os.sep + "HDATA_L_LICENCE_0","HDATA_L_LICENCE_0_view")

gp.AddJoin_management("HDATA_L_LICENCE_0_view", "HDATA_L_LICENCE_LICENCE_TYPE", "HDATA_L_LICENCE_TYPE_view", "LICENCE_TYPE", "KEEP_COMMON")
gp.copyrows_management("HDATA_L_LICENCE_0_view", myTempWorkspace + os.sep + "HDATA_L_LICENCE_1")
gp.maketableview_management(myTempWorkspace + os.sep + "HDATA_L_LICENCE_1","HDATA_L_LICENCE_1_view")

gp.AddJoin_management("HDATA_L_LICENCE_1_view", "HDATA_L_LICENCE_0_HDATA_L_LICENCE_REGSITE", "GAZETTEER_layer", "NSG_REF", "KEEP_COMMON")
**gp.copyrows_management("HDATA_L_LICENCE_1_view", myTempWorkspace + os.sep + "HDATA_L_LICENCE_2")**
gp.maketableview_management(myTempWorkspace + os.sep + "HDATA_L_LICENCE_2","HDATA_L_LICENCE_2_view")

The first two blocks work, but when the code tries to run on the line

    gp.copyrows_management("HDATA_L_LICENCE_1_view", myTempWorkspace + os.sep + "HDATA_L_LICENCE_2")

I get the following error:

    Error in creating output D:\ScheduledTasks\Temp\TempWorkspace_HyMS_Licences.gdb\HDATA_L_LICENCE_2
The name of the Field is invalid: valid names may contain letters, numbers or underscores.
The table was not found. [HDATA_L_LICENCE_2]
The table was not found. [GDB_Locators]
Failed to execute (CopyRows).

Has anyone got any ideas, as I've looked at the code for a few hours but I can't seem to figure out why it is not working... and my head hurts :(


texts:
I would like to display a Message Box Dialog when the user activates the -esri selection tool- and double clicks anywhere on the map. or when they doubleclick any feature on the map

please note that I am not creating a new tool. MapControl is an option but I am not able to access the object, or cast anything to it. so please help

EDIT:-
I have tried to cast my Focus map into a Mapcontrol and get that object and cast into ImapControlEvents2. and threw it I can handle the OndoubleClick event

I have tried getting the selection tool from the active tool and caste it into an Itool object but that didnt work since I am not creating a new tool,

I researched the API Reference provided on esri website to see what options I have and I found nothing that might help me. the interfaces that I need are the following,


texts:
How can I tell mapserver to ignore a filter that is set in a mapfile?

I am displaying a wms layer with point features using mapserver. In one of the layers in my mapfile I have defined a filter similar to this:

FILTER          "mycolumn='%myvalue%'"

such that I can use runtime substitution in a mapserver request similar to this one:

http://someurl.com/mapserv?map=mymap.map&layer=mylayer&mode=map&myvalue=thisvalue

This works just fine.

However, sometimes I would like to display the entire dataset and thus forcing mapserver to ignore the filter.

Is it possible to deactivate the filter in a request?


texts: Where can I find shapefile (or other vector) data for US private parcel boundaries (such as for a plat map)?


texts: Convert Google Earth Engine images to histograms without downloading them first

texts:
I want to download a large collection of MODIS satellite data. The problem is that the data will be upwards of 200GB and this is not feasible for me to download. So I want to access the images and convert them to histogram of pixel counts instead without having to download the images. This way I only store the histogram representation of the images.

How can I access the images and convert them on the fly without having to export them to Google Drive or Google Cloud?


texts: Geoserver REST API: Duplicate layer name

texts:
Using the web admin page, I created two layers in different workspace but with the same name:

In the admin view I see something as this

Workspace  Store        Layer Name
WK1        PostGIS1     mylayer
WK2        PostGIS2     mylayer

All works very well, except when I try to access the properties of these layers with Geoserver REST API.

/layers.xml give me something like:

<layer><name>mylayer</name><atom:link rel="alternate" href="http://localhost:8080/geoserver/rest/layers/mylayer.xml" type="application/xml"/></layer>
<layer><name>mylayer</name><atom:link rel="alternate" href="http://localhost:8080/geoserver/rest/layers/mylayer.xml" type="application/xml"/></layer>

Two identical entries! And when I ask for the layer with /layers/mylayer.xml I only get the properties of the first...

How can I access the properties of the second layer with the REST API? Can it be done or it's a limitation or a bug of the REST API?


texts: Is it possible to load 2 base maps in Qfield (both of the 2 rasters are in geopackage format), and choose wich one I want to see and maybe switch during the field work ? I tried with the map themes but it does'nt seems to work with rasters.


texts: My raster values are always changed to zero

texts: I'm uploading a series of XYZ data points onto GIS and the co-ordinates are coming up fine, but some of the Z values are changed to zero - and in some cases, all! I've interpolated the majority of the data, but cannot for those with 0's across the board. Trawled these pages for answers but I just can't get my head around it. Would someone please put me out of my misery?


texts: How to filter (clip) a WFS before download by Shapefile?

texts: I have one polygon stored in a shapefile and an WFS Layer with a lot of Polygons. Is it possible to filter the WFS by the shapefile-polygon before the Download begins, in QGIS or ArcGIS? I just want to Download these polygons which are within the big area (.shp).


texts:
for a science project at university I have the task to create a web application based on Geoserver. The data is already stored, I have to visualize it now on the web for different users. Each user can only see certain layers.
I only have basic programming skills but experience in Drupal.

Would Cartaro be a good solution for this?


texts: vertexInfo properties in ArcGIS JS 4x API for event when reshaping polygon

texts:
When dragging a polygon vertex using the Sketch class, the update event is fired, with a type of "reshape".
4X:
sketch.on("update", onGraphicUpdate);

function onGraphicUpdate(event) {
      if (event.toolEventInfo && event.toolEventInfo.type === "reshape") {
        console.warn(toolEventInfo);
        //this does not seem to have any info about the vertex being dragged, even though we have the
        //polygon
      }
}

Previously in 3x, we could grab a vertexInfo object off of that type of event (vertex-move), which contained the pointIndex and segmentIndex of that vertex.
3X:
editToolbar.on('vertex-move',function(evt){
   console.log(evt.vertexInfo.pointIndex);
});

Does anyone know if pointIndex and segmentIndex exist anywhere, or how to decipher/identify which vertex within the polygon is actively being moved?


texts: Correct SQL expression when looping different fields

texts:
I'm having trouble with the follow lines of code. Practically what I am doing is a 10 shapefile loop in which field represents the name of the field that I want to be calculated( the name differs by every shapefile and the name of that column is also the name of the file). Since it's 10 shapefiles, I want each one to be calculated individually to the new field I created. Everything works fine except the SQL expression. I have no idea how to format the expression in so that field is divided by the [MinorityS] field. If I were to be doing it without the loop it would just be the expression of '[Cuban] / [MinorityS]'. I just want Cuban to be the 10 shapefiles which in this case is field representing all the field names that I want to be in that first bracket of Cuban.

import arcpy
arcpy.env.workspace = r"myworkspace"
arcpy.env.overwriteOutput = True
MetroCT =  r"metrotracts"

Fieldslist = arcpy.ListFields(MetroCT)
Metrofields = [f.name for f in arcpy.ListFields(MetroCT)]
Clusterslist = ["SouthAmeri", "SpanishSpa", "Mexican", "PuertoRica", "Cuban", "Dominican", "CentralAme", "European", "Arab", "Africans"]

for field in Clusterslist:
    out = field + ".shp"
    Cluster = arcpy.ClustersOutliers_stats(MetroCT, field, out)
    output = field + "1"+".shp"
    HHClusters = arcpy.FeatureClassToFeatureClass_conversion(Cluster,r"mylocation", output, "\"COType\" = 'HH'")
    arcpy.AddField_management(HHClusters, "MinorityS", "LONG", 10)
    cursor = arcpy.da.SearchCursor(Cluster, field)
    summed_total = 0
    for row in cursor:
        summed_total = summed_total + row[0]
        arcpy.CalculateField_management(HHClusters, "MinorityS", summed_total)
    arcpy.AddField_management(HHClusters, "MinbyTot", "LONG", 10)
    expressionmin = expressionmin = '[{}] / [MinorityS]'.format(field)
    arcpy.CalculateField_management(HHClusters, "MinbyTot", expressionmin, "VB")


texts: Does anyone know why there are multi-year gaps in the Landsat record in some places? This screenshot is from earthexplorer, but I found the whole when my Earth Engine code wouldn't run because it couldnt' find imagery. I was skeptical but, lo and behold, there is a >6 year gap in coverage of this path/row (and those around it that are part of my Earth Engine analysis). I had a similar error in an area of Brazil that I suspect also has a multi-year coverage gap.


texts: Update existing field using CASE -statement

texts:
When I run this statement in qgis fieldcalc it returns correct reults only from the first line ('Skog') and the rest gets NULL-values.

Any ideas what might be wrong here? Is it not pocssible to run CASE-statement inlcuding several lines? It seems to work when I run it with the three middle rows excluded.

CASE
when "KATEGORI"  LIKE '%skog' then 'Skog'
when "KATEGORI"  LIKE '%vatten%' then 'Vattenomrade'
when "KATEGORI"  LIKE '%Annan%' then 'Oppen mark'
when "KATEGORI"  LIKE '%ker%' then 'Oppen mark'
Else 'Bebyggd och Hardgjord mark'
END


texts: Issues with projection system in R (raster package)

texts:
I'm new with GIS in R, and I'm facing a problem I really do not understand...

The aim is easy : I have dataset of GPS points (decimal) located on a Mediterranean island, and I want to extract Bioclim data for each point.

  My dataset looks like this

  loc

    code    municipality        x        y status balea
1   AHO1         Alghero 8.243606 40.61025      m   bal
2   AHO2         Alghero 8.175231 40.62431      m   bic
3   AHO3         Alghero 8.236722 40.61369      m   bal
4   AHO4         Alghero 8.236989 40.62351      m   bal
5   AHO5         Alghero 8.242231 40.61663      m   bal
6   AHO6         Alghero 8.239211 40.61161      m   bal
7   AHO7         Alghero 8.249450 40.61243      m   bal
8   ASI1    Porto Torres 8.305914 41.09242      m   bal
9   MAC1         Macomer 8.705692 40.24535      m   bal
10   MAD  Santulussurgiu 8.619711 40.16222      m   bal
11   MON       Montresta 8.526969 40.37630      m   bal

First, I create a map of Sardinia using function "map"

zone.etude <-map("world", regions="Italy:Sardinia", exact=TRUE,fill=T)
IDs <- sapply(strsplit(zone.etude$names, ":"), function(x) x[1])
zone.etudeSpp <- map2SpatialPolygons(zone.etude, IDs=IDs, proj4string=CRS("+proj=longlat +datum=WGS84"))
e <- extent(zone.etudeSpp)

Then, I create a raster stack with bicolor variables

for (i in 1:19) {
  r<-raster(paste("bio_",i,sep=""), proj4string=CRS("+proj=longlat +datum=WGS84"))
  r.crop=crop(r, e)
  r.mask=mask(r.crop,zone.etudeSpp)
  assign(x=paste("r.crop_",i,sep=""),value=r.mask)
}
predictors=stack(r.crop_1,r.crop_2,r.crop_3,r.crop_4,r.crop_5,r.crop_6,r.crop_7,r.crop_8,r.crop_9,r.crop_10,r.crop_11,r.crop_12,r.crop_13,r.crop_14,r.crop_15,r.crop_16,r.crop_17,r.crop_18,r.crop_19)
projection(predictors) <- "+proj=longlat +datum=WGS84"

The issue appears now : when I use extract to get data for GPS points, the first 8 locations get NA...

val.bioclim = extract(predictors, loc[,3:4],cellnumbers=TRUE)

> val.bioclim
      cells bio_1 bio_2 bio_3 bio_4 bio_5 bio_6 bio_7 bio_8 bio_9 bio_10 bio_11
 [1,]   587    NA    NA    NA    NA    NA    NA    NA    NA    NA     NA     NA
 [2,]   586    NA    NA    NA    NA    NA    NA    NA    NA    NA     NA     NA
 [3,]   587    NA    NA    NA    NA    NA    NA    NA    NA    NA     NA     NA

Exactly as if they would "fall into the sea"... [I know 2 of them are out of the raster extent, as they come from small islets not included in Bioclim. But the remaining 6 are definitely far from shore]

So I decided to first check both projections

> proj4string(predictors)
[1] "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"
> proj4string(zone.etudeSpp)
[1] "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"

They seem quite coherent, as they are all expressed in longlat and WGS84.

Then I made a map using plot, and that's where things went wrong

plot(predictors$bio_1)
plot(zone.etudeSpp, add = T)
points(loc[,3:4], pch = 16, cex = 0.6)

If I check the plot in my Rstudio window, everything looks OK, all point fall "inland" within the extent of the "predictor" raster, and the edge of both layer (predictors$bio1 and zone.etudeSpp are superposed)... but when I export
the map in PDF, my points and zone.etudeSpp seems to be in the same projection, but are completely stretched toward N and S... which is not the case for the predictor plot. Thus I can see my 8 first points "falling into the sea".

[I'm struggling to add picture but I've screen shot if this to make it clearer... just don't know how to add them]

I'm quite lost,

Is it a problem with CRS? I thought both maps and my GPS locations were expressed in the same system
Is it a problem with my computer, or R version? I have updated all my packages this afternoon.

Guillaume


texts:
I am building a cadastral fabric and it's really nice to have the bearing and distance come up on the lines as I traverse around. I just found a great way to do this with the Add Geometry ArcGIS tool. It works most excellent:

But the issue I am having here is that the field does not auto-populate. I have rerun the Add Geometry tool everytime. Is there anyway to automate calculation? Its a pain to have to traverse a line and then run the tool.

Using the answer of @FelixIP, I was able to get it going, just add this into the label expression, set the parser to python and customise your attributes:

arcpy.AddGeometryAttributes_management(Input_Features="Cadastral Reference Lines", Geometry_Properties="LENGTH;LINE_BEARING", Length_Unit="FEET_US", Area_Unit="SQUARE_METERS", Coordinate_System="PROJCS['NAD_1983_UTM_Zone_17N',GEOGCS['GCS_North_American_1983',DATUM['D_North_American_1983',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Transverse_Mercator'],PARAMETER['False_Easting',500000.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',-81.0],PARAMETER['Scale_Factor',0.9996],PARAMETER['Latitude_Of_Origin',0.0],UNIT['Meter',1.0]]")
def FindLabel ([BEARING], [LENGTH]):
  return [BEARING] + " - " + [LENGTH] + "ft"


texts: Reading elevation data from WMS response?

texts:
I am trying to get elevation data from the tile received from USGS server.

I query the data using Image Format image/Png. When I read the data I get color values.
I need to know how to read the elevation values from it.

Is there any documentation or sample available?


texts:
I have to create multiple maps using the same layout, but different map views (more than 100 maps).  There is a separate rectangles layer created for each view at a specific scale using the map frame from the layout window. Within the layout window, I would like to automate the process by updating automatically attribute information for layers displayed within the map frame within text boxes on the layout.

Is there any possibility to do this process through a MapBasic script? As far I have researched the topic, I didn't find any solution. Also, I know that Mapinfo 15.2 has the capability to support multiple pages layouts, but the amount of information that needs to be displayed increases the resources usage, which make the program run really really slow.


texts: Obtaining areas from overlapping polygons

texts:
Using ArcView 9.
I have a data set of about ten thousand parcel polygons and a similar number of driveway polygons extending from a roadway polygon into the parcel polygons. My ultimate goal is to get the driveway area within each parcel.
Ideally this could be done through the creation of a new shapefile so that I can have a concise plan for each parcel, but any method that works is welcome.


texts: Doing definition query on layer and zooming to it using arcpy.mapping?

texts:
I am trying to do definition query to a layer and to zoom to the same layer, but I got an error ("SyntaxError: can't assign to operator") ??

dataFrame = arcpy.mapping.ListDataFrames(mxd)[0] # the first data frame
for lyr in arcpy.mapping.ListLayers(mxd,"lyr", dataFrame)[0]:
        Current_Layer=lyr.name
        if lyr.name=="Polygons":
            lyr.definitionQuery = "[Floor] ="+ Floor_num  and "[Physical_number_1] = "= physical_num
            dataFrame.extent=lyr.getSelectedExtent()


texts: How to use osmfilter to drop features with non-empty values on a key?

texts:
I know one can use --drop="key=value" option with osmfilter to filter out OSM features for which key has a value of value.

For example, I can use

osmfilter input.o5m --drop="service=parking_aisle" --drop="service=driveway"  -o=output.o5m

to drop out service roads with the service tag and parking_aisle or driveway values.

How can drop all features for which the service have a non-empty value?

I tried --drop="service=", and it seems to drop features with an empty service tag too.

Any official documentation on how to do this?


texts:
I have successful created a shapefile from an island with contours by downloading the data from http://srtm.csi.cgiar.org, converted it with gdal and imported it in TileMill.

Everything is working very well, except that the contour lines don't stop at the coast. It looks like this coast from OpenTopoMap https://opentopomap.org/#map=17/28.52831/-16.13842

Is there any way to cut or stop the contour lines at the coast?


texts:
Link to GEE

function BAI(Im){
var pcr = 0.1
var pir = Im.select('B3')
var pcnir = 0.06
var pinir = Im.select('B4')
var one = 1
var bai = one.divide(pcr.subtract(pir).pow(2).add(pcnir.subtract(pinir).pow(2)))
return bai
}
var BAIColle = LSatTMCf.map(BAI)

When running this code, the Error "pcr.subtract is not a function". I guess the issue is with pcr beeing a number and pir being a Pixel Value from an Image


texts: Why do I get invalid elements (polygons), when smoothing a vector layer (smooth geometry)?

texts:
The problem occurs, when using the commands "Smooth geometry" or "v.generalize.smooth" in QGIS.

The shapefile, which I want to smooth, shows a flood water level and looks like this:

The water level contains some hard edges, which wouldn't look like this in reality, so I want to smooth those edges a little bit, using the commands I mentioned above. It gives a nice looking result, but unfortunately includes invalid elements.

I uploaded the files at wetransfer: https://we.tl/45WV0KpTLx

water_level_HQ30_diff.shp: original file without validity problems
water_level_HQ30_smooth.shp: smoothed with v.generalize.smooth (Grass Gis 7 command)
water_level_HQ30_smooth_native.shp: smoothed with "Smooth geometry" (QGIS geoalgorithm)

How can I get this vector layer smoothed without getting invalid elements?


texts: How do I prevent QGIS error messages from vanishing?

texts:
I get error messages with some layers, but I can only see their first line and they disappear before I get the chance to read them.
How can I prevent this issue?


texts: NNJoin, MMQGIS and Nearest Neighbor QGIS failed to measure shortest distance

texts: I have two layers, one point and one line.  I need the distance to the closest line from the point layer.  Not NNJoin, MMQGIS nor the distance hub in QGIS provided accurate measures.  Any ideas?


texts: Fill NoData raster cells with floating point value of closest cell

texts:
I am working in ArcMap 10.3.1 with Spatial Analyst tools.  I have a raster with some large contiguous areas with NoData.  I want to assign the floating point value of the closest cell to each NoData cell.  This would be similar (as I understand it) to using the Euclidean allocation tool, but the return value (i.e., output raster) for that tool is integer type, not the floating point type which I need.  I want to do this only for areas with NoData, w/o affecting values in cells surrounding the NoData areas. The attached image shows a small part of the study area (colored areas have values, areas w/ DEM background have NoData). A similar question was asked here, but again, I want results to be floating point.
Any ideas?


texts: MapServer won't open multiple layer with STATUS on

texts:
So I made this mapfile and I tried to load multiple layers with status set to on but that doesn't work.  The layers with "STATUS DEFAULT" load but the other two won't when I change the URL.

#### places LAYER
  LAYER   NAME  "places"
          TYPE           polygon
          DATA           "places"
          STATUS         DEFAULT
            CLASS
                  STYLE
                          COLOR           "#b2b266"
                          OUTLINECOLOR    "#331a00"
                  END
          END
  END
  #### water LAYER
  LAYER  NAME "water"
         TYPE         polygon
         DATA         "water"
         STATUS       ON
         CLASS
                STYLE
                        COLOR               "#1a8cff"
                        OUTLINECOLOR        "#00264d"
                END
        END
  END
  #### railroads LAYER
  LAYER  NAME "railways"
          TYPE             LINE
          DATA            "railways"
          STATUS        ON
          CLASS
                  STYLE
                          COLOR              "#ff6600"
                          SIZE                    3
                          SYMBOL              "line7"
                  END
          END
  END

If I set "LAYERS=water" or "LAYERS=railways" the layers will load but "LAYERS=water,railways" won't with the URL:

http://localhost/cgi-bin/mapserv.exe?map=/ms4w/apps/MyProject/project.map&layers=water,railways&mode=map


texts: Plotting approximate geospatial data on historical map?

texts:
Our library is undertaking a project to map historical data regarding 18th and 19th century British and American printers and publishers. One of the problems that we have encountered is the different granularity of location data. For example, one publisher may list a street number, street name, and city whereas another publisher may only mention a street name or even only the city. In such cases, we are wondering the best way to georeference the location and present this approximate data.

Does the spatial humanities community have rules or conventions for dealing with these cases?

Also, what are the most effective ways for visualizing “approximate vs precise” geospatial data on the same map?


texts: Split GEE FeatureCollection according to latitude

texts:
I have an ee.FeatureCollection with points across the globe. As a pre-processing step, I need to split it according to whether the points are in the Northern or Southern hemisphere.
Conceptually it's not a big deal, but since I am new to GEE and its JavaScript syntax, I can't get it to work.
Consider we have the FeatureCollection myPoints:
var myPoints = ee.FeatureCollection([
  ee.Feature(ee.Geometry.Point(-60,-20), {'id': 0}),
  ee.Feature(ee.Geometry.Point(140,-20), {'id': 1}),
  ee.Feature(ee.Geometry.Point(-70,10), {'id': 2}),
  ee.Feature(ee.Geometry.Point(80,40), {'id': 3}),
  ee.Feature(ee.Geometry.Point(10,50), {'id': 4}),
  ee.Feature(ee.Geometry.Point(20, -10), {'id': 5})
]);

I think it should work something like:

Define a function "GetLatitude" which gets the latitude of each point.
I thought of myPoints.geometry() but I don't succeed at indexing the resulting list of lists
myPoints.geometry()[0][1] is not the way to go.
Define a function "FilterHemisphere" that serves as a handle to split positive latitudes to the
Northern hemisphere, and negative latitudes to Southern hemisphere. JavaScript has a function
Math.sign() to tell positive numbers from negative numbers, but I don't see any equivalent under
Docs->ee.Number.
Apply the functions to myPoints. I thought in proper GEE JavaScript it should be something like
myPoints.map(GetLatitude).filter(FilterHemisphere)

I wonder if anyone could help me with defining the functions and how to apply them in proper JavaScript?


texts: Find and Replace python script to replaces values in attribute table?

texts: I am working on a stand alone script in python that will allow to me to specify a shapefile and being able to find a value and replace it with a new value. Does anyone know what code will allow me to do this?


texts: Erasing line features based on intersecting polygons using GeoPandas?

texts:
I have two interesecting geopandas.GeoDataFrame datasets, one containing line features (lines_gdf) and one containing polygons (polygons_gdf):

import matplotlib.pyplot as plt
import geopandas as gpd
from shapely.geometry import Polygon, LineString

# Create sample polygon and line data
polygons_gdf = gpd.GeoDataFrame(geometry=[Polygon(zip([150, 151, 151, 150],
                                                      [30, 30, 31, 31]))])
lines_gdf = gpd.GeoDataFrame(geometry=[LineString(zip([148, 150, 152],
                                                      [28, 30, 32]))])

# Plot example data
polygons_gdf.plot()
lines_gdf.plot(ax=plt.gca(), edgecolor='red')

Desired output

I want to use my polygons dataset to remove or erase any part of my lines dataset it intersects with (essentially a GIS "erase" operation), and return only the non-overlapping parts of my lines (i.e. black lines below):

Using the geopandas.clip tool I've been able to do the exact opposite of what I need to do: clip the data to keep (not remove) the intersecting portion of my lines:

# Clip lines to keep only area within `polygons_gdf`
clipped_gdf = gpd.clip(lines_gdf, polygons_gdf)

# Plot clipped lines over the top of input layers
polygons_gdf.plot(alpha=0.1)
lines_gdf.plot(ax=plt.gca(), edgecolor='red', alpha=0.1)
clipped_gdf.plot(ax=plt.gca(), edgecolor='black')

My question: How do I use geopandas to remove/erase part of a line dataset that intersects with a polygon dataset?

My geopandas version is v0.6.3.


texts: Get a mask for MODIS (250m, MOD09GQ) using MODIS (500m, MOD09GA) in Google Earth Engine?

texts:
How can I use the quality information stored in MODIS (500m) to create a mask for my MODIS (250m) image - e.g. cloud mask? When I merged the two image collections and run a function to detect the cloudy pixels, I got an error:

Error in map(ID=1_2018_01_01):
Image.select: Pattern 'state_1km' did not match any bands. My code:

// A function to mask out cloudy pixels.
var maskClouds = function(image) {
  // Select the QA band.
  var QA = image.select('state_1km')
  // Make a mask to get bit 10, the internal_cloud_algorithm_flag bit.
  var bitMask = 1 << 10;
  // Return an image masking out cloudy areas.
  return image.updateMask(QA.bitwiseAnd(bitMask).eq(0))
}
// Get a feature collection with Mexico boundary
var countries = ee.FeatureCollection('ft:1tdSwUL7MVpOauSgRzqVTOwdfy17KDbw-1d9omPw')
var Mexico = countries.filter(ee.Filter.eq('Country', 'Mexico'));
var dataset = ee.ImageCollection('MODIS/006/MOD09GQ')
                  .filter(ee.Filter.date('2018-01-01', '2018-05-01'));
var terra = ee.ImageCollection("MODIS/006/MOD09GA");
var test = dataset.merge(terra)
.filterDate("2018-01-01", "2018-05-01");
var image = test
.filterBounds(Vietnam)
.map(maskClouds);


texts:
I have 300 m resolution  ESA Land cover raster that I want to clip using
the shapefile of Asia. (The shapefile has a border for each country, so therefore not dissolved).

Now I ran this modified code in the clipper:

gdalwarp -q -cutline /data/Forest_Extent_Asia_pacific_Ocean_Threats/asian_countries.shp -tr 0.00277777777778 0.00277777777778  -dstalpha -crop_to_cutline
-of GTiff /data/Forest_Extent_Asia_pacific_Ocean_Threats/glob_cover.tif /data/Forest_Extent_Asia_pacific_Ocean_Threats/glob_cover_v2.tif

My question is, why does it takes too much time to finish the gdal processing? Is there a bug that caused it?


texts: ArcMap: Output to XML Including Shape Data?

texts:
I need to create a solution for non-technical people to export data from an ArcMAP document (.mxd) to XML files. I will need each feature's data to be exported in such a way that the output xml file (or could be an Excel file) contains both spatial and non-spatial attributes.

The only I put I see which has the required data is a Google Earth .kmz file. Now, I have been able to convert the .kmz to .kml using some .NET code but the final .kml file is too unstructured to be parsed by my .NET code.

I have tried xToolsPro and they do directly output to .kml file but that doesn't good either--and costs $250 for the license (the export function will expire after the trial).

So what are my options? Could I program some Python script and create some extension to Arctoolbox?

I am stuck! Can't believe this kind of functionally is not built-in (export to xml including Shape data).


texts: How does one modify selected geometry with a function in qgis python console?

texts:
My example use-case is that I want to redistribute vertices along lines as in:

Redistribute linestring vertices to specified distance
How to create points in a specified distance along the line in QGIS？
How to automatically add vertices to lines every x meters?
How to densify a polyline according to a specific length?

I've written some functions to do some of the work:

import math

def redistribute_vertices_num(line,num):
    return LineString([line.interpolate(i/float(num),normalized=True) for i in range(num+1)])

def redistribute_vertices_dist(line,dist):
    return LineString([line.interpolate(i*float(dist)) for i in range(int(math.ceil(line.length/dist)))])

I would like to apply these functions (or other, modified functions) to the selected objects in the current layer to modify the layer.  I'm looking for a general function like this:

edit_layer_selection_geometry(redistribute_vertices_dist,500)
# or
def my_function(line):
    redistribute_vertices_dist(line,500)
edit_layer_selection_geometry(my_function)

...that might do something like the processing in How can I create a line with three points with Python in QGIS? or How can I switch line direction in QGIS? on the active layer with the user-defined function.

def edit_layer_selection_geometry(myfun, extra_arguments...):
    layer = qgis.utils.iface.mapCanvas().currentLayer()
    # check for selected_layer, editability ?
    for feature in layer.selectedFeatures():
        newgeom=myfun(feature,extra_arguments)
        # update feature with geometry
        layer.changeGeometry(feature.id(),newgeom)
    ...
    # update layer stuff (extents?)

I'm looking for the edit_layer_selection_geometry() or an equivalent functionality, and I'm happy to modify the potential 'myfun' functions to fit requirements of the general function.

I've put this in my .qgis2/mycode.py directory:

from qgis.core import (QgsFeature, QgsGeometry,
                       QgsVectorLayer, QgsMapLayerRegistry,
                       QgsField)
from PyQt4.QtCore import QVariant
from qgis.utils import iface
import math

def redistribute_geometry_vertices(geom,dist):
    mylen=geom.length()
    line=geom.asPolyline()
    return QgsGeometry.fromPolyline([geom.interpolate(i*float(dist)).asPoint() for i in range(int(math.ceil(mylen/dist)))])

def redistribute_selection(dist):
    """Redistribute the vertices along a polyline """
    layer = iface.mapCanvas().currentLayer()
    for feature in layer.selectedFeatures():
        geom = feature.geometry()
        newgeom=redistribute_geometry_vertices(geom,dist)
        layer.changeGeometry(feature.id(),newgeom)

Then I can solve the use-case by selecting the layer and feature, opening the python console and typing:

import mycode
mycode.redistribute_selection(0.1)


texts: Validating polygons for no gaps/overlaps/slivers using ArcGIS Desktop?

texts:
I would like to assess the spatial accuracy of a load of polygon shapefiles which are supposed to show air quality for a bunch of countries.

Firstly I want to assess whether each shapefile has a complete coverage of the country it is supposed to represent.

Secondly I want to assess whether there are any slithers in the shapefiles, areas which fall out of the country boundary.

My initial method is using the rasterize polygon method, using the extract by mask tool to exract the supplied data by a known country boundary shapefile, and perhaps counting the pixels in the original file, and the extracted file, for each case.

Can you think of a more efficient method for doing this?

I tried to use the split function but it didn't work as I expected it would.


texts: QGIS: changing style not possible any more

texts:
The setting: I encountered a very strange behaviour using QGIS on Win 10, standalone installer. It both happens with 3.16 and 3.18 (I also tested with older, still installed versions like 3.4: same problem).
The problem: I was able to work normally with both QGIS versions in the last days and weeks. But today, suddenly, I was not able any more to change layer styling - see screenshot at the bottom. Whatever I tried, no way to set another style, QGIS always sticks to the styles attributed by default when creating new layers or by the one defined in the project when opening existing projects. Assigning one of the pre-defined styles does not work either. The preview in the layer styling panel however works correctly.
What I tried:

I restarted QGIS several times and my computer as well. When I finally had a look to the QGIS version installed, I remarked that it was QGIS 3.18.0 instead of 3.18.1. I had installed 3.18.1 before, but for some mysterious reason it seems to have gone. The same, by the way, with 3.16.5 - it's gone, when I start 3.16, the pointrelease 3.16.3 (the last I installed before 3.16.5) opens. In both versions, I'm not able to change styles.

I re-installed QGIS 3.18.1, the newst standalone release, downloaded directly from https://www.qgis.org. However, when I start the newly installed version, the problem persists.
Screenshot: Installation folder with the newly installed 3.18.1 version - you can see as creation date today's datum:

The context: The problem happened after QGIS crashed on the Processing plugin was deactivated. I re-activated it. Don't know if that has some connection to the problem.
The Question: Does anyone have an idea what happened and how to solve the problem?
Screenshot: opening an empty project, adding a temporary polygon layer in EPSG:3857, drawing a single polygon with valid extent. As you can see on the right in the layer styling panel, changing the style does not affect how the polygon is rendered. Adding additional symbol layers also does not change anything. Same for lines and points:


texts: Is there a service or library that takes a route (with unlimited points) as input and returns text turn-by-turn instructions?

texts:
I am making an application that generates routes based on atypical heuristics so I would like a service much like Mapbox Directions API but I could:

Disable the internal routing. (i.e. if I give a path A -> B -> C I don't want the api to change my route to the shortest/fastest path of A -> C much less waste the resources to do so).
Since i'll be inserting the entire route I would need the capability of inputting a virtually limitless number of points instead of the 2 - 26 that are offered by most directions apis.

Essentially I am looking for a turn-by-turn instruction generator that is decoupled from the routing portion of the library.

The closest thing I can find would be to set up my own OSRM backend as that would allow me to achieve my second goal of unlimited waypoints but still runs its own routing and is too much hassle to set up and maintain.

Does anyone have any ideas of a library that may achieve what I'm looking for?

My question is similar to some existing questions like  Seeking route calculation service to get route between several points? but the main difference being I do not want any routing so please don't refer me to libraries that have the routing portion like OSRM.

For example, I have a list of coords that represent the exact route that I have calculated. I give the api these coords.

        "coordinates": [
      [
        -87.62252748012543,
        41.88330112135344
      ],
      [
        -87.62416362762451,
        41.88329313366847
      ],
      [
        -87.62420117855072,
        41.880840867169994
      ],
      [
        -87.6208108663559,
        41.88088480111109
      ],
      [
        -87.62087523937225,
        41.881004620796936
      ],
      [
        -87.62073040008545,
        41.88429558033205
      ],
      [
        -87.62248456478119,
        41.884275611428535
      ],
      [
        -87.62250602245331,
        41.883297127511085
      ]
    ]

And the api would return directions without changing the route.

       "steps": [
        {
          "points": [
            {
              "location": [
                87.426579,
                41.508068
              ]
            },
            {
              "location": [
                -87.6206579,
                 41.88308068
              ]
          ],
          "geometry": "mnn_Ick}pAHUlAqDNa@",
          "mode": "walking",
          "duration": 18.6,
          "maneuver": {
            "bearing_after": 125,
            "type": "depart",
            "modifier": "right",
            "bearing_before": 0,
            "location": [
              -87.6206579,
              41.88308068
            ],
            "instruction": "Head southeast on Köpenicker Straße (L 1066)"
          },
          "ref": "L 1066",
          "weight": 18.6,
          "distance": 98.1,
          "name": "The Bean"
        },
        {
        .....
    }
  ],


texts: Adding Point Graphic to ArcGIS API for JavaScript Map?

texts:
Can you please take a look at this JSFiddle and let me know why I am not able to add the Marker into the Map? I am getting this error

  Uncaught TypeError: Cannot read property 'add' of null

Here is the code I have

var map;
var graphicsArray = [];
require(["esri/map",
    "esri/geometry/Geometry",
    "esri/geometry/Point",
    "esri/geometry/Polyline",
    "esri/geometry/Polygon",
    "esri/graphic",
    "esri/symbols/SimpleMarkerSymbol",
    "esri/symbols/SimpleLineSymbol",
    "esri/symbols/SimpleFillSymbol",
    "esri/Color",
    "esri/InfoTemplate",
    "dojo/domReady!",
    "esri/geometry"], function (Map,

Geometry,
Point,
Polyline,
Polygon,
Graphic,
SimpleMarkerSymbol,
SimpleLineSymbol,
SimpleFillSymbol,
Color,
InfoTemplate) {
    map = new Map("map", {
        basemap: "topo",
        center: [-106.61, 35.1107],
        zoom: 13
    });

    var point = new Point(-106.61, 35.1107);
    var pointSymbol = new SimpleMarkerSymbol();
    var pointAttributes = { city: "Albuquerque", state: "New Mexico" };
    var pointInfoTemplate = new InfoTemplate("Albuquerque");
    var pointGraphic = new Graphic(point, pointSymbol, pointAttributes).setInfoTemplate(pointInfoTemplate);

    graphicsArray.push(pointGraphic);
    for (i = 0; i < graphicsArray.length; ++i) {
        map.graphics.add(graphicsArray[i]);

    }

});

      html, body, #map {
          height: 100%;
          width: 100%;
          margin: 0;
          padding: 0;
      }
      body {
          background-color: #FFF;
          overflow: hidden;
          font-family:"Trebuchet MS";
      }

<link rel="stylesheet" href="http://js.arcgis.com/3.14/esri/css/esri.css">


texts: QGIS, how do you break segments of a polygon, ie break 1 polygon into 2 or more?

texts: I have a polygon but would like it to be in broken up to put into different layers. I can't work out how to do this.


texts:
Is there a way to continuously rotate an icon from the CalciteWebCoreIcons?

I have a layer with the single icon on it, that shows up correctly.

    view.graphics.add({
      id: "centers",
      symbol: {
        type: "text",
        color: getRandomColor(),
        text: "\ue613",
        angle: 180,
        font: {
          size: 30,
          family: "CalciteWebCoreIcons"
        }
      },
      geometry: {
        type: "point",
        latitude: parseFloat(_lat),
        longitude: parseFloat(_lon)
      },
      popupTemplate: {
        title: "",
        content: ""
      }
    });

Basically, i would like the icon to rotate 360 degrees every 5 seconds, without clearing and reloading the entire layer.


texts: How to remove/disable last row/column of SRTM GeoTIFF in QGIS?

texts:
I've started using QGIS yesterday, and one of my main reasons to use it is to apply colormapped SRTM v4 over Google aerial imagery.

Everything is working fine, except that SRTM v4 elevation tiles have overlapping last row and last column (6001 x 6001 pixels), and that creates duplicated raster over the edges, which means an ugly, "unacceptable" border over the displayed image.

So before I consider editing the files to remove the last row and column, I wonder if there is a way to tell QGIS to ignore this part of the GeoTiff.

The image below (usind soft light blending mode) shows the problem over the latitude of -30°


texts: Multipolygon supposed to be contained by another but returned False with Shapely

texts:
There are two multi-polygons a and A where a is supposed to be contained inside A. Using shapely, the following difference operator
b = a.difference(A)
print(b)

gives
POLYGON EMPTY

However testing A.contains(a) returns False; a.touches(A) also gives False.
How could this be possible with Shapely?
The wkt descriptions for A and a are hosted in the following Gist link (due to characters limit)
https://gist.github.com/vodp/169b12f507a9730c0dd2ac86ca2589db
Edits
So what I am trying to test is assessing the equivalence between
a.difference(A) -> False

and
A.contains(a) -> True

But A.contains(a) returns False instead. If I were buffering a bit, for example
A.contains(a.buffer(-1))

then I get True. Notice that both
aa = a.buffer(-1)
aa.difference(A) -> EMPTY

and
a.difference(A) -> EMPTY

Does that mean for contains to work, none of the point of a lying in the boundary of A ?
Any chance this is because of rounding errors? I try to get the minimal buffering value and the smallest one I could get before .contains() return False is that:
A.contains(a.buffer(-0.000000001)) # 1e-9


texts: weighting scheme of gstat for fitting models to variograms

texts: I use R and the gstat package to construct variograms of temperature, and I would like to ask If anyone knows what weighting scheme does the gstat package use to fit a model to the sample variogram?


texts: User query results surpass "Maximum Number of Records Returned by Server"

texts:
I have a feature server running on ArcGIS Server 10.61. The MaxRecordCount parameter is set to 500 for the feature server. Nonetheless, with normal 'User' role credentials, I am able to select and export well in excess of 500 records from individual layers when using ArcGIS Pro/Desktop 10.8/QGIS 3.8 to connect to the feature server. The expected limitation is observed when connecting to the data with ArcGIS Desktop 10.2.
Is there a secondary parameter that could be causing this behavior in the newer software?


texts: MapInfo Tab Files - Which feature types are contained?

texts:
Without opening the TAB within a GIS package, is it possible to determine which feature class it contains?

Could the TAB file be open in notepad and assessed for example?


texts: Creating random point polygon fill symbology in QGIS?

texts:
Using QGIS 2.18.3, I would like to apply a random point pattern to my polygon symbology.  The default Style Manager has a regularly-spaced point pattern, called Cemetery-25-50k (please see image below).  This is a good starting point, but obviously not random.

Alternatively, I have reviewed the Point Pattern fill options within the Layer Properties (below), but nothing appears that would generate a random pattern.

What I need to recreate is the ESRI Arcmap random point symbology, shown here:


texts: Pass variable value from Geojson file to filter another Geojson layer Leaflet

texts:
I have an leaflet map with 2 geojson file loaded locally using the leaflet ajax filter.
I am trying to filter the second geojson (polyline) based on 'click' of any of the first geojson point feature. both of them had the same SiteID field
I am able to get the feature ID,

           layer.on("click", function(e){
                 pollid= feature.properties.SiteID; // ID from first layer geojson

                 console.log(pollid);

                lyrSpider.addTo(map); // second geojson file
                console.log(pollid);
            });

how can I pass the (pollid)  to filter the second geojson

 function spiderfilter(feature) {

            if (feature.properties.SiteID == pollid) return true; // problem,
        }

I could not figure out the variable scope. here.


texts: Understanding the differences in PostGIS full version string from two Docker images

texts:
I have to migrate a PostGIS database from mdillon/postgis:11 to postgis/postgis:13-3.1.
From mdillon/postgis:11 docker image I can read:
# select postgis_full_version();

                        postgis_full_version
----------------------------------------------------------------------------
 POSTGIS="2.5.2 r17328" [EXTENSION] PGSQL="110"
 GEOS="3.5.1-CAPI-1.9.1 r4246"
 PROJ="Rel. 4.9.3, 15 August 2016"
 GDAL="GDAL 2.1.2, released 2016/10/24"
 LIBXML="2.9.4"
 LIBJSON="0.12.1"
 LIBPROTOBUF="1.2.1"
 TOPOLOGY
 RASTER
 (1 row)

But from postgis/postgis:13-3.1 docker image I can read:
# select postgis_full_version();

                        postgis_full_version
----------------------------------------------------------------------------
 POSTGIS="3.1.0 5e2af69" [EXTENSION] PGSQL="130"
 GEOS="3.7.1-CAPI-1.11.1 27a5e771"
 PROJ="Rel. 5.2.0, September 15th, 2018"
 LIBXML="2.9.4"
 LIBJSON="0.12.1"
 LIBPROTOBUF="1.3.1"
 WAGYU="0.5.0 (Internal)"
 TOPOLOGY
 (1 row)

So, why the differences?
Especially, why the drop of GDAL, and RASTER in the, yet official, PostGIS image?
Does it mean it doesn't have gdal and raster support?


texts: Missing entities when JOIN between differencce layer to origianl layer in SQL QGIS3

texts:
I have two layers :
1. layer of fake plots
2. layer that shows only the safe area inside those plots (generated with the Difference tool), and contains data about the area of each safe polygon.

this is how the layers look like:

layer no. 1:

Layer 2:

I want to join the information from layer 2 (safe area) to layer 1, so for each plot I  can know the size of the safe area,including plots that are not safe at all (that don't have safe area)
for that I have tried to do left outer join in qgis as following:

SELECT *
FROM "Lease" LEFT OUTER JOIN "safe"
WHERE "Lease"."OBJECTID"="safe"."OBJECTID"

the join worked but in the result table I don't have all my original plots- I have 30 plots but here you can see in the result that some plots are missing( 16,17 and 21):

and the plots that are missing are the ones that don't have "safe" area.

how can I get back those lost plots so they will appear in my table, even if they get null or 0?

My endgoal: to have one table with al the entities from the original plot layer, with the data about the  size of the safe area, even if there is no safe area, so I can claculate the precentage of sade area in each plot.


texts: Is a way to find lost geocode result file?

texts:
I geocoded a list of postal codes into the World Geocode Service and all of them matched properly, with the output file being put as the .gdb. The geocode result was then renamed and I saved it as a layer file. I saved the map previously, but my computer died. I then went to go re-load the program and map, but the shapefile was not under the .gdb, and the layer file was showing up with a missing reference.

Is there any way I can find the lost source files? I looked everywhere on the computer and could only find previous geocoded results in the same database. Would really like to find out as the file was quite large and took up quite a few service credits to geocode.


texts: Fixing sharp edges between GeoTIFF tiles?

texts:
I've downloaded a couple SRTM DEMS from http://srtm.csi.cgiar.org and when I load some of them into QGIS, there is a noticeable edge between 2 of the tiles (see image).

Why does it happen? These are both from the same datasource, why don't they match up? Is there a good way to correct for this in QGIS?


texts: OGR OSM-driver unable to open datasource

texts:
I want to read an OSM file with OGR. I use Version GDAL 1.9.1, released 2012/05/15.

To read a SQLITE file works fine:

$ ogrinfo test.sqlite
INFO: Open of `test.sqlite'
      using driver `SQLite' successful.
1: aaa (Point)

When I try to read a OSM File it fails to open the datasource:

$ ogrinfo test.osm
FAILURE:
Unable to open datasource `test.osm' with the following drivers.
  -> ESRI Shapefile
  -> MapInfo File
  ...

The OSM file itself is not broken (I opened it in QGIS).

What am I doing wrong?


texts: Why does my mosaic not render anything at deeper zoom levels, when the individual rasters are perfectly visible?

texts: We currently generate a mosaic dataset from rasters. We can zoom in as far as we can into the individual rasters. However, nothing is rendered when we attempt to zoom at the same depth on the mosaic. Why is that? What can we do to prevent this behavior?


texts:
I have a DTM with 40cm resolution and a DSM with 2cm. Now I would like to compare the two so I can identify vegetation that is higher than a certain distance from the ground. I changed the pixel size from the DSM to 40 cm and wanted to use the raster calculator with a the simple calculation:
Output = DSM - DTM
I get the message "insufficient memory available for this operation"
I have 16gb of ram and would consider my laptop as good. Even after splitting up the area, it would not calculate the difference.

What can be the issue? Is my equation correct?


texts: How can I make a BaseCustomLayer Draw() fire more often?

texts:
I'm experimenting with the idea of replacing my dynamic display moving graphics with the static (DynamicDisplayEnabled = false) moving graphics.

public class CustomLayer : BaseCustomLayer
{
        public override void Draw(esriDrawPhase drawPhase, IDisplay Display, ITrackCancel trackCancel)
        {
            int hdc = Display.hDC;
            IntPtr p = new IntPtr(hdc);
            Graphics graphics = Graphics.FromHdc(p);
            Image newImage = Image.FromFile(carImage);
            System.Drawing.Point ulCorner = new System.Drawing.Point(x, y);
            graphics.DrawImage(newImage, ulCorner);
        }
}

This code draws the image, but it only invokes draw once, unless I move the map with the pan tool or change the scale.   How can I make the Draw function fire multiple times a second?


texts: R MODISTools downloading from web service currently not working?

texts:
I'm trying to download MODIS FPAR data (and other products) using the R package MODISTools.

It has always worked for me but since about two weeks (beginning of March 2017), I'm always getting the error:

  Downloading from the web service is currently not working. Please try
  again later.

Has anyone found a solution?


texts: Issue with dissolving a polygons layer in QGIS 2.4

texts:
I downloaded a dataset of roads and parks from OSM. I then buffered the roads and merged with the parks.
Now - all i want is to dissolve / union all the features in the merged layer.

I fail to do it using QGIS, even though a test of the same layer in ArcGIS 10.0 worked.

By failing i mean: it takes ages and creates an empty layer at the end.

When trying on a much smaller dataset from the same layer in QGIS it works, so i assume the dataset is valid.

Did anyone else face this issue?

Do you have any idea how to approach it?


texts:
I'm creating a plugin and trying to get the duplicate geometry in my layer and I did it by putting the query in the Python program:
query = """
select * from
(select *, count(geometry) as num_geom from "Mainline" group by geometry) as a
where a.num_geom > 1
"""

vlayer = QgsVectorLayer( "?query={}".format(query), "vlayer", "virtual" )
QgsProject.instance().addMapLayer(vlayer)

This will create a new layer for the duplicate feature.
Instead, I want a code to print out all the duplicate without creating an additional layer. Will also create a GUI using Qt Designer and put the output on a table. Is there a way to still use the query to do that or do you have some code on how to do this without SQL?


texts: Snapping tool issues on QGIS 3.12: Virtual polygons between two edges

texts:
I am digitizing land use inside some cadastral parcels. I have two polygon layers, one for the parcels, and one for land use. Since I don't want to digitize outside the parcels, I used the snapping tool, with the parcel layer checked, snapping to segments and vertices, tolerance 12.

I also checked land use layer, same tolerance, snapping to segments and vertices, and checked 'avoid intersect' to be able to go faster. So I filled out my parcels, making sure to snap each time, etc... At the end though, to check that I had left no holes, I copied my features from my parcel layers, and pasted them in my landuse layer. Since the 'avoid intersect' is checked, it shouldn't be able to copy anything if there is no hole.

However it did copy stuff. When I filter for the new features added, this is what I see, in a lot of places:

Lines, with two or more vertices (red segments. here the cross at the top means there is two vertices at this end). When I zoom in to the maximum, I still see a line, not a polygon. And this line correspond to a shared edge between two adjacent polygons from my land use layer.

Has anyone already experienced this? It is not that bothering for what I intend to do, but for the future, I'd like to know if I'm using the tool wrong, or if it's just a known issue...


texts: Generated GeoTiff in GDAL not opening in ImageMagick

texts:
I have netCDF data I convert to geoTiff using the following code:

import numpy as np
from osgeo import gdal
from osgeo import osr
from netCDF4 import Dataset
import glob

def open_netCDF(filename, path):

    """
    Opens netCDF and returns dataset
    """

    file = Dataset(path + filename + ".nc")

    # get variables
    time = file.variables["time"]
    lon = file.variables["lon"]
    lat = file.variables["lat"]
    data = file.variables["data"]

    # three variables as numpy arrays
    lonarray = np.array(lon)
    latarray = np.array(lat)
    dataarray = np.array(data)

    return dataarray, latarray, lonarray

def write_GeoTiff(data, lats, lons, filename):

    # get bounds
    xmin,ymin,xmax,ymax = [lons.min(),lats.min(),lons.max(),lats.max()]

    # no. rows and columns
    nrows,ncols = np.shape(data)

    # resolution of x/y
    xres = (xmax-xmin)/float(ncols)
    yres = (ymax-ymin)/float(nrows)

    geotransform=(xmin, xres, 0,ymax ,0, -yres)

    output_raster = gdal.GetDriverByName('GTiff').Create(filename + ".tiff",ncols, nrows, 1 ,
                                                         gdal.GDT_CInt16)  # Open the file

    output_raster.SetGeoTransform(geotransform)            # Specify its coordinates
    srs = osr.SpatialReference()                           # Establish its coordinate encoding
    srs.ImportFromEPSG(4326)                               # specifies WGS84 lat lon
    output_raster.SetProjection(srs.ExportToWkt())         # Exports coordinate system
                                                           # to the file
    output_raster.GetRasterBand(1).WriteArray(data)        # Writes my array to the raster
    output_raster.FlushCache()
    del output_raster

if __name__ == '__main__':

    ncfiles = glob.glob("/example/path/*.nc")

    for ncfile in ncfiles:

        ncfile = ncfile.split('/')
        ncfile = ncfile[-1]
        ncfile = ncfile.split('.')
        ncfile = ncfile[0]

        dataarray, latarray, lonarray = open_netCDF(ncfile, "example/path/")

        write_GeoTiff(dataarray, latarray, lonarray, ncfile)

My code runs without complaint. But, when I open in a image viewer (ImageMagick for example) it is black.

I can open with ArcGIS no problem, it works fine.

What is ArcGIS reading/defining that ImageMagick can't from the source GeoTiff?


texts: What are the plugins installed with QGIS 2.18.14, at the beginning?

texts:
I want to organize the plugins to make QGIS open faster and get more memory from PC, but I don´t know what plugins are crucial for QGIS.

The Core plugins seems to be different...

Is there a place to see what are the plugins installed with QGIS at the very beginning?


texts: Raw Data Instead of Zonal Statistics in ArcGIS

texts:
I have 2 layers:

1) Shapefile contains about 15 features.
2) Satellite Imagery.

I would like to get pixel values of satellite imagery (reflectance) in each zone. However unlike 'Zonal Statistics' tool, I do not want statistical summary, rather all the pixel values (So I can further analyze them outside of ArcGIS environment).

Current Approach:
1) Select the feature and create new layer.
2) Clip image using this new layer.
3) Copy from attribute table.
4) Repeat for all features.

This is very time consuming and need to repeat the whole process each time I get new imagery.

Is it possible to automate this process? Specifically, is there a tool, which outputs all the values instead of statical summary?


texts:
I am trying to code save/load function. I can say my save and load function works properly because it displays on console that for example when I load "Canada" file from database, it shows that Canada's properties.color has been changed from default color "green" to red. But it seems like I ran into an issue of not updating. How can this be fixed?
//Map.jsx
import React from 'react';
import { MapContainer, GeoJSON,} from "react-leaflet";
const Map = ({countries}) => {

  const onEachCountry = (country, layer) => {
    console.log(country.properties.color) //This shows me that some of the country colors has been changed
    layer.options.fillColor = country.properties.color; // this should change the color of the map according to country.properties.color. But even after an update in country.properties.color, the color of the map does not change.
  }

  return (
    <div>
          <MapContainer>
              <GeoJSON style = {countryStyle} data = {countries} onEachFeature={onEachCountry}/>{/*Displays the map */}
          </MapContainer>
    </div>
  );
}


texts: Display all field of one record in mapbasic

texts: I have a table with many columns, one of that is lon/lat. After I draw all points on a map, I want to display all other columns contents if I right click on a specific point. (It looks like to change a row to a column). I need some ideas to do this.


texts:
I'm using OpenLayers to draw a polygon to my users. I want the ability to implement through a line, split the polygon into two, with the cutoff line that the user drew.

Does anyone know if this is possible to do with OpenLayers? Maybe indicate some link or sample stuff?


texts: Spatial index on PostGIS with ST_MakeEnvelope testing against point geometries (before vacuum full analyse)

texts:
I just updated this question as I was able to fix the problem by running a vacuum full analyse on the tables. I tested vacuum alone and vacuum analyse before which did not solve the problem. Not sure what vacuum full analyse does but it fixed this kind of problem.

OP>>>

I have several large PostGIS tables with millions of records. Each record has a two geometry fields. A polygon geometry (Polygon, 3857) and a centroid of the polygon (Point, 4326). I also created a (Point, 3857) geometry for testing.

All geometry fields have a spatial index (GIST) and all tables are vacuumed.

I am trying to select records that are within a lat, lon bounding box which I create with ST_MakeEnvelop(). I test ST_Intersect as well as ST_DWithin with a buffer radius of 0. Testing a bounding box against point geometries does never use the spatial index. While testing bounding box against the polygon geometry does use a spatial index. This seems counter intuitive. As I am not even interested in displaying the geometries but only their centroids.

Is there any way to tweak my index or query to get points that are within / intersect a bounding box?

This query uses a spatial index:

SELECT lon, lat, pop
FROM cj_hx_16k
WHERE
  ST_DWithin(
    ST_MakeEnvelope(106.3, 26.4, 130.8, 35.7, 4326),
    geomcntr, 0);

This is the corresponding description for the explain plan.

Node Type = Seq Scan; Relation Name = cj_hx_16k; Alias = cj_hx_16k; Startup Cost = 0.0; Total Cost = 28343.77; Plan Rows = 1127; Plan Width = 23; Filter = ((geomcntr && '0103000020E610000001000000050000003333333333935A406666666666663A403333333333935A409A99999999D941409A999999995960409A99999999D941409A999999995960406666666666663A403333333333935A406666666666663A40'::geometry) AND ('0103000020E610000001000000050000003333333333935A406666666666663A403333333333935A409A99999999D941409A999999995960409A99999999D941409A999999995960406666666666663A403333333333935A406666666666663A40'::geometry && st_expand(geomcntr, '0'::double precision)) AND _st_dwithin('0103000020E610000001000000050000003333333333935A406666666666663A403333333333935A409A99999999D941409A999999995960409A99999999D941409A999999995960406666666666663A403333333333935A406666666666663A40'::geometry, geomcntr, '0'::double precision));

This query does use the spatial index on the geom field.

SELECT lon, lat, pop
FROM cj_hx_32k
WHERE
  ST_DWithin(
    ST_Transform(
        ST_MakeEnvelope(106.3, 26.4, 130.8, 35.7, 4326),
                 3857), geom, 0);

This is the corresponding explain plan:

Node Type = Bitmap Heap Scan; Relation Name = cj_hx_32k; Alias = cj_hx_32k; Startup Cost = 147.26; Total Cost = 2861.79; Plan Rows = 204; Plan Width = 24; Recheck Cond = (geom && '0103000020110F00000100000005000000EDE4E1BBF5916641444475CC81424741EDE4E1BBF59166418189CBFD963F50414311AAACA9C56B418189CBFD963F50414311AAACA9C56B41444475CC81424741EDE4E1BBF5916641444475CC81424741'::geometry); Filter = (('0103000020110F00000100000005000000EDE4E1BBF5916641444475CC81424741EDE4E1BBF59166418189CBFD963F50414311AAACA9C56B418189CBFD963F50414311AAACA9C56B41444475CC81424741EDE4E1BBF5916641444475CC81424741'::geometry && st_expand(geom, '0'::double precision)) AND _st_dwithin('0103000020110F00000100000005000000EDE4E1BBF5916641444475CC81424741EDE4E1BBF59166418189CBFD963F50414311AAACA9C56B418189CBFD963F50414311AAACA9C56B41444475CC81424741EDE4E1BBF5916641444475CC81424741'::geometry, geom, '0'::double precision));
Node Type = Bitmap Index Scan; Parent Relationship = Outer; Index Name = sidx_cj_hx_32k_geom; Startup Cost = 0.0; Total Cost = 147.21; Plan Rows = 3058; Plan Width = 0; Index Cond = (geom && '0103000020110F00000100000005000000EDE4E1BBF5916641444475CC81424741EDE4E1BBF59166418189CBFD963F50414311AAACA9C56B418189CBFD963F50414311AAACA9C56B41444475CC81424741EDE4E1BBF5916641444475CC81424741'::geometry);


texts: Why does the file size increase after I clipped the polygons to be smaller?

texts:
I have a polygon layer that shows properties along a road.

I want to create a kml-file of the polygons, and to downsize the file, I tried clipping the properties with a buffer. But, the filesize of the shapefile increased from 370 kB to 1.97 MB, and when I saved the clipped file to kml it ended up being almost 5 MB?

Why does the filesize increase and what can I do differently to get a smaller kml-file?


texts: Exporting a single feature out of featureCollection as shapefile in Google Earth Engine

texts:
Does any one know, how to export a single feature out of featureCollection in Google Earth Engine?

I am using below code to export every features separately. The fetureCollection has 33 different features which represents different 33 geometry of regions. I want to get those geometry separately as shapefile.

The blow code gives me an error in task section: Error: Invalid argument: 'collection' must be a FeatureCollection. So, tried to get the sub_region(in the code) as ee.FeatureCollection instead of ee.Feature , but I am still getting same error.

var refReg = ee.FeatureCollection("users/username/referenceRegions");

var referenceRegions={
  0:'NTP',1:'STP',2:'ETP',3:'WIO',4:'CAR',5:'NAU',6:'NEB',7:'NEU',8:'SAF',
  9:'SAH',10:'SAU',11:'SSA',12:'TIB',13:'WAF',14:'WAS',15:'WNA',16:'WSA',
  17:'ALA',18:'AMZ',19:'CAM',20:'CAS',21:'CEU',22:'CGI',23:'CNA',24:'EAF',
  25:'EAS',26:'ENA',27:'MED',28:'NAS',29:'ARC',30:'SAS',31:'SEA',32:'ANT'
};

var region_list = refReg.toList(refReg.size());

var reg_size = refReg.size().getInfo();

for (var i = 0; i<reg_size;i++){
   var sub_region = ee.Feature(region_list.get(i));

   Export.table.toDrive({
   collection: sub_region,
   description:referenceRegions[i],
   folder:"GEE_Features",
   fileFormat: "SHP",
   })
}

If you see the below image, the FeatureCollection refReg has all these 33 regions, I want these regions to export individually.

Thanks for the help!! :)


texts: How to properly use JOIN or GROUP By for ST_Union in PostgreSQL with PostGIS extension

texts:
There was already a related question, but it is not answered yet and doesn't provide the necessary information for this case...
I want to cut a raster with a polygon. The raster is stored in tiles in the PostGIS database, as is the polygon.

To find out the overlapping areas I wanted to use ST_Clip or ST_Intersection. Since the raster is stored in tiles and no empty raster should be returned, I want to join the tiles of the raster with ST_Union first and then apply ST_Clip or ST_Intersection.
The problem is that there is not enough information on the PostGIS website about the syntax of these functions or about further necessary steps.
Version 1:
SELECT ST_SetSRID(ST_Clip (ST_Union(landcover_raster.rast), polygon.geom),3035)
FROM landcover_raster, polygon;

Version 2:
WITH raster_together AS( ST_Union (rast) FROM landcover_raster)
SELECT ST_SetSRID(ST_Clip (raster_all, polygon.geom)),3035)
FROM landcover_ausschnitt_polygon;