texts: Generating continuous aspect surfaces from DEM with QGIS

texts:
I'm trying to extract aspect surfaces from different DEM sources at different resolutions in QGIS. I've had success working with SRTM-derived grids (3-arcseconds) and am trying to perform the same gdal-based functions on Canadian DEMs (0.75-arcseconds). My issue is that there are very distinct boundaries between elevation levels and it's causing my aspect calculation to have plateaus.

SRTM Aspect Grid

Canadian DEM Aspect Grid

I'm wondering what I'm doing wrong here or if I'm forgetting something. I'm not sure if I should be smoothing my DEM out or if there's anything the can be done to the Aspect grid to make it continuous.


texts:
I am trying to access the file geodatabase and I get a message of General Function Failure: failed to connect to database. I can access other file geodatasabses, the problem is only with this one. I was doing Network Analysis in this file  geodatabase, in case it's relevant.
I could not find any relevant solutions. Do you have any idea what could be the reason?
I am using ArcGIS Desktop 10.3.


texts:
I need to invoke python plug-in in the QGIS startup. So that the python script will do some functions repeatedly in the periodic intervals until the QGIS software gets closed.

Is there any way to get this done in QGIS?


texts: In Mapinfo it´s easy to link a text-layer with a geometric layer by attributes with outcome as a new layer. Where and how to do this in Q-GIS?


texts: GeoDjango query raster pixel values in PostGIS

texts:
I am building a REST API in GeoDjango that queries pixel values on a multiband raster held in PostGIS (out of db).

I have defined the following model to match the raster held on the database:

from django.contrib.gis.db import models

class Brt2(models.Model):
    rid = models.IntegerField(primary_key=True)
    rast = models.RasterField()

class Meta:
    db_table = 'brt2'
    managed = False

When I try to query the raster model using a point, it's taking around 30 seconds and returning an empty Queryset.

from django.contrib.gis.geos import Point

>> pt = Point(-3.43, 52.946)
>> qs = Brt2.objects.filter(rast__1__contains=pt)
>> qs
>> Traceback (most recent call last):

File "<stdin>", line 1, in <module>
  File "/home.local/share/virtualenvs/AQM_API-RewsMqZ2/lib/python3.6/site-packages/django/db/models/query.py", line 258, in __len__
    self._fetch_all()
  File "/home/.local/share/virtualenvs/AQM_API-RewsMqZ2/lib/python3.6/site-packages/django/db/models/query.py", line 1261, in _fetch_all
    self._result_cache = list(self._iterable_class(self))
  File "/home/.local/share/virtualenvs/AQM_API-RewsMqZ2/lib/python3.6/site-packages/django/db/models/query.py", line 74, in __iter__
    for row in compiler.results_iter(results):
  File "/home/.local/share/virtualenvs/AQM_API-RewsMqZ2/lib/python3.6/site-packages/django/db/models/sql/compiler.py", line 1081, in apply_converters
    value = converter(value, expression, connection)
  File "/home/.local/share/virtualenvs/AQM_API-RewsMqZ2/lib/python3.6/site-packages/django/contrib/gis/db/models/fields.py", line 362, in from_db_value
    return connection.ops.parse_raster(value)
  File "/home/.local/share/virtualenvs/AQM_API-RewsMqZ2/lib/python3.6/site-packages/django/contrib/gis/db/backends/postgis/operations.py", line 362, in parse_raster
    return from_pgraster(value)
  File "/home/.local/share/virtualenvs/AQM_API-RewsMqZ2/lib/python3.6/site-packages/django/contrib/gis/db/backends/postgis/pgraster.py", line 57, in from_pgraster
    pixeltype = POSTGIS_TO_GDAL[pixeltype]
IndexError: list index out of range

My hope was that using the above code would return a value "qs.rast = 34.1". Could anyone shed any light on this?

The query I am trying to replicate from PostGIS is:

#select ST_SetSRID(ST_MakePoint(-3.43, 52.946),4326) as loc,

ST_Value(Brt2.rast,1,  ST_SetSRID(ST_MakePoint(-3.43, 52.946),4326)) as pix_val from Brt2
WHERE ST_Intersects(Brt2.rast, ST_SetSRID(ST_MakePoint(-3.43, 52.946),4326))

I've come up with the following solution as:

qs = Brt2.objects.raw('select rid, ST_Value(brt2.rast,1,  ST_SetSRID(ST_MakePoint(-3.43, 52.946),4326)) as pix_val from brt2 WHERE ST_Intersects(brt2.rast, ST_SetSRID(ST_MakePoint(-3.43, 52.946),4326))')

for i in qs:
    print(i.pix_val)

This returns the desired value, but as far as I understand, raw sql queries are open to injection attacks, and I am unsure of how to protect against this. Ideally, I would like to use the GeoDjango ORM.


texts: Creating Locator using 2010 Census Block Group data

texts:
When I attempt to create an Address Locator in ArcMap using the shapefile (El Paso County, Texas) I download from the census website (t1_2010_48141_bg10.shp) I am not sure which alias names to match with the required field names.

Do I need to download a different shapefile?


texts: OpenLayers 3 WMS Styling Using SLD_BODY and POST request

texts:
I'm using OpenLayers 3 and GeoServer to create a spatial web application in which I want to enable users to choose a colour ramp (using Colorbrewer) and a number of breaks for styling an attribute in their analysis layer (i.e. to style a choropleth). I have been able to get this working by dynamically generating an SLD XML string and passing it in as the value to the SLD_BODY parameter in the TileWMS source, but only with 3 breaks/classes. When I select 4 or more classes nothing is returned and I believe the issue is that the SLD becomes too long for the OpenLayers TileWMS GET request. I have read a variety of posts that explain you can do a POST request in OpenLayers 2 to get around this issue, but I haven't found anything on how to do this using OpenLayers 3. Only a response in this post hints at a solution, but it's not enough for me to work out an answer. The following is how I'm currently defining the layer and adding the SLD. Does anyone know how to get around character limit with the OpenLayers 3 GET request by using POST instead?

window.layer = new ol.layer.Tile({
  source: new ol.source.TileWMS({
   url: 'http://' + window.location.hostname + ':8080/geoserver/wms',
   params: {LAYERS: geoserverLayer, STYLES: undefined, SLD_BODY: sld, CQL_FILER: geoserverFilter, TILED: true},
   serverType: 'geoserver'
  })
});


texts: What would cause Dojo in later versions of ArcGIS Javascript (v3.12) from working in Internet Explorer?

texts:
On my testing machine I experience no problems with using

http://developers.arcgis.com/javascript/samples/featuretable/

but if I run the same example on my client's machine I get "'require' is undefined".

I also experience the same issue with all of their 3.12 examples. Only 3.1 examples seem to work on my client's machine.

Update: It also works just fine in FF and Chrome on client machine


texts:
I'd like to use USGS data in other software to produce 3d terrains. Ideally I would use IMG format as its most common from what I see. But how can I convert it to a 32 bit exr format for usage as displacement or height map etc etc?

Found few converters and other programs but none of them work.


texts: How to load a data table as a 'vector layer' using PyQGIS?

texts:
This might be a simple problem, but how does one load a simple two column CSV with no associated geometric data?

I have given a couple go's via QgsVectorlayer[delimited text, spatialite] types - but none seem to work.

for example:

    theMinistryOf = ':\...\sillyWalks.csv'
    uri = QgsDataSourceURI()
    uri.setDatabase(theMinistryOf)
    uri.setDataSource('','','')
    notQuiteSillyEnough = QgsVectorLayer(uri.uri(), 'skip_step_slide', 'spatialite')

or:

    uri = theMinistryOf+'?delimiter=%s' % (",")
    sillyWalks = QgsVectorLayer(uri, 'sillyWalks', "delimitedtext")
    QgsMapLayerRegistry.instance().addMapLayer(sillyWalks)

This need is easily fulfilled manually by selecting load vectorlayer, but programmatically?

Can a CSV data-table be loaded via pyqgis without geometry?


texts:
I am building the plug-in for QGIS by following this tutorial.

When I am loading the plug-in, I am receiving the error The plug-in is broken, no module named resources.

Couldn't load plugin SavingAttributes due to an

  error when calling its classFactory() method.


texts: I have installed GDAL and checked on the command line, it is working and when I run the gdal2tiles.py on the command line it works but when I run the same file gdal2tiles.py with python IDLE, I get "No module named 'osgeo'" error. Can you assist.


texts: Fetch lat lon from eachLayer that uses coordinates

texts:
fairly new to leaflet and such and have been experimenting with a little project.
Essentially I have a csv file with latlng information and a function in the code that reads it and places markers accordingly.
My question is: How do I retrieve individual marker coordinates, convert them to latlng and apply a polyline in the order of marker appearance?
I've attempted using coordToLatLng(despite there being minimal documentation on its use) with negative results.
Below is my current code I'm working on:

<script>
var map = L.mapbox.map('map', 'faraday2.hb90naie', {
detectRetina: true,
zoomControl: true,
//take away the zoom controls on the top left
infoControl: false,
maxZoom: 15
});
map.setView([-34.0153, 93.1673], 4);

function eachLayer(layer) {
    var feature = layer.toGeoJSON();
    //if (feature.properties && feature.properties.data.title) {
        //layer.bindPopup("Name - " + feature.properties.title);
    //}
}

var moderateIcon = L.icon({
    iconUrl : './img/moderate.png',
    iconSize : [25, 25]
})

var highIcon = L.icon({
    iconUrl : './img/high.png',
    iconSize : [25, 25],

})

var extremeIcon = L.icon({
    iconUrl : './img/extreme.png',
    iconSize : [25, 25]
})

var points = omnivore.csv('./data/data.csv')

.on('ready', function() {
    map.fitBounds(points.getBounds());

    points.eachLayer(function(layer){
        if(layer.feature.properties.priority == 'high'){
            layer.setIcon(highIcon);
        }
        else if(layer.feature.properties.priority == 'extreme'){
            layer.setIcon(extremeIcon);
        }
        else if(layer.feature.properties.priority == 'moderate'){
            layer.setIcon(moderateIcon);
        }
        layer.bindPopup(layer.feature.geometry.coordinates);
        coords = [];
        coords.push(layer.feature.geometry.coordinates);
        //coordinates = coordToLatLng(coords);
        //window.alert(coords);

        /*var poly = L.polyline([], {color: 'red',
                weight : 5,
                smoothFactor : 1})
        for(var i = 0; i < coords.length; i++){
            poly.addLatLng(this.coords[i]);

        }*/

        var latlngs = Array();
        marker1 = L.marker([layer.feature.geometry.coordinates[1], layer.feature.geometry.coordinates[0]]);
        marker2 = L.marker([-30.4676, 115.7739]);

        latlngs.push(marker1.getLatLng());
        latlngs.push(marker2.getLatLng());

        var polyline = L.polyline(latlngs, {color: 'red'}).addTo(map);
    });

    var markers = L.markerClusterGroup({
        showCoverageOnHover: false,
        maxClusterRadius: 25
    });
    markers.addLayer(points);
    //window.alert(markers.length);
    map.addLayer(markers);
    points.eachLayer(eachLayer);

}).addTo(map);

    </script>

EDIT:
Sorry I didn't realise I didn't complete my question.
How would I also iterate through every coordinate I have and create a new variable accordingly for a marker so that polyline would work?

EDIT 2:
I can add a second variable for marker and give it hard coded location information and it will then render a line on the map, however it turns into a tree of lines rather than a line going from marker to the next to result in one long line.


texts: What do you call this combination of two district maps?

texts:
My real-world problem is that I want to combine state House and Senate legislative districts so that a user can pick a region (not a point) on a map that is included in both their House and Senate districts, to identify who their legislators are.

I have shapefiles for both House and Senate districts.  I'm looking for a library that will, given the two sets of shapes, emit a minimal set of shapes for which, by selecting one, a user may identify both their House and Senate district.  I can imagine how to do it myself by doing lots of intersection and difference operations, but also can imagine problems that others hopefully have solved.

I'm working in Javascript, have looked in turf.js, and not found this operation there.

My first problem is that I can't search for something to do this operation because I don't know what it's called.

Does it have a commonly-used name?


texts: Is it possible to calculate flow accumulation using google earth engine?

texts: I've done this with PCRaster and ArcGIS in the past but would like to know if there's a way to get it done in Google Earth Engine with the python API? Can you use other python libraries like pyDEM or pygeoprocessing with google earth engine objects?


texts: ogr2ogr importing geojson into PostgreSQL Error "Segmentation fault (core dumped)"

texts:
Having read everywhere as to the cause of these errors I'm at a loss. Trying to import geojson Point data into a new PostgreSQL database with the POSTGIS extension enabled.

PostgreSQL-9.4.17
Postgis 2.4.3
GDAL 1.10.1
Ubuntu 14.04
I've edited my pg_hba.conf to "trust"
I've created a table, granted all permissions, created password for superuser postgres. I'm able to import other SQL files
Edited permissions on the geojson file I want to import
GeoJSON is valid

The issue comes when I try to import a geojson file using ogr2ogr

ogr2ogr -f "PostgreSQL" PG:"dbname=dbname user=postgres" "file.geosjon" -nln tableName

And all i get in response is:

Segmentation fault (core dumped)

Using pgadmin if I inspect the table it imports all the correct column headers but no actual data.

I'm thinking its a permission issue or version issue after reading this .


texts:
I am working in Ubuntu Xenial with R v 3.2.3. I am trying to updated to the last version of R since I need to install the rgdal package to work with raster data in RStudio but I can't manage. I have tried already those options. After completing all the steps I runn R --version in terminal and it keeps showing R v.3.2.3

https://www.biostars.org/p/88651/

https://stackoverflow.com/questions/10476713/how-to-upgrade-r-in-ubuntu

Any idea about how to do it?

$ sudo apt-get update
    Reading package lists... Done
    E: Could not get lock /var/lib/apt/lists/lock - open (11: Resource temporarily unavailable)
    E: Unable to lock directory /var/lib/apt/lists/
    W: Target Packages (main/binary-amd64/Packages) is configured multiple times in /etc/apt/sources.list:50 and /etc/apt/sources.list:51
    W: Target Packages (main/binary-all/Packages) is configured multiple times in /etc/apt/sources.list:50 and /etc/apt/sources.list:51
    W: Target Packages (main/binary-amd64/Packages) is configured multiple times in /etc/apt/sources.list:50 and /etc/apt/sources.list:52
    W: Target Packages (main/binary-all/Packages) is configured multiple times in /etc/apt/sources.list:50 and /etc/apt/sources.list:52
    W: Target Packages (main/binary-amd64/Packages) is configured multiple times in /etc/apt/sources.list:50 and /etc/apt/sources.list:53
    W: Target Packages (main/binary-all/Packages) is configured multiple times in /etc/apt/sources.list:50 and /etc/apt/sources.list:53

I pasted the sources on the sources.list

The output of

 $ sudo apt-get update

Reading package lists... Done
E: Could not get lock /var/lib/apt/lists/lock - open (11: Resource temporarily unavailable)
E: Unable to lock directory /var/lib/apt/lists/


texts: Canada FSALDU 6-digit Postal Code Data Providers

texts:
I have been looking into finding point data for the Canadian 6-digit postal code (FSALDU). I have come up with 3 sources so far, 2 being far more expensive than the other. After speaking with the sales reps they all seem to be only offering me just the data. Does anyone know what is a causing the large difference in price? Are the more expensive ones more accurate or something?

Cheapest (by far)
Statistics Canada

Most Expensive
DMTI Spatial , www.maponics.com/products/gis-map-data/canadian-postal-codes-database/overview


texts: Narrow viewport => not showing map in browser?

texts:
I have a webapplication using the arcgis javascript API and it works great. Unfortunately I ran into a small peculiarity:

If the viewport for the map is smaller than 930 pixels it refuses to load my map. I get the zoom in/zoom out button and the scalebar, but the map area stays empty. If I resize my browser slightly and reload the page the map shows.

The annoying in all this is that, as soon as the initial load is complete, I can perfectly resize the map area to 100 pixels or less and it still works. It's just the loading step where it breaks.

I've tested this in the latest firefox, but chrome exhibits the same behavior.

Any ideas/pointers where to look?


texts: What is the most stable implementation of GeoPandas for sjoin and Overlay?

texts:
I am working with GeoPandas in VS Code on a Windows 10 machine in Python 3.7-32 bit and GeoPandas 0.8.1.
For the script I am working on I also need to get data from a SQL server, and need functions from shapely that are not in the main GeoPandas so I also have

pyodbc(4.0.26)
pandas (0.24.2)
shapely (1.7.0)

I need to read one data set from a .gdb file, so that goes straight into geopandas and one from a SQLServer db. Which I read in the WKT and then use shapely to convert it to a geopandas dataframe.
gdb = geopandas.read_file(gdb_file)
gdb = gdb.to_crs("EPSG:2913")
gdb ['buffer'] = gdb.geometry.buffer(25,cap_style=2)
gdb = gdb.set_geometry('buffer')

cnxn  = sqlserver_con(constr)
sql = '''SELECT OBJECTID, SHAPE.STAsText() as WKT
        FROM dbo.Source'''
df = pd.read_sql(sql,cnxn,coerce_float=False)
cnxn.close()
df['LINE'] = oh.WKT.astype(str)
df['LINE'] = oh.LINE.apply(shapely.wkt.loads)
gdf = geopandas.GeoDataFrame(oh[['OBJECTID','LINE']],geometry='SEGMENT', crs="EPSG:2913")

I have more data than I can hold in memory so I am looping through chunks from both the gdb and the SQL connection (WHERE clause is not shown).  Every loop there is a spatial join. I have noticed there are some times it works flawlessly and others it fails. I have already tested:

reset_index on both sides
making sure there will be at least one geometry result for geopands

I started testing the spatial index gdf.sindex and ran into a weird problem:
gdb = geopandas.read_file(gdb_file)
gdb.sindex
gbd.sindex
gdb.sindex

This crashed Python every time. It works the first two times, but not the third. This is true even if the .gdb file has one record.


texts:
I would like to install Geoserver silently on Windows machines. I tried running the windows installer from command line, passing to it /S argument:

.\geoserver-2.13.4.exe /S

It works, but I need to pass also some installer's parameters such as the port number. Also the path of the Java VM need to be passed, as the silent installer didn't detect it automatically.

How can I do that?


texts: I have AIS Data (MMSI, timestamp, latitude, longitude) in Excel files that I want to visualise using density plots/heatmaps, particularly being able to observe which regions have higher number of ships going through them. Since I am completely new in this area, any insight on the initial steps and which tool would be suitable for this would be of great help.


texts:
I am making a map for my thesis. While using ArcGIS, I must use Import71 to convert some .e00 files (in this particular case, political boundaries and Hydrography networks). When I select my output folder, I receive the message:

  E:\ File Path This file cannot be accessed. Make sure you have security privileges on the drive.

I use Windows 7. Is this simply a Windows/administrator problem or am I doing something wrong? If someone knows, I would appreciate any suggestions.


texts: Trying to add ECW raster format to geoserver 2.9.3 stores list

texts:
I am trying to add the ECW raster format into geoserver (2.9.3) stores list as a raster data source.

I followed the "http://docs.geoserver.org/latest/en/user/data/raster/gdal.html" tutorial and i managed to add gdal extension and now i can see more raster data formats but not ecw.

I also tried to follow the installation of ecw sdk from this link:
"https://trac.osgeo.org/osgeo4w/wiki/pkg-gdal-ecw"

But i have still an error occuring when starting geoserver :

ERROR 1: Can't load requested DLL: C:\Pgrogram Files\GDAL\gdalplugins\gdal_ECW_JP2ECW.dll
126: The specified module could not be found.

spec : windows 7 64 bits, java 8 64 bits, MSVC2010, geoserver 2.9.3, imageio-ext 1.1.16, gdal-1.9.2


texts:
I am new to GIS software. I've played with Google earth a bit and am familiar with downloading Canadian topographic maps. Anyways, I've been playing around with qgis. I have imported a few georefferenced tiff files and a gpx file.

I am having difficulty displaying the raster maps correctly though. The maps looks squished in the vertical dimension. If I zoom in each pixel is not square, it is a short/fat rectangle. I feel like it has to do with the coordinate system?

Does anyone have any ideas what i might do/read up on in order to figure out how to make the pixels square and the maps not stretched?


texts: Pulling a CSV file from website into QGIS with ability to automatically refresh?

texts:
I am attempting to pull a CSV file from a website that is updated daily. I would like to have a process that can automatically perform the process and save daily as a new layer. The CSV file is related to tracking forest fires and has multiple attributes. The CSV download url is http://www.wildfire.alberta.ca/reports/activedd.csv. I attempted to use the a vrt driver to load the file into QGIS as a layer but I was unable to pull the attributes as well as the lat/long locations.

I am fairly new to QGIS.


texts: Increase each polygon area to the same area

texts:
I'm using ArcGIS Desktop 10.2.

I have many polygons (buildings) with various areas and I want the same area for all my polygons: 100m². The final goal is to aggregate them and build the urban agglomeration.

I was thinking of creating a buffer with a value for each polygon and using the field value to generate the buffer, but the challenge is to find the value.


texts: How to preprocess and remove noisy data from Landsat-8 image Bands for ndvi time-series analysis

texts:
I'm working on satellite image processing, I'm calculating NDVI for time series analysis and vegetation prediction. I'm having a tough time removing the noisy data like,

Cloud coverage
Haze
Cloud shadow

on Landsat-8 satellite image bands.

Is there any way to do this in Python?

I came across fmask package in python to mask cloud coverage but wasn't clear how to use it.

Then used pymakser but was only able to create a separate cloud mask file using BQA band of Landsat-8 but couldn't apply the mask on the image bands.

Finally tried using rio-cloud mask package but wasn't able to create the pre-requisite top of atmosphere reflectance for landsat-8 bands.

Also it'll be great if there's a method to just use the BQA band of Landsat-8 to remove noisy data from the TIF files using python.


texts: I'm working on forest biomass estimation with Alos-Palsar and Landsat-8 data for 2015. For this I need to coregister the scenes from Alos-Palsar and Landsat-8. I would like to know which open source software can I use to perform this task.


texts: Leaflet Custom TileLayer for ArcGIS Cache Tiles in Folder

texts:
This is ArcGIS Server tile cache folder structure on a server.

I need to get these tiles on my Leaflet map. I need to create a new TileLayer. I will extend existing TileLayer layer object's getTileUrl() method.

L.TileLayer.MyCustomLayer = L.TileLayer.extend({
    getTileUrl: function(coords){
        ?????
    }
});

I need to create new tile url address by coords parameter of getTileUrl() method.


texts: ArcGIS symmetrical difference in arcgis desktop standard

texts:
I have an ArcGis Desktop Standard license and I would like to perform an operation to find out the difference between 2 features classes geometries (Something like symmetrical difference in Arcgis advanced version). Unfortunately symmetrical difference is not available for Arcgis desktop standard. The operation should be performed using Arcpy. Has anyone a tip for me, how to perform such an operation on Arcgis Standard?

Thanks


texts: How can one compute & render velocity from kml track in qgis?

texts: I have a bunch of KML files recorded from activity trackers.  These record a stream of coordinate.  These lack individual timestamps, but one can assume that the samples were taken at regular intervals.  Is there a way of having qgis compute vector deltas between adjacent coordinate-pairs, so as to estimate a velocity?  Then I'd like to have the velocity value influence rendering of the shape.


texts: How to extract coordinates & elevations of boundary cells of DEM?

texts:
I am using a 50mx50m grid irregular shaped DEM file.

I have to work on the boundary conditions.

How do i extract the cell locations ( in x,y coordinate ) along the boundary edge and consequently find the elevation of the respective boundary cells.

I am using ArcGIS 9.3 for my work.


texts: QGIS SAGA's Principal Component Analysis (PCA) makes output rasters (components) bigger

texts:
I have multispectral TIFF bands/features (dtype = float32) that I passed onto QGIS SAGA Principal Component Analysis tool to perform PCA. The bands are each a separate file but have exactly the same extent and CRS (WGS UTM Zone 33N - EPSG:32633). I used the method [0] correlation matrix and chose 3 components to output.
The issue I am facing is that the resulting components have enlarged their extent. In the picture below, there is a pseudocolor NDVI feature-image (one of the inputs to the PCA). The greyscale image below is one of the resulting components. As you can see, the component's area has grown, to me, for no reason.

Has anyone met with such issue? Or can reproduce? Attached you find the sample data: https://drive.google.com/file/d/1vSkaPsQUqE_M-uIKnYSZzTkiqAC12ZKZ/view?usp=sharing
One reason could be that SAGA creates exclusively SDAT files. I am not familiar with SDAT, could they cause that because of wrongly handling the CRS...? I nevertheless suspect a potential bug.
Software used:
QGIS v3.18.2
SAGA v7.9.0 using Next Gen Provider that should facilitate compatibility with QGIS
Info on data:
PCA inputs: Sentinel-2 imagery and extracted features, many times modified and rewritten using Python rasterio, always as float32, consistent CRS EPSG:32633. The features are temporal characteristics with complicated logic behind, probably not key for the issue.
Once difference for SAGA PCA result was that the SDAT components had no-data values as -99999and no longer nan (such as Python's numpy.nan).


texts:
I am having the same issue on QGIS 3.10 or the long term release, 3.4

If I go to Settings > Custom Projection and add a name and parameters, the "OK" button does nothing.

I have tried clicking the plus button before and after adding the name and parameters, it makes no difference.

I have read similar questions:

QGIS custom CRS won't save
Why can't I save my custom CRS?
Can't define Custom CRS in QGIS

but none of their answers work for me.

I have tried uninstalling and reinstalling different versions of QGIS, and have tried granting write permissions to folders within the QGIS directory.


texts: Programatically populate field_4 after choosing values from other fields?

texts:
I'm working with QGis 2.8, editing the attribute table of a SHP file, named  "myTable".
I need to populate field_4 with the conditions listed below:

For 1 to n (it runs across records of myTable)
  select record n
    set var1 = field_1.value (it puts the value of field_1 in var1)
      if field_2 >1 then field_4.value = 'OK'
      else
      SELECT MAX(field_2) FROM myTable WHERE field_1.value = var1 (it selects the record with the higher value of those who have field_1=var1)
      set var3 = field_3.value (it copies field_3 value into var3)
      end
    select record n (it come's back to current record)
    set field_4.value = var3 (it populates field_4 with var3 value)
next n

myTable will be like this after running:

geometry, field_1, field_2, field_3, field_4

1       1       10          a       ok
2       1       1           b       c
3       1       80          c       ok
4       2       60          d       ok
5       2       20          e       ok
6       2       1           f       d
7       2       5           g       ok
8       1       1           h       c
9       1       8           i       ok
10      2       5           l       ok
11      2       5           m       ok
12      2       4           n       ok
13      3       39          o       ok
14      3       1           p       r
15      3       10          q       ok
16      3       50          r       ok

I'd like write this script in Python, but my skill is very low !!
(note that it calls some SQL functions)
I don't know how to use "cursor" tecnique into Python.
Can someone help me writing this procedure?


texts: How can I update the label of a Polygon Vector after it has been drawn?

texts:
I am quite new to OpenLayers. Right now, I have a polygon vector with some styling applied and a label.

var style = $.extend(true, {}, OpenLayers.Feature.Vector.style['default']);
style.pointRadius = 15;
style.label = "My Polygon";
style.fillColor = #f00;
style.strokeColor = #000;

var styleMap = new OpenLayers.StyleMap({"default" : style});

var polygonLayer = new OpenLayers.Layer.Vector("Polygon Layer", {styleMap: styleMap});

At some point after doing some processing, I want to display the result as a label. How can I update the label? I figure it would be something like this

polygonLayer.options.styleMap.styles.label = "Updated label";


texts: Choose an attribute dropdown for ArcGIS webbuilder app

texts: I'd like to build a simple web mapping application using ArcGIS Online and Web App Builder.  I would like a user to be able to choose an attribute via a dropdown menu and see the symbolized results on a polygon chloropleth style map.  Is the only way to do this to have multiple layers, each representing one attribute?  Or is there a widget I am missing that allows the interactive selection of attributes to visualize?


texts: Change TIFF image projection without losing layers in QGIS 3.2?

texts: I have a .tif file (see color image below) that has 10 bandsets corresponding to a different color in the image. When I import the image to QGIS I can see all 10 layers. I need to reproject the image so it has units of degrees, not meters. I use  Raster -> Projections -> Warp (reproject) and then choose target CRS as EPSG:4326-WGS 84. However, when I run the program I get my reprojected image but I lost all of my layers. All I get is a binary image (see bw image below). How do I reproject my .tif  image so it maintains its color scheme and layers?


texts:
I have a shapefile of the country. I need to select only percentage of that polygon for demonstration purposes. As an example, how to select 60% of polygon?

I am using ArcGIS 9.3.


texts: I'm a beginner. I selected 9 polygons. I want to change their color to green. They should stay green after the selection is removed.


texts:
I am writing a processing algorithm and I would like to return a feature sink as an output, but without treating it as a parameter. The point is I want to just add it to the table of content as a temporary layer, without showing it in the parameters window.
However, I can't find any other way to create the feature sink other than by declaring parameter QgsProcessingParameterFeatureSink and referencing it by self.parameterAsSink.

Is there a way to return a feature sink as an output ? Or maybe it should be done with another approach?


texts: Delete polygonal features in one layer outside polygon feature of another layer in PostGIS

texts:
For deleting polygonal features in one layer (~1000000 features) outside polygon feature of another layer (one feature) in PostGIS I use query:

DELETE FROM layer WHERE layer.id IN (SELECT a.id FROM layer a, boundary b WHERE ST_Disjoint(a.geom, b.geom));

Now the process takes more than 6 hours. How can I speed up the process of deleteng big count features one layer outside feature in another layer?


texts:
I know this will sound like a foolish question once I find out what I'm doing wrong but here goes:

I had downloaded an outdated version of my data collection map onto Collector, but I need to switch it out for the more current version. After deleting the content of the previous map I tried to download the newer version, but to no avail! Its not that the download didn't work, its just that there was no download option. I tried looking this problem up online, but I couldn't find any similar cases (which lead me to believe I'm probably doing something stupid). Would this be a problem with my ArcGIS account or my iPhone that I am downloading the map onto?


texts: Select a polygon layer where a line touches its boundary more than 1 time

texts:
My problem is that I don't know how to select the layers (polygons) where the lines touches the boundarys twice or more. I have tried to illustrate this in the picture.

(see image)
I am using arcmap desktop advanced

Any tip or workaround?


texts: ArcGIS Pro dynamic text from attribute table in layout. Get first value

texts:
I have feature class with many objects in extent and I need to add dynamic text to layout with only one value from Text field. All object in current extent has same field value.
Has anyone figure out how to do it?
Tag in text element
"<"dyn type="table" property="value" mapFrame="Фрейм карты" mapMemberUri="CIMPATH=__________/______2.xml" isDynamic="true" field="LH_mery.Lesnich" delimiter=" "/> –
What i have

What i need


texts: Getting List elements from .reduceRegion() results in Google Earth Engine?

texts:
I have a FeatureCollection that I am using .reduceRegions() over the Hansen Forest Loss dataset using ee.Reducer.fixedHistogram() to get the count of pixels per year/per feature.

The result is a feature which I can see the results in the 'histogram' property but can't select those results to store in another variable to do further analysis.

How can I select the results from the reducer (ex: how can I select 3.011764705882353 in the first element)?

I keep getting the following error:

ComputedObject (Error)
List.get, argument 'list': Invalid type. Expected: List<Object>. Actual:
Float<dimensions=2>.

Here is the code I have:

var forest_years = forest.select('lossyear')

var test = forest_years.reduceRegions({
  collection:GB_patches,
  reducer: ee.Reducer.fixedHistogram(0,18,18),
  scale:30,
})

var first = ee.Feature(test.first())
var patch_years = ee.List(first.get('histogram'))
print (first)
print (patch_years)
print (patch_years.get(1))

Link to GEE code: https://code.earthengine.google.com/9ab9b35f4cc49b45f29386c70a0cb2cc


texts: Export DBF from shapefile with ModelBuilder (to make a XY event layer)

texts:
I'm using ModelBuilder and I need to export a shapefile to a DBF. Usually I go to "table option" and I click to "Export", how can I do it with ModelBuilder?

I need to do this process because I have to use the XY Event Layer which requires a DBF to work correctly (if you put a shapefile, the function duplicate the geometric attributes).


texts:
I'm trying to convert this file to EPSG:3857 with this command:

gdalwarp -t_srs EPSG:3857 S2A_rhow_l2r.nc test-proc.nc

Warning 1: No UNIDATA NC_GLOBAL:Conventions attribute
Warning 1: Several drivers matching nc extension. Using NETCDF
ERROR 1: Input file S2A_rhow_l2r.nc has no raster bands.

I have no idea how to solve this problem.


texts: Mean coordinates of polygon in PostgreSQL

texts:
I need a way to get the mean (or median) coordinate of a polygon. With QGIS this is quite simple. But I need to do it with PostgreSQL.

There is the expression "st_GeometricMedian". But this only works with Points.

I thought that I am going to transform the polygon into the points with and get the median:

select
sourceeid,
ST_GeometricMedian((st_dumppoints(the_geom)).geom)
from
polygons_view
where
sourceid=67475

But the result is that I get a line for every Point and then also a median for each line. But I need for one polygon one line.

How can I do this?


texts: Python GDAL open dataset in PostGIS with quotations

texts:
Error trying to use gdal to open a dataset in PostGIS with the following code:

ds = gdal.Open("PG: dbname=nyc host=localhost user=postgres password=mypassword port=5432 mode=2 schema=public column=rast table=mytable where='collected_date > 2014-04-04 07:00:00'", gdal.GA_ReadOnly)

The log file in PostGIS shows:

  2016-04-18 21:09:36 MYT ERROR:  syntax error at or near "07" at
  character 336

  2016-04-18 21:09:36 MYT STATEMENT:  select srid, nbband, ST_XMin(geom)
  as xmin, ST_XMax(geom) as xmax, ST_YMin(geom) as ymin, ST_YMax(geom)
  as ymax, scale_x, scale_y from (select ST_SRID(rast) srid,
  ST_Extent(rast::geometry) geom, max(ST_NumBands(rast)) nbband,
  avg(ST_ScaleX(rast)) scale_x, avg(ST_ScaleY(rast)) scale_y from
  public.cea where collected_date > 2014-04-04 07:00:00 group by
  ST_SRID(rast)) foo

Apparently there is a syntax error in the date condition (in bold), where it needs to be enclosed with quotations. But I can't seem to be able to insert the quotations at all as they are always removed when the string is parsed.

I guess we can replace the date with mktime function instead, but I'll looking for a solution that can be extended to string datatype as well.


texts: Ogr2ogr unable to export S57's Symbol to Shapefile

texts:
I am using ogr2ogr tool to convert the S57 files to shapefile. However in the resulted shapefile, I am unable to show the original icons from the S57. For example, the SOUNDG layer be represented by number (depth), but my shapefile just shows a circles when imported to Arcgis map.

I would like to know is there any way to keep the depth information of the S57 while convert to shapefile or even kml file.


texts: Vector lines as blockage for viewshed analysis in QGIS

texts:
I have a SRTM elevation model with a point from which I would like to calculate the aproximate visibility. Using the line-of-sight tool in QGIs, I can achieve that without any difficulty. However, I'd like the tool to take in consideration the vector lines I have drawn around the point (the lines should serve as a blockage and visibility shouldn't pass through them).

Do I have to modify the raster DEM where this vector lines are (and set them to a height of 99999) or there is a way to make QGIS consider the vector lines while doing the viewshed analyses?


texts: Merging two pyramid tile sets (directory)?

texts:
I was searching on the Internet about the merging of GeoTIFF tiles (pyramid tiles) generated by gdal2tiles.

Below are the steps I performed

Downloaded tif-files.
Converted tif to vrt (below command)

gdalbuildvrt temp_intermediate/merged.vrt *.tif

gdal_translate -of vrt -expand rgba temp_intermediate/merged.vrt temp_intermediate/merged-4326.vrt

Translated vrt to pyramid tiles (below command)

gdal2tiles.py -z '14-21'  temp_intermediate/merged-4326.vrt  output/

Now if I get another set of tif-files and performed the same operation I'll be having two pyramid tiles directory (lets say outout1 and outout2 directories) and I want to merge them. How could this be done?


texts: iface.addVectorLayer on a gml file causes "Select vector layers to add" dialog (pyqgis)

texts:
I am trying to load a directory full of zipped GML files (*.gz) into QGIS using pyqgis' iface.addVectorLayer(), but as the GML file contains multiple separate vector layers (Area, Text, Symbol, Point and Line), the "Select vector layers to add" dialog box comes up for each GML in the folder. This means that rather than saving me time by loading all the files in at once, I would still have to click on the dialog box for every layer in the folder.

The ideal situation is that all of the vector layers within each GML are added into QGIS without any manual intervention. How does one select all of the vector layers in the dialog and add them to the gui without manual intervention?

I'm running this from the python console within QGIS.

import os
for r, d, fi in os.walk(path):
    for f in fi:
        iface.addVectorLayer(os.path.join(r, f), f, "ogr")


texts:
Is there a way to manually place labels in QGIS?

I have a map which has almost-duplicate labels in some places (e.g. E 1st St and 1st St) and I'd like to remove some of them.

Also, I'd like to move some of the labels around and maybe label some features that aren't currently.

I'm using new labeling, with the labels curved along the features.


texts: How to view URL of GetFeatureInfo request in QGIS?

texts: I have a WMS layer in QGIS. I am using the Identify Features tool to perform a GetFeatureInfo request. It retrieves the expected attribute data in the Identify Results panel. However I would like to see the URL of the request too. How can I achieve that?


texts: Broken colors or black background with OpenLayers Plugin in QGIS

texts:
Approximately last few months I am experiencing problems with OpenLayers plugin in QGIS - especially with OpenStreetMap data.

The layer loads sometimes terribly slowly although the connection is quick enough. Sometimes during map panning I get damaged colors or black background.

I had similar issues both at work (latest 2.18 LTR, 64bit Windows 7) and at home (latest 2.18 LTR, 64bit Kubuntu linux).

Is there any issue in OpenLayers plugin or it is just caused by increased traffic on OSM servers which are temporarily overloaded?


texts:
There are a lot of POIS in the database with fileds id name address lon lat.

When we tried to query some points for a given bounding box, it is simple to build the accordingly sql and get the results, however we found that the distribution of the result are not desired(see the left image below):

While we want the results of every page distributed evenly in the view port like the image in the right.

Is this possible?


texts: Which query to use to find closest items to a lat, long pair using PostGIS?

texts:
I have clients (mobile and web) who can request nearest groups to them by providing a lat, lng pair.

Here's a basic snippet of my query, and I am not sure which one I should use for finding nearest groups.  I believe the 2nd one would be easier for me to use for keyset pagination, but anyway here they are:

SELECT
  *
FROM
  groups
ORDER BY
  groups.geography <-> st_setsrid(st_makepoint(-122.057608, 37.37761), 4326)

OR

SELECT
  *
FROM
  groups
WHERE st_dwithin(groups.geography, st_setsrid(st_makepoint(-122.057608, 37.37761), 4326), 5000)

Which approach would be better to use to find the nearest results in order?

The second one is not ordered though, but the documentation says to use st_dwithin for finding results within a range of a lat, lng.


texts: Easiest way to programmatically clip raster to rectangular extent

texts:
I am looking for the easiest way to programmatically clip a raster to a rectangular extent for testing purposes. I have GDAL and Python installed on a virtual machine and would like to run the code from Jupyter (python/command line) or command line directly. I do not have a polygon with the extent yet. I am not looking for answers that involve QGIS/ ArcGIS or some other GUI.

I found this resource but I was hoping for something more succinct like:

some-command inputraster.tif outputraster.tif -extent [lower left corner, upper right corner]

or

Create polygon
Clip using polygon

I am pretty sure there must be some gdal command out there.


texts:
I'm trying to download all features from this WFS: capabilities, describefeaturetype

It's working fine with the most simple requests like: getfeatures

But now, I need >all< features in the dataset, but the service won't allow me to do so since the maxFeatures is set to 100 (serverside). Now, I'd want to create some sort of quadtree to iterate through the data. But therefore I need a boundingbox filter. I found some examples like;ogcfilter

and one using a cql_filter:cqlfilter

None of them is working. Is there someone who can guide me in the right direction? I'd like to know how to use the boundingbox filters, but still, ultimately I'm looking in downloading >all< features within a certain boundingbox.

Edit:
I also tested the WFS server in QGIS. Without enabling the viewport as a boundingbox, it works fine. As soon as I say I want only the data within the viewport, it gives an error on the bbox... Is this maybe a server-side problem? Or is it just the QGIS method doesnt work on this specific Service?


texts:
I have an access database which holds several sites. This is linked to a MapInfo TAB file which holds geometry for each site.

I would like for our non-GIS users to be able to view the sites in an interactive way (i.e. a slippy map with layers, not just a static image - I've already got static images in there) rather than having to boot up MapInfo or, worse, ask someone else to...

I have embedded a MapInfo mapper inside a form but distributing this has proved to be more expensive than I originally expected (though to be fair my expectations may have been unreasonable!).

As far as I know it is not possible to embed HTML in a form - otherwise I would use OpenLayers or similar.

I could include a link to a web-based map but that isn't ideal...

Any suggestions? What software is available for this sort of thing?


texts: Getting intersection of circles using shapely?

texts:
I have a List of Points which I transform to circles. Code from here

circles = [Point(random.random(), random.random()).buffer(random.random() * 0.1) for i in range(3)]

I want to know which circles within this List intersect each other. Is there a way to achieve that with a native shapely function, or do I have to iterate over the all possible combinations?


texts: Error 400 when cloning GEE script locally

texts: I need to modify google earth engine (GEE) scripts on my mac. I was hoping to run git clone https://earthengine.googlesource.com/users/... to get a local copy. Unfortunately, I get this response after a query for username and password unable to access 'https://earthengine.googlesource.com/users/...': The requested URL returned error: 400


texts: Got "No such layer" response through RESTful API

texts:
We publish a feature type and layer called "ice" through GeoServer RESTful API and then we close the computer without exit the GeoServer instance normally. After that, we restart the computer and GeoServer service, found this layer 'ice' could not be found through restful api 'geoserver/rest/layers/ice'. But the 'featuretype.xml' and 'layer.xml' could be found in the 'ice' dir under data_dir.  After that，I found the problem is the DataStoreInfoImpl value in featureType.xml is not consistent with that in datastore.xml. After I correct featureType.xml and restart GeoServer , the layer could be loaded as expected.
But I wonder what is the ‘DataStoreImpl' meaning, and if the datastore for example PostGIS, or the GeoServer was crashed unexpectedly and restart later, those layers published before would not work since 'DataStoreImpl’ has been changed?
The datastore is PostGIS. And the log in GeoServer start process shows:
Failed to load actual store for FeatureTypeInfoImpl[ice]
25 八月 15:26:26 WARN [org.geoserver] - Failed to load resource
java.lang.IllegalArgumentException: Resource must be part of a store


texts: Moving rotated polygons horizontally to avoid overlap

texts:
To rotate rectangular polygons, I use the following (courtesy from this answer):
rotation = 50
vlayer = iface.activeLayer()
provider = vlayer.dataProvider()

couples_id_geom = []
for feature in vlayer.getFeatures():
    geom = feature.geometry()
    centroid = feature.geometry().centroid().asPoint()
    geom.rotate(rotation, centroid)
    # accumulate args to avoid rotation feature by feature
    couples_id_geom.append([feature.id(), geom])

# Change the layer features rotation in one go
provider.changeGeometryValues({
  couple_id_geom[0]: couple_id_geom[1] for couple_id_geom in couples_id_geom
})

# Refresh to see the changes
vlayer.triggerRepaint()

However, when rotating adjacent polygons, overlaps occur. For example at 50 degrees:

How can I set it like the last pair of rectangles shown in the image so that when the polygons are rotated, they are moved at a minimum horizontal distance to avoid overlap and print the distance required?
EDIT
I asked a question on Math SE about calculating the minimum horizontal distance but I am unsure how to program this into Python.


texts:
Is anyone aware of a script or tool that would allow me to view what users have a paticular mxd open at one time?

We are using ArcGIS Desktop 9.3.1.


texts: Geotools creating polyhedral surfaces from polygons

texts:
I have polygon geometries that represent building footprints in a shapefile and I want to extrude them according to an attribute "height" and turn them 3B objects as polyhedral surfaces via extrusion. For this purpose I am writing an extrusion function. I am using geotools to access and extract geometries and attributes from shapefile. I am creating top and bottom faces of buildings as polygon geometries but I could not generate vertical faces. I could not write a loop that traverse positions and create vertical faces as polygon geometries. I am stuck with it.

Here is the rest of my code after I extract polygons ;

List<Polygon> polygons = bb.getBboxFeature(shapeFilePath);

List<PolyhedralSurface> mesh = new ArrayList<PolyhedralSurface>();

for(Polygon poly : polygons) {

    Coordinate[] coords = poly.getCoordinates();
    List<Polygon> faces = new ArrayList<Polygon>();

    for (int i =0; i<coords.length; i++) {

        Double h = 15.00;
        GeometryBuilder builder = new GeometryBuilder(DefaultGeographicCRS.WGS84_3D);
        List<Position> dpsBottom = new ArrayList<Position>();
        List<Position> dpsTop = new ArrayList<Position>();

        //bottom positions created using coordinates

        dpsBottom.add(builder.createDirectPosition(new double[] {coords[i].x, coords[i].y, 0}));

        //top positions created using coordinates

        dpsTop.add(builder.createDirectPosition(new double[] {coords[i].x, coords[i].y, h}));

My final purpose is create vertical faces and add them to a polygon list. Then, finaly giving this list to createPolyhedralSurface() function as input and create building as polyhedral geometry.

Is there a more efficient solution for my purpose?


texts:
I have a map in QGIS displaying many languages being spoken in the displayed area. Each one of them has a designated marker with a HTML frame that shows up on hover.

However, in this frame I would like to display data like the ISO-code of that language. This data is stored in a text file that I added to the project via layers.

It contains the X Y coordinates, the language name (displayed on the map) and the ISO-code, each tab separated.

The HTML frame code is this (http://pastebin.com/xZYCk6S0 if it's not working properly):

<!DOCTYPE html>

Search this ISO code on

WALS

This language is found here:

I haven't found out yet how to reference the data on the map so that the frame would say:
"language_name" with ISO-Code "ISO_code" is found here:
So that each frame adaptively has the right information in it.


texts: Is there a bug in Georeferencing in QGIS 2.0?

texts:
I recently updated to QGIS 2.0. But now, when I need to georeference some pictures, I have serious problems.

When loading the files to be georeferenced, I can see only a part of it in the window, so that working with it is not possible. Zooming out or in, zooming to the layer doesn't help anything:

First I thought that the size 2362 x 2637 is not suitable to QGIS, so I reduced it to 567 x 633 at 72 dpi, but the problem occurs.

What can I do? Thanks a lot for your hint!!!

Harald


texts: Dividing or splitting polygon to smaller polygons with PyQGIS?


texts: How to get the lat/long from highest point in a polygon, on earth engine?

texts:
I have a polygon (fusion table) and I need to discover which is the highest point from point. Here's what I have done until now:

var elevation = ee.Image('CGIAR/SRTM90_V4');
var rainierWaypoints =
    ee.FeatureCollection('ft:1oCYHTwkJhzDdJi3sabCm-BFyTwhu7N2HuQPpceBv');

Map.addLayer(elevation.clip(rainierWaypoints), {min:0, max:300, palette:['000000',"ffffff"]},"elevation2");
var chart = ui.Chart.image.byRegion({
  image: elevation,
  regions: rainierWaypoints,
  scale: 200,
  xProperty: 'name'
});
chart.setOptions({
  title: 'Propriedade',
  vAxis: {
    title: 'Elevation (meters)'
  },
  legend: 'none',
  lineWidth: 1,
  pointSize: 4
});
var options = {
  title: 'Altitude',
  fontSize: 12,
  hAxis: {title: 'Frequência'},
  vAxis: {title: 'Altitude'},
  series: {
    0: {color: 'blue'}
  }
};

var histogram = ui.Chart.image.histogram(elevation, rainierWaypoints, 30)
    .setSeriesNames(['blue'])
    .setOptions(options);
print(histogram);

print(chart);

Map.centerObject(rainierWaypoints);
Map.addLayer(rainierWaypoints, {color: 'FF0000'});

Here's the fusion table that I'm using as example.
https://fusiontables.google.com/data?docid=1oCYHTwkJhzDdJi3sabCm-BFyTwhu7N2HuQPpceBv#map:id=3


texts: Quick way to do repetitive task? Point feature class (paste one field value from one feature class to another)

texts:
I am identifying certain points in one feature class within ArcMap and want to quickly make a separate point feature class in the same location and paste a value from the first feature class into the new point feature class.

Is there a quicker way to do this than my current method, which is to copy the feature and paste it (into the new feature class then paste the attribute from one to the other)?


texts:
Following this tutorial (Rails Port, OpenStreetMap website) https://github.com/openstreetmap/openstreetmap-website/blob/master/INSTALL.md , I successfully installed my local server and it is working. When I go to my server adress, the website just looks like OpenStreetMap's website.

Now I want to use in JOSM my local server instead of api.openstreetmap.org/api alt text But I have this error :
Failed to retrieve a list of changesets from the OSM API server at 'http: //localhost'. The server responded with the return code 404 instead of 200. Please check the spelling and validate again


texts:
I've set up a WFS layer in OpenLayers:

 var v1497_shp = new OpenLayers.Layer.Vector("v1497_shp", {
     strategies: [new OpenLayers.Strategy.BBOX(), saveStrategy],
     protocol: new OpenLayers.Protocol.WFS({
         version: "1.1.0",
         // loading data through localhost url path
         url: "http://localhost/geoserver/wfs",
         featureNS :  "http://www.opengeospatial.net/cite",
         maxExtent: mapextent,
         // layer name
         featureType: "v1497_db_wfs",
         // geometry column name
         geometryName: "the_geom",
         schema: "http://localhost/geoserver/wfs/DescribeFeatureType?version=1.1.0&typename=cite:v1497_db_wfs"
     })
});

Now I want to highlight the layer when I pass over with the mouse.

Can I do this?


texts: Transfering from pandas dataframe to Geopandas?

texts:
I am not able to transfer a pandas dataframe column with geolocation coordinates into a Geopandas Point format.

location = pd.read_csv('City_addresses')
print location['Geolocation'].head()

Outputs the following:

0    (50.673675, -120.298973)
1    (50.678354, -120.329258)
2    (50.672496, -120.333317)
3    (50.673359, -120.332912)
4    (50.673411, -120.32978)
Name: Geolocation, dtype: object

Then I attempt to transfer it to Geopandas Point coordinates.

from shapely.geometry import Point
from geopandas import GeoDataFrame

location['Geolocation'] = location['Geolocation'].apply(Point)
gpd = gpd.GeoDataFrame(location, geometry='Geolocation', crs={'init': 'epsg:4326'})

print gpd

To which I get:

  TypeError: a float is required

But it won't allow me to change the datatype to a float in pandas either. What should I do?


texts: What dataset contains a complete location hierarchy from continent down to city (or perhaps postal code)?

texts:
I'm looking for a dataset that contains a complete hierarchical geographic representation of the world. Basically, I want names, boundaries, relationships ("is a member of"), ISO codes, and more for:

continent (optional but useful)
country
state/province/territory
district/county
city

Postal code would be nice, too, but that can complicate things since postal codes can sometimes span multiple regions (e.g., ZIP code 02861 covers area in both Pawtucket, Rhode Island and Attleboro, Massachusetts).

Some sort of globally unique ID for each location would be really useful. ISO 3166 specifies country codes (3166-1) and sub-country (state/province/territory generally) codes (3166-2), but it doesn't seem to contain any finer administrative units. Is there an ISO standard for finer administrative areas than state/province/territory? For example, is there a standard (ISO?) that labels the following administrative areas?:

Orange County, CA, USA
Anaheim, CA, USA (a city within Orange County)
Bengbu District, Anhui Province, China

Basically, I want a full geographic hierarchy of the entire world. Does such a dataset exist?

Things I've explored:

Natural Earth - Great open source dataset that contains good data on countries, states, and cities. It doesn't offer county/district data, though, and there doesn't appear to be any sort of standard global location ID. It also doesn't offer boundary data for cities; instead, cities are represented by a point (centroid?).
Universal Postal Union POST*CODE Database - From what I can tell, this has global postal code data, but I'm not sure where/how I can get the full dataset; all I see are a few sample files.
UN/LOCODE - This is a really extensive dataset that seems to have lots of cities and special locations (Air Force bases, airports, etc.) and even offers a unique 3-letter identifier for each, but it doesn't seem to have counties.

What else should I be looking at? Surely such a dataset exists, right?


texts: Multiple choice for table list in ArcGIS Survey123

texts:
I am working with ArcGIS Survey123 appearance. I want the question to be listed as in table form and user able to select multiple answer

I been trying to create something similar to the above in ArcGIS Survey however I ending with something like. I am not sure why

However documentation did state: The table-list appearance does not support select_multiple questions and does not support media in choice lists.
So does this mean it is some sort of limitation?


texts: Given 3 or 4 known points can I generate transformation equations to wgs-84. How would I do this.  I wont know what coordinate system the original 3-4 points are in.


texts: Arcpy - handling intermediate raster files

texts:
I have 8,760 (24*365) GRID rasters in a directory. I'd like to use use FocalStatistics analysis to reassign cells to the maximum of a 5x5 window/kernel and extract the resultant values at points specified in a points feature class.

from arcpy import CheckOutExtension
    from arcpy.sa import FocalStatistics, NbrRectangle

    if CheckOutExtension("spatial") != 'CheckedOut':
        raise Exception('Failed to check out spatial analyst extension.)

    # Testing subset of raster data. For simplicity.
    HRRR_rasters = [r'E:\HRRR_Data\2014\GRID\1 - JAN_2014\g01010000',
                    r'E:\HRRR_Data\2014\GRID\1 - JAN_2014\g01010100',
                    r'E:\HRRR_Data\2014\GRID\1 - JAN_2014\g01010200',
                    r'E:\HRRR_Data\2014\GRID\1 - JAN_2014\g01010300',
                    r'E:\HRRR_Data\2014\GRID\1 - JAN_2014\g01010400']

    neighborhood = NbrRectangle(width=5, height=5, units="CELL")
    extract_points = r'C:\Users\jking\Documents\ArcGIS\Default.gdb\sample_points'

    for raster in HRRR_rasters:
        raster_obj = FocalStatistics(raster, neighborhood, 'MAXIMUM', "DATA")

        # ------------------------Psudocode ------------------------------------
        # ExtractValuesToPoints(in_point_features, in_raster, out_point_features)

        # arcpy.da.InsertCursor to append to master table.
        # Allow raster object to be garbage collected or deleted.

The main problem that I'm having is that ExtractValuesToPoints takes a raster_layer as input rather than a raster object. I have looked into saving the raster to the 'in_memory' workspace without any luck. I've also looked into the properties of the raster object finding that raster_obj.catalogPath & raster_obj.path are None. MakeRasterLayer takes the path and name of the raster dataset which also doesn't work. I can, albeit reluctantly, save the raster to disc, use it for processing and delete it manually. I would like to avoid this write heavy approach though.


texts: Joining csv table and shapefile via ArcPy?

texts:
I want to join the attribute table of a shapefile with an externally created csv table via ArcPy.

The format of the columns are:

for attribute table: "Variable1" "Variable2" "Variable3"  "CATCHID"  "Variable 4"...

for csv table:       "CATCHID" "Results1"  "Results2"...

I want to append the attribute table with the contents of the csv table with respect to the "CATCHID" column. Since my csv table doesn't have OID's, I couldn't use JoinField_management function. Hopefully, there's another function named AddJoin_management, it is really useful when used in ModelBuilder, but it is a mess in ArcPy since it requires a table view/feature layer as an input. I made a feature layer, joined two tables and copied the feature layer to a brand new shapefile but this time the names of the columns were changed as "[attributetablename].Variable1", "[attributetablename].Variable2" etc. and they are truncated while creating the new shapefile.

Is there a way to join two tables PROPERLY using ArcPy? (which means I don't want to see the layer names predicate the variable names in the columns.)

inputShapefile = r"C:\\input.shp"
inputPath = r"C:\\input.csv"
outputPath = r"C:\\output.shp"

arcpy.MakeFeatureLayer_management(inputShapefile, "tempLayer")
arcpy.AddJoin_management("tempLayer", "CATCHID", inputPath, "CATCHID", "KEEP_COMMON")
arcpy.CopyFeatures_management("tempLayer", outputPath)

Here's the code snippet that I have problem with. This join operation produces a new shapefile that includes "[input].CATCHID" as a column header, for instance. However I need just "CATCHID" as the header.


texts: How can I measure the distance between many polygon features and either another polygon or a line in QGIS?

texts:
I often have to manually measure the distance between a single site and other sites around it, but obviously this can get long-winded and I feel like there must be a much more efficient way of doing it using QGIS.

I need to be able to measure the distance between the closest points of each feature within a polygon vector layer to my site boundary, which is also normally a polygon layer.


texts: Controlling transparency for basemap layer in ArcGIS for Android?

texts:
The options seem sparse, but maybe I'm missing something (I'm new to both ArcGIS for Android and Android Honeycomb). Is there a way to make the basemap transparent?

Update: this question is on the Esri application available in the android market.


texts: PostGIS query works as SELECT but not UPDATE

texts:
I want to clip the geometries of my table by a polygon.
I have a query that does this just fine as a SELECT query:
WITH
    clip AS
    (
        SELECT ST_GeomFromEWKT('SRID=26915;POLYGON ((379308.073398988 5232934.49102113,377638.690726783 5102722.64258913,537899.427258473 5097714.49457251,541238.192602883 5236273.25636554,379308.073398988 5232934.49102113))') AS geom
    ),
    new_geom AS
    (
        SELECT
            input.gid,
            CASE
                WHEN ST_CoveredBy(input.geom, clip.geom) THEN input.geom
                ELSE ST_Multi( ST_Intersection(input.geom,clip.geom) )
            END AS geom
        FROM "rail_lines" AS input
        INNER JOIN clip
        ON ( ST_Intersects(input.geom, clip.geom)
        AND NOT ST_Touches(input.geom, clip.geom) )
    )
SELECT * FROM new_geom;

However, when I use this same query to UPDATE the geometry of the input table the query either does not complete or does nothing (ie. does not alter the geometries of the input table):
WITH
    clip AS
    (
        SELECT ST_GeomFromEWKT('SRID=26915;POLYGON ((379308.073398988 5232934.49102113,377638.690726783 5102722.64258913,537899.427258473 5097714.49457251,541238.192602883 5236273.25636554,379308.073398988 5232934.49102113))') AS geom
    ),
    new_geom AS
    (
        SELECT
            input.gid,
            CASE
                WHEN ST_CoveredBy(input.geom, clip.geom) THEN input.geom
                ELSE ST_Multi( ST_Intersection(input.geom,clip.geom) )
            END AS geom
        FROM "rail_lines" AS input
        INNER JOIN clip
        ON ( ST_Intersects(input.geom, clip.geom)
        AND NOT ST_Touches(input.geom, clip.geom) )
    )
UPDATE rail_lines
SET geom = new_geom.geom
FROM new_geom
WHERE new_geom.gid = rail_lines.gid;

Does anyone know what is happening here?
I'm also open to other (simpler) means to clipping a table's geometry in-place without creating a new table.
Update
I'll update my question based on the suggestion by @JohnPowellBarca.
Updating the table still doesn't work but it might help to illustrate the issue.
Here's the query I'm using. First, I create a new table based on the intersection with my area of interest. Then I use that new table to attempt to update the geometry of the original table.
CREATE TABLE rail_lines_clip AS
SELECT in_table.gid, in_table.geom
FROM
    rail_lines AS in_table,
    aoi AS clip
WHERE ST_Intersects(in_table.geom, clip.geom);

UPDATE rail_lines
SET geom = rail_lines_clip.geom
FROM rail_lines_clip
WHERE rail_lines_clip.gid = rail_lines.gid;

And here's a visual of the result after running the query:

The area of interest (ie. aoi) is in black.
The new table (ie. rail_lines_clip) is in green.
And the original table (ie. rail_lines) is in red.
From my understanding, the geometry of rail_lines should be identical to that of rail_lines_clip.
Maybe I should just use that new table, and then drop the old table, and rename the new to old.  Such as:
CREATE TABLE rail_lines_clip AS
SELECT in_table.gid, in_table.geom
FROM
    rail_lines AS in_table,
    aoi AS clip
WHERE ST_Intersects(in_table.geom, clip.geom);

DROP TABLE IF EXISTS rail_lines;

ALTER TABLE IF EXISTS rail_lines_clip
RENAME TO rail_lines;

(Which works as expected, by the way).
But, I'm still curious: Has anyone else had success updating the geometry of a table in place?


texts: Calculating TWI after getting the flow direction, flow accumulation and slope

texts: There are many tutorials online with this, but it seems that people use different equations and/or what they put into the general ln(a/(tanb)) equation. Do I need to know the cell size of the raster DEM to complete the equation, and if so, where does it go?


texts: How can I read MrSID files in Python (using gdal)?

texts:
I have to read a large number of MrSID images and process them in Python. I usually access GIS file formats using the gdal application and the gdal Python library. However, the gdal package for Python does not come with MrSID support (I use the portable WinPython and the precompiled gdal package from here). Then I also use QGIS that comes with gdal support: here, gdal does support MrSID out of the box (command: gdalinfo.exe --formats).

I was wondering if there is a gdal package for Python that does support MrSID out of the box, or do I have to compile it from the source (using Windows 8.1 64-bit).

I tried to just copy the QGIS dll "...\gdalplugins\gdal_MrSID.dll" to the same folder in my Python gdal library, but this makes the command "gdalinfo.exe --formats" to stop with an error.

This is the code that I use in Python to see if gdal supports MrSID:

import gdal

gdal.AllRegister()
for i in range(1, gdal.GetDriverCount()):
    drv = gdal.GetDriver(i)
    print drv.GetDescription(),

Using this after implementing Kersten's solution, the list does not include the entry "MrSID". Calling the command "gdalinfo --formats" however, does list the entry "MrSID (rov): Multi-resolution Seamless Image Database (MrSID)".


texts: Landsat / MODIS spatio-temporal fusion: STARFM and software implementations

texts:
My main goal is to create Landsat/MODIS image pairs, perform an spatial-temporal fusion which leads to synthetic composites on dates between Landsat acquisition dates.

This brought me to the paper by Gao, and the STARFM software. However, the only software implementation of STARFM or derivatives seems to be the executable from the USDA (https://www.ars.usda.gov/research/software/download/?softwareid=432).

I cannot find any other software implementation of STARFM or derivative, neither on GitHub, R/Python, or elsewhere. Does anyone know of any other software implementations of these Spatio-temporal fusion models?


texts:
Is there a way to add a combobox to a toolbar for QGIS? So far I can only find examples that add buttons.

I found this tutorial which seems to explain how to do this with general PyQT but haven't had any luck at all doing this with QGIS. So any help or guidance would be greatly appreciated.


texts: Earth Engine: feature collection will not export even in chunks

texts:
I am trying to export zonal statistics from the Worldpop dataset for about 450,000 0.01° tiles, but I can't get it to work. When I break the data into chunks of 50,000 observations, the computation still times out on each.
My tiles were imported into GEE (anyone can read):
var tiles = ee.FeatureCollection('users/salem043/deforSample');

I mosaicked the country-year data in Worldpop together to produce annual images:
var years = ee.List.sequence(2000, 2012)
// Make some band names from the years numbers.
var yearBandNames = years.map(function(n) { return ee.Number(n).format("pop_%d") })

var yearImages = ee.ImageCollection(
  years.map(function(year) {
    var oneYear = ee.ImageCollection("WorldPop/GP/100m/pop")
      .filter(ee.Filter.calendarRange(year, year, 'year'))
      .filterBounds(tiles)
    return oneYear.mosaic()
  }))
  .toBands()
  .rename(yearBandNames)

Then to produce the data frame, I select the year image and perform a reducer function to get the zonal statistics. I perform the same function as defined below for all years 2001-2012:
var pop_2001 = yearImages.select("pop_2001");

var sumPop_2001 = tiles.map(function(feature) {
  return feature.set(pop_2001.reduceRegion({
    reducer: ee.Reducer.sum(), geometry: feature.geometry(), scale: 30
  }));
});

I tried exporting this data in chunks of 50,000 tiles at a time, but each computation timed out. Can anyone offer advice as to how to make this program more efficient so that it exports? Full code here


texts: Get ArcGIS Pro to show which feature is causing copy/paste to fail

texts:
I have an Enterprise geodatabase feature class in a feature dataset that I need to copy over to a gdb. I've tried

right click to copy and paste,
the copy tool in ArcGIS Pro,
copy management in ArcPy.

All raise the same error:

SG_POLYGON_HAS_VERTICAL_LINE (-2033)     Polygon shell contains a
vertical line

I'm having trouble tracking down, which feature is causing the error so I can fix it. It's an Oracle db and I've tried to run validate geom with context and got no results that stand out. None of the features are 3D and the feature class does not store z values.
I have no idea how this feature made it into the db to begin with. All features are loaded to the db using copy paste so I would think this would have been caught then, but alas it did not and I am at a loss.
I'm wondering if anyone knows a way to get the tool to print the name of the feature it's trying to copy so I can get a hint of where to look.
EDIT: Wanted to add that I have tried to use the export tools as well. They run without error, but less than half of the features end up being exported to the new geodatabase.


texts: R: Using st_intersects to classify points inside, outside and within a buffer

texts:
I have two objects, a shapefile of location points, points, and a shapefile of polygons, world_buffer (both +proj=utm +datum=WGS84). The world_buffer is the world's coastlines, found here, with a 5 km buffer around it, computing using gBuffer().
I am trying to classify the "points" as either within the world polygons (land), within the buffer (coast), or outside of the buffer (ocean).
I have tried using st_intersects() to get the points within the buffer:
library(sf)
intersect <- st_intersects(points, world_buffer, sparse = TRUE, prepared = FALSE)

But this returns a large sgbp that has two columns, a column with an intersect number and then a column full of integer [0].


texts: Populate b.field into a.field when 2 vector points not always coincident but should be coincident

texts:
Using QGIS/pgAdmin3, I have a situation where two vector points, should be coincident but I can't guarantee they always are, I need to take feature a.field and populate a field in feature b, upon both insert and update to the geometry, update being let's say feature a has been moved spatially to be coincident to a new/different b feature.

A lot of my attempts have resulted in this error:

Could not commit changes to layer term

Errors: ERROR: 1 feature(s) not added.

  Provider errors:
      PostGIS error while adding features: ERROR:  more than one row returned by a subquery used as an expression
    CONTEXT:  SQL statement "update term set
    station_id = (select id from struct
    where
    st_within(term.geom, st_buffer(struct.geom, 0.1)))"
    PL/pgSQL function term_station_id_update() line 5 at SQL statement

I am trying to create the right function to use with insert and update triggers, function attempt is as follows:

update term set

station_id = (select id from struct

where

st_within(term.geom, st_buffer(struct.geom, 0.1)));"

I have also tried similar functions with st_dwithin, st_within, and just with the geometry being = between the two layers if truly coincident, does not work however, similar error messages return.


texts: Adding a Tiles/wms? (openStreetMap) layer with python in QGIS

texts:
In my QGIS python plugin I want to add a background map (if the user haven't added one), hence I think I missing something in the url for the Tiles/wms service.. (the function returns invalid layer) Any suggestions?

    sources = [layer.source() for layer in QgsProject.instance().mapLayers().values()]
    print(sources)
    source_found = False
    for source in sources:
        if 'xyz&url' in source:
            source_found = True
            print('found')
    if not source_found:
        print('adding')
        urlWithParams = 'type=xyz&url=http://a.tile.openstreetmap.org/%7Bz%7D/%7Bx%7D/%7By%7D.png&zmax=19&zmin=0&crs=EPSG3857'
        rlayer = QgsVectorLayer(urlWithParams, 'OpenStreetMap', 'wms')
        if rlayer.isValid():
            QgsProject.instance().addMapLayer(rlayer)
        else:
            print('invalid layer')


texts:
I georeferenced a map and it posted to canvas as a mirror image of the original map. I set my project CRS before starting and also made sure to set the same SRS for the Geotiff being created and saved. I made sure that my Eastings and Northings are correct.

Please see that Texas is on the wrong side of Louisiana, all text from the map is mirror imaged from the original way the map was georeferenced. North is still pointing north. I am using Version 2.82 and I entered the latitude longitude values in dd mm ss.ss

Do you know what I am doing wrong?

I re-attempted the georeference by adding points in clockwise motion, starting with top right corner. Output shows that georeferenced points & the map itself are mirror imaged. CRS is NAD27 Louisiana South EPSG: 26782. Transform type: Thin Plate Spline, Resample method: Near neighbor, Compress: None


texts: Is there a way to edit a column in an attribute table at once within an Arc geodatabase when there is subtype and coded domain?

texts:
In a regular shapefile with no coded properties, I can use 'Field Calculator' on my selected features to change all attribute rows at once. However with subtypes and coded domains there are drop-downs.

Is there a script I can run or similar to do this batch at once. There are 500+ records I need to change, and that would be very tedious work.


texts: How to check are there any OGRGeometry (polygons) within an OGRGeometry (polygon) ?

texts:
Please see my image:

I used the below statements to get Geometry:

OGRGeometry *poGeometry;

poGeometry = poFeature->GetGeometryRef();

Do you know how to check whether the Polygon 1 (geometry) is empty (or include a small polygon) or not ?

If it is possible, how to get or identify the polygon 2 ?

Thanks and regards,

Tai


texts:
New to laspy and encountering an error. The first code block describes the file header, the second (In[24]) fails to open the file with laspy:

In [21]:
import numpy
import laspy
from laspy.file import File

In [22]:
import struct
headerstruct = ( ('filesig', 4,'c',4) ,
         . . .  omitted for brevity, see result for fields . . .
                   ('zmin'         , 8,'d',1) )

def parseHeader(filename, verbose=True):
    # courtesy Phil Wilkes
    fh = open(filename,'rb')
    header = {'infile':filename}
    with open(filename, 'rb') as fh:
        for i in headerstruct:
            if i[2] == 'c':
                value = fh.read(i[1])
            elif i[3] > 1:
                value = struct.unpack( '=' + str(i[3]) + i[2] , fh.read(i[1]) )
            else:
                value = struct.unpack( '=' + i[2] , fh.read(i[1]) )[0]

            if verbose:
                print i[0] + '\t', i[2] + '\t', value

            header[i[0]] = value

In [23]:
filename = r'./sf_shared/40095_Flightline_002_Z.las'
parseHeader(filename, verbose=True)

filesig   c   LASF
filesourceid  H   0
reserved  H   0
guid1 L   0
guid2 H   0
guid3 H   0
guid4 B   (0, 0, 0, 0, 0, 0, 0, 0)
vermajor  B   1
verminor  B   2
sysid c
gensoftware   c   TerraScan
fileday   H   349
fileyear  H   2011
headersize    H   227
offset    L   229
numvlrecords  L   0
pointformat   B   1
pointreclen   H   28
numptrecords  L   8917629
numptbyreturn L   (6711217, 2002830, 196359, 7135, 88)
xscale    d   0.01
yscale    d   0.01
zscale    d   0.01
xoffset   d   -0.0
yoffset   d   -0.0
zoffset   d   -0.0
xmax  d   606205.36
xmin  d   605000.63
ymax  d   6599629.15
ymin  d   6596922.71
zmax  d   129.27
zmin  d   43.74

In [24]:
infile = File(filename, mode='r')  #laspy.file.File

---------------------------------------------------------------------------
LaspyException                            Traceback (most recent call last)
<ipython-input-24-19025bf517df> in <module>()
----> 1 infile = File(filename, mode='r')

/usr/local/lib/python2.7/dist-packages/laspy/file.pyc in __init__(self, filename, header, vlrs, mode, in_srs, out_srs, evlrs)
     60         self.out_srs = out_srs
     61
---> 62         self.open()
     63
     64     def open(self):

/usr/local/lib/python2.7/dist-packages/laspy/file.pyc in open(self)
     71             ## Make sure we have a header
     72             if self._header == None:
---> 73                 self._reader = base.Reader(self.filename,mode= self._mode)
     74                 self._header = self._reader.get_header()
     75             else:

/usr/local/lib/python2.7/dist-packages/laspy/base.pyc in __init__(self, filename, mode, header, vlrs, evlrs)
    151         self.padded = False
    152         if self.mode in ("r", "rw"):
--> 153             self.setup_read_write(vlrs, evlrs)
    154             return
    155

/usr/local/lib/python2.7/dist-packages/laspy/base.pyc in setup_read_write(self, vlrs, evlrs)
    162     def setup_read_write(self, vlrs, evlrs):
    163         self.data_provider.open("r+b")
--> 164         self.data_provider.map()
    165         self.header_format = Format("h" + self.grab_file_version())
    166         self.get_header(self.grab_file_version())

/usr/local/lib/python2.7/dist-packages/laspy/base.pyc in map(self)
     93                 raise LaspyException("Invalid Mode: " + str(self.mode))
     94         except(Exception):
---> 95             raise LaspyException("Error mapping file.")
     96
     97     def remap(self,flush = True, point_map = False):

LaspyException: Error mapping file.

Unfortunately, the laspy documentation leaves you at "When a file is opened in read mode, laspy first reads the header, processes any VLR and EVLR records, and then maps the point records with numpy. If no errors are produced when calling the File constructor, you’re ready to read data!"  But it doesn't suggest what to do if there are errors!

Can you see an error in my code or suggest the source of the problem?


texts: Batch flip direction of selected lines in ArcGIS 10.3

texts: Is there a way to flip the direction of all selected lines at once?


texts: Adding/Removing GeoJSON layers using Leaflet

texts:
I am trying to show different GeoJSON layers at different zoom layers using the Leaflet API.  I can load and display all three layers at once (though I do not actually want them all to show at once).  I can load and display them at different zoom levels.

I have the code set up so that at Zoom levels 1-6, the map will show one GeoJSON layer.  At levels 7-10, it will show another, and at levels 11+ it will show a third.  Displaying them works.  What I am trying to get working now is turning the others off if one is displayed.  Going from 1-6 to 7-10 works (meaning that it turns off the 1-6 layer correctly), but not from 7-10 to 11+ (meaning that the 7-10 layer sticks around) and I cannot figure out why (it uses the same code).

Here is the ajax for the GeoJSON layers:

function getJson(defaultStyle, map, simp, geojsonLayer){
var url = 'file' + simp + '.json';
map.removeLayer(geojsonLayer);
geojsonLayer.clearLayers();
$.getJSON(url, function(data){
    geojsonLayer = L.geoJson(data, {
        style: defaultStyle,
        onEachFeature: onEachFeature
    });
    geojsonLayer.addTo(map);
});
}

And here is the main function that calls the ajax depending on the zoom. simpCounter is set to 0 initially.

map.on('zoomend', function(e) {
if (map.getZoom() >= 7 && map.getZoom() <= 10) {
    if (simpCounter == 0 || simpCounter == 2) {
    getJson(defaultStyle, map, 60, geojsonLayer);
    simpCounter = 1;
    }
} else if (map.getZoom() >= 11) {
    if (simpCounter == 0 || simpCounter == 1) {
    getJson(defaultStyle, map, 35, geojsonLayer);
    simpCounter = 2;
    }
}
});

So again, the first transition turns off the old layer correctly, the but second transition does not.


texts: Calculate nearest town for 2900 lot centroids

texts:
I have 16 town point shapefiles.
I have a layer of 2900 centroids for lots.

I need to determine the nearest town for each of the centroids.

The option I came up with is to Split by Attribute by name, thus creating 2900 point layers for the centroids, then doing the Near tool for each centroid. This will take a long time and may even crash ArcMap. There is probably a script for it, but I am very new to python.


texts: can Well Known Text tables be used effectively, stored in Access or Base?

texts:
My question is:  How do I import WKT text data back into QGIS, and how useful could Base or Access be in storing the data using WKT rather than spatially enabling it within a proper GIS database, like Spatialite?

I get an error when I use the "Add Delimited Text Layer" (within the Manage Layers toolbar) in QGIS 1.8 on Windows XP, to add data with multipolygon geometries from a WKT field.  I select the WKT[] and Delimiters[] radio buttons.  Before I press [OK], I check that the Sample text box shows the columns parsed correctly.

A CSV file with commas scrambles the columns together by the tool, even though Excel reads it perfectly (it was created using QGIS, exporting a Spatialite table the OGR LAYER option to add WKT geometries column).

And so I used Excel to SAVE AS a TAB delimited .txt file.  The txt file, TAB delimited, pre-processing sample text shows that the data is being parsed into the proper columns perfectly.  This looks promising, then I hit [OK] to process it and a dialog box pops up reporting it failed, "The layer...is not a valid layer and can not be added to the map [OK]"

Currently I'm using Spatialite with great success.  But I'm experimentally trying to learn how to create and use WKT text files so that I can perhaps import them into Access or Base, to see how feasible it is to deliver data in alternative, consumer level database formats.  The users would have experience with them, but I currently don't know very much about either.


texts: WFS Transaction Post Method 403 OpenLayers 3 + GeoServer

texts:
I am having difficulty understanding the topic when it comes to XML generated from GeoServer and how can we interact with it. I want to to implement WFS transaction service and updating my attributes via web browser.

.

I have the snippet code from book title Mastering OpenLayers about WFST

        map.getLayers().item(1).getSource().on('change', function (evt) {
        if (this.getState() === 'ready') {
        map.getLayers().item(1).buildHeaders();
        var geoJSONSerializer = new ol.format.GeoJSON();
        var featString = geoJSONSerializer.writeFeatures(this.getFeatures());
        var request = new XMLHttpRequest();
        request.open('POST', 'http://127.0.0.1:8080/geoserver/wfs');
        request.send(featString);

        var WFSTSerializer = new ol.format.WFS();
        var featObject = WFSTSerializer.writeTransaction(this.getFeatures,
        null, null, {
            featureType: 'admin_kel_batu',
            featureNS: 'http://sigweb_cafe',
            srsName: 'EPSG:4326'
        });
        var serializer = new XMLSerializer();
        var featString = serializer.serializeToString(featObject);
        var request = new XMLHttpRequest();
        request.open('POST', 'http://127.0.0.1:8080/geoserver/wfs');
        request.setRequestHeader('Content-Type', 'text/xml');
        request.send(featString);
        console.log(featString);
    }
});

But still I got error and I have no idea how to fix it. What am I doing wrong or what is missing?

this is the xml being sent :

<Transaction xmlns="http://www.opengis.net/wfs" service="WFS" version="1.1.0" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://www.opengis.net/wfs http://schemas.opengis.net/wfs/1.1.0/wfs.xsd"/>


texts: How to retrieve historical imagery from web mapping services?

texts:
Is there a way to retrieve historical imagery from Google/Bing/Yandex mapping services through their APIs?

For example, I'd like to show on my website imagery from 2005.


texts: ArcMap symbols exported to PDF don't draw when viewed on another computer

texts:
I have maps of highway data. They include highway shields to identify highways. When I export them to PDF they look fine on my computer, but when viewed on computers outside of work, the highway shields turn into letters. See below

This is how the highway shields look on my work computer.

This is the shields viewed on another computer. This image is not the same location geographically, but I think you get the picture.

The shields themselves are Marker Text Backgrounds of Character Marker Symbol type. The blue and yellow symbols, for example, are made of 3 pentagons laid on top of one another.

My solution for this is to check "Convert Marker Symbols to polygons" in the export to PDF commands. This makes the shields viewable on other computers, but also causes them to come out a little squiggly. It's not bad enough to make the solution unviable, so I won't share that image or ask about how to fix it. So I'm not actually really looking for help fixing my maps, since I have an acceptable solution.

Am I right in saying that Adobe is actually trying to populate the maps with other files stored locally?

I would have thought everything was in the map file itself.

If I'm wrong, what is actually causing this error?


texts:
I need to publish WMS service with GetFeatureInfo but I want to hide some fields.  I've made them not visible using layer properties in ArcMap but it resolved the problem partially. When I add WMS service to our www geoportal I can't see hidden field in html popup. It's correct. But when I connect with WMS service through ArcCatalog and  add the service to ArcMap, I can see still all the fields in Identify window.

Is it possible to hide the fields entirely by modifying the GetFeatureInfo XSL templates?


texts: Dealing with multiple input layer in QGIS modeler

texts:
I would like to add atrribute with layer name using field calculator for many input layers and later merge updated layers into one.
I am able to use multiple input layer in modeler and later merge them, but I have no clue how to add field calculator algorithm before as it seems not to allow multiple input.

I need field calculator on place of red square to add layer name before merging.
There is also issue as data types of inputs are not fully consistent for ID (Int and Int64) taht is not passing merge tool, but merge from MMQGIS plugin is dealing with it alright. Plugin algorithms are not available in modeler.


texts: Getting the imageID within a set of images (imageCollection) using the targeted coordinate

texts:
I am trying to download a high solution satellite image of a targeted area (roughly 30m * 30m) with (-80.562459, 43.470784) as the center.
This is my approach.
I am accessing the skysat dataset:
https://developers.google.com/earth-engine/datasets/catalog/SKYSAT_GEN-A_PUBLIC_ORTHO_RGB#image-properties
This is the code in the code editor:
var dataset = ee.ImageCollection('SKYSAT/GEN-A/PUBLIC/ORTHO/RGB');
var rgb = dataset.select(['R', 'G', 'B']);
var rgbVis = {
      min: 11.0,
      max: 190.0,
    };
Map.setCenter(-80.562459, 43.470784, 24);
Map.addLayer(rgb, rgbVis, 'RGB');
Export.image.toDrive({
      image: rgb,
      description: 'imageToDrive',
      scale: 30,
      region: geometry
    });

geometry here is a hand-drawn rectangle. Then I run the code. However, the image is not downloaded correctly. Error message is: "Error: Image.clipToBoundsAndScale, argument 'input': Invalid type. Expected type: Image<unknown bands>. Actual type: ImageCollection."
So I think the solution is to get the image of the targeted area using the coordinate, (-80.562459, 43.470784). But how to find out the image ID of this satellite image?


texts:
I have MEDIA column in my attribute table contains video URLs (YouTube iframe & HTML5 video tag).

By Clicking on feature on map popup display and play the video.

I would like to stop video play when I click on close button of esriPopup window (InfoWindow).

I know the hide event of infoWindow, but How do I stop play video on close of the InfoWindow

 app.map.infoWindow.on('hide', function () {
var info = app.map.infoWindow.features[0].attributes.MEDIA;
});

Is there a way to stop video on close event of infoWindow?


texts: Generating points with specific distance properties

texts:
I have a layer with multiple lines. I want to generate points at each intersection which is easy with the Geoprocessing tools. The next step however is giving me issues;

From each intersection of two or more lines I want to generate points with equal distance to the intersection on all lines intersecting.

I can create the point (x,y) but want to make sure that (x_1, y_1), (x_2, y_2). (x_3, y_3), (x_4, y_4) are the exact same distance from the intersection point no matter the shape length.

Is there any smart way of doing these steps?


texts:
I have two datasets:

A point shapefile of companies (each row/record is a unique company)
A list of awards given to those companies (each row is a unique award, but since companies can win multiple awards, the same company can exist in multiple rows).

I want to create a Carto map where each point on the map represents one company, but when you click on that company point's pop-up, it shows all the awards that that company has won.

I've tried multiple variations of the "Link second layer" and "Add Columns from Second Dataset" analysis options, but I'm not getting the result I want. I'm thinking that the join created in "Add Columns from Second Dataset" is just a one-to-one join, and it's only joining on the first record it comes across that has a shared company, OR the pop-up itself is restricted to showing only one award.


texts: Creating local accurate projection for proj4?

texts:
We recently discovered that the UTM projection we have been using is not accurate enough for our application when doing calculations that require distances.

Instead of rewriting our code base to use great circle calculation I would like to understand how proj4 projections work but havent been able to google what i needed. Possible because I dont yet know exactly what to google for.

I have a lat,lng of a point that I want to place my local projection around and the area it should cover is a radius of 25km. How would I create a proj4 string for this local projection with units of meters.


texts: How to list data in the attribute table in the same order as I selected the data?

texts:
I have a map that consists of a number of pipes that are connected to each other. Is there a way I can list data in the attribute table in the same order as I selected the data?

For example: I select first pipe1, then pipe5, then pipe7 and then in the attribute table the first item on the list is pipe1, second pipe5 and third pipe7.


texts: In QGIS combine multiple polygon parts of a single feature

texts: This is a repost of my unsolved QGIS dilema. The image shows a single selected feature that consists of two polygons and the corresponding attribute table. Some nodes of the left poly have been added and moved to illustrate the boundary/ overlap problem. Before the nodes were moved and added the two polys were snapped together and I could not get rid of the line between them. Hovering over the nodes, pop ups identify the polys as either 0 or 1. How does one create a single polygon from the two?


texts:
I am iterating through ~300 raster tiles, using a polygon feature class to extract by mask. The model successfully processes the first ~100 tiles and then I get the following error on the same tile each time:

  ERROR 999999: Error executing function.
  The table was not found. [VAT_T2001]

  ERROR 010092: Invalid output extent.
  Failed to execute (ExtractByMask).

I changed the extent parameter of the extract by mask tool from default to "intersection of inputs" and got the same error for the same tile. The tool seems to work fine when I change the extent parameter to "union of inputs" but it takes an extremely long time to process each tile.

I thought that the error might be occurring because there is no overlap between the error tile and the mask feature class, but I confirmed that the tool successfully processes other tiles where there is no overlap.

The tiles and the feature class have the same spatial reference. The error also occurs for a handful of other tiles in the data set.

Does anyone know how to resolve these errors while maintaining efficient processing?

The model is just an iterator sitting in front of the extract by mask tool. So this is a question about why extract by mask works for some raster datasets and not others when the datasets have the same attributes and the tool has the same parameters.


texts: QGIS3 - Clipping polygons to a hemisphere for azimuthal orthographic projection

texts:
I'm a bit of a novice with QGIS, so this is really throwing me for a loop.

For the maps I'm intending on making, an azimuthal orthographic projection would be my preference.

For example, one map I'm trying to recreate is this map of the Pacific islands from Wikipedia:

I've been using the NaturalEarth Quick Start package, which has been a really good launching-off point, but I'm running into some issues at the horizon of the orthographic projection.

I've tried adding a custom projection:

+proj=ortho +lat_0=0 +lon_0=170 +x_0=0 +y_0=0 +a=6371000 +b=6371000 +units=m +no_defs

But that gives the following map:

Which certainly is close to being right, but the ocean is obviously being cut off in the Arctic and on the southern horizon.

I understand the issue is that naively mapping layers described in WGS 84 to an orthographic projection leaves you with artifacts at the horizon, but how do I fix it?

I first tried to use the Clip to Hemisphere plugin described in Ortho Projection produces artifacts, but it seems incompatible with QGIS 3; when I try adding the plugin, QGIS3 gives an error message and places it in the "Invalid" section of the Plugins window.

I then tried both solutions presented at Where did the polygons go after projecting a map in QGIS? Again the issue is with QGIS 2/3 incompatibility.

For Solution 1, both Numerical Digitize and CAD Tools are again incompatible with QGIS 3 and give the same result as Clip to Hemisphere.

For Solution 2, I can't disable on-the-fly-reprojection in settings since that setting was removed in QGIS 3.


texts: Performing unsupervised classification using ArcGIS Desktop?

texts:
How to distinguish between water bodies and clouds using unsupervised classification in sentinel 2?

When i used NDWI and NDSI both are same objects.


texts:
As a part of my university project I am developing a web service that provides information about coordinate reference systems that were created in USSR.

It will show its Proj.4 and WKT definitions, history of creation and use. Also there will be a visualization of base layer and correct graticule for the user chosen projection.

I need to find some kind of similar services to make comparison. Maybe borrow good features and avoid doubtful ones in my project.

All I could find by so far are epsg.io, spatialreference.org, en.wikipedia.org/wiki/List_of_map_projections mand QGIS CRS project property option window.

Are there any other similar web/desktop services of that kind?


texts: I cannot use Point entry into Graphical Modeler in QGIS

texts:
I would like to use a Point entry to get the point coordinates  in the map to add in the network analysis as a start point.
I have added the Point entry into Graphical Modeler, I can click on the map but I doesn't get any coordenates.
When I use the analysis network tool "Shortest path (point to layer)" I haven't any problem.


texts: Inaccurate difference operations in GEOS