texts: Warning: skipping unsupported geometry type 'Compound Curve' when using v.transform

texts:
I am trying to transform data in QGIS from a DXF file.
I keep getting the error:

Warning: skipping unsupported geometry type. Compound Curve

when I do a v.transform.
Eventually, I want to store the data in a shapefile.


texts:
I'm trying to read a folder containing ~80,000 .kmz files of vector data using sf::st_read.

The top-level folder contains a single large .kml file that I expect ties all of the other folders/data together.

These are arranged into a group of OSGB folders, e.g. folder  contains four sub-folders    and  and four .kmz files that I presume describe the folder contents.

Each of these folders then contains a cluster of .kmz files where I presume the vector data is contained for each grid square (e.g. )

So the structure is:
1x top-level.kml

~400x  grid square file.kmz

~400x  grid square folder

~100x     smaller grid square files.kmz

Or...do I simply take all of the small grid square .kmz files, unzip them and import each 1-by-1?
I've not yet tried this in the hope that there's a more elegant solution somewhere

Apologies, if I've missed any earlier advice about this but I've searched StackExchange for a solution to this)
Further apologies but the data isn't public domain so I'm not sure how best to share, either!


texts: pyqgis: AttributeError: 'QgsRasterLayer' object has no attribute 'drawingStyle'

texts:
I am following http://www.qgis.org/en/docs/pyqgis_developer_cookbook/intro.html. In python console after a raster layer is loaded to a variable "rlayer"

>>>rlayer
  <qgis.core.QgsRasterLayer object at 0x113faadf8>

>>>rlayer.drawingStyle()
  Traceback (most recent call last):
         File "<input>", line 1, in <module>
         AttributeError: 'QgsRasterLayer' object has no attribute 'drawingStyle'

>>>rlayer.setDrawingStyle(QgsRasterLayer.SingleBandPseudoColor)
  Traceback (most recent call last):
         File "<input>", line 1, in <module>
         AttributeError: type object 'QgsRasterLayer' has no attribute 'SingleBandPseudoColor'

>>>rlayer.setColorShadingAlgorithm(QgsRasterLayer.PseudoColorShader)
   Traceback (most recent call last):
         File "<input>", line 1, in <module>
         AttributeError: 'QgsRasterLayer' object has no attribute 'setColorShadingAlgorithm'

what did I do wrong?


texts: Android App to help digitizing point layers on a Tablet

texts:
In QGIS, using a pen on my touchscreen PC I can I can digitize a lot of points (e.g. 14000 pts in 2 hours) faster than with a mouse: without having the attribute form to show up I simply need to touch the screen where the point should be rather than moving the pointer with the mouse and clicking.

Unfortunately, I have a collaborator that does not have a touchscreen PC.

Does anybody knows if:

Is there any app for android tablets that can be used to achieve
the same goal?
It is possible in some way to take advantage of the touch screen
of an Adroid tablet to help digitizing?


texts:
How can I use Geotrellis (distributed format) to Rasterize LiDAR (or XYZ DATA ) in HDFS?

I converted las file to spark RDD with the IQMULUS library and I want to use this RDD ( [x1,y1,z1] [x2,y2,z2] ....) to create a DEM with Geotrellis.


texts:
I have my OSM data in a local server.
I already done with the contour lines using Phyghtmap and achieve a good result.

Now I'm trying to show the relief "shadows" like this beautiful map:

I know it is not OSM data, just some kind of transparent PNG layer, but I don't know how to create it...

EDIT

By reading this I'm now able to convert my HGT files downloaded by Phyghtmap  to GeoTIFF but it have a very bad resolution...


texts: How to Track Individual Layer Usage Within Map Service?

texts: I have a public-facing web mapping application (flex-based), and I need to get the usage/statistics on the layers that people are using, however, a lot of these layers are one of many within a map service (e.g. "Utilities" is a service that has many layers within it).  Besides splitting them apart into separate services, what are some ways I can go about tracking this information?  As far as I know, ArcGIS Server stats only show usage by service (not the layers within) and Google Analytics only displays information on page views/widget usage.


texts: Overpass-API request does not get executed completely : particular results lead to mistakes

texts:
I have some issues with the following request via overpass-turbo.eu. Only the last expression works. In other words, only the very last expression will be performed and all the others are ignored.
[out:csv(::id,::type,::lon,::lat,amenity,name,"addr:postcode","addr:city","addr:street","addr:housenumber","contact:website",website,"contact:email")]
[timeout:600];
rel[boundary=administrative][admin_level=6][name="München"] -> .city;
(nwr[amenity=hospital][name](around.city:2000);
nwr[amenity=school][name](around.city:2000););
nwr[amenity=childcare][name](around.city:2000);
nwr[amenity=kindergarten][name](around.city:2000);
nwr[amenity=nursing_home][name](around.city:2000);
nwr[amenity=college][name](around.city:2000);

out center;

In fact only this term will be executed:
nwr[amenity=college][name](around.city:2000);

This is given back:
@id @type   @lon    @lat    amenity name    addr:postcode   addr:city   addr:street addr:housenumber    contact:website website contact:email
2354666181  node    11.5557656  48.1274889  college EAM-Europa Akademie München - School Of International Business  80337   München Adlzreiterstraße    6
3324407853  node    11.5551400  48.1264588  college Akademie Deutsche POP   80337   München Adlzreiterstraße    18a     http://www.deutsche-pop.de/standort/muenchen
3324407854  node    11.5713030  48.1285751  college Akademie Deutsche POP   80469   München Klenzestraße    67      http://www.deutsche-pop.de/standort/muenchen
4379246720  node    11.5629942  48.1469822  college Hochschule für Politik München  80333   München Richard-Wagner-Straße   1
4448465264  node    11.4461165  48.1488143  college DEKRA Akademie GmbH                     https://www.dekra-akademie.de/de/muenchen/
4984543723  node    11.5998193  48.1353510  college Münchner Volkshochschule – Einstein 28  81675   München Einsteinstraße  28      https://www.mvhs.de/ueber-uns/standorte/gasteig-einstein-28-stadtbereich-mitte/
5542020160  node    11.6368555  48.0790856  college Fakultät für Betriebswirtschaftf
5542020161  node    11.6348063  48.0806016  college Fakultät für Betriebswirtschaft
5542020162  node    11.6336959  48.0810352  college Fakultät für Betriebswirtschaft
18940172    way 11.4648639  48.1127381  college Haus A Genzentrum der LMU und Institut für Biochemie    81377   München Feodor-Lynen-Straße 25      https://www.genzentrum.lmu.de/
19329687    way 11.5817581  48.1500346  college Herzogliches Georgianum 80539   München Professor-Huber-Platz   1       http://www.herzoglichesgeorgianum.de/Herzogliches_Georgianum/Home.html
27959435    way 11.5929977  48.1424114  college Bayerische EliteAkademie    80538   München Prinzregentenstraße 7       http://www.eliteakademie.de
44615795    way 11.5806579  48.1550310  college Priesterseminar St. Johannes der Täufer                 http://www.priesterseminar-muenchen.de/
166433609   way 11.5888108  48.1304498  college Collegium Carolinum
867729755   way 11.6027411  48.1324959  college Katholische Stiftungshochschule München                     https://www.ksh-muenchen.de

Well I wonder if there must be missing something in the request.


texts: Measuring Distances between Points with the same timestamps using QGIS

texts:
I'm facing the problem that I want to calculate flight distances of Red Deer.
I have a Layer with point Locations of animals; One attribute is the timestamp. I have another layer with point locations for humans who were trying to disturb the deer.

The problem is, I can't just calculate the distance between the points, because if the animal and the human are moving, the results will be wrong: if the animal moves away from the human, their paths will appear like they crossed although they both were not there at the same time.

So long story short:

I want to calculate the distances between human and animal simultaneously. What I want QGIS to do is to select all locations with the same/matching timestamp and calculate their distance.


texts: Point in Polygon: Aggregation does not work

texts:
I'm working with these data: click here.
The task is actually pretty simple but somehow it won't work properly.

I have a point and a polygon layer. In most polygons there are a couple of points. The point layer has an attribute called EWZ. I just want to add EWZ to each polygon and it shall be the sum of EWZ of all points which lie within the polygon.

I used the point to polygone function where I chose under statistical methods or aggregation method (not sure how it is called in the english version) sum.

But if I generate the now layer, the values are not correctly aggregated. I also tried max and standard deviation. However, in the new field EWZ_Sum always appears just the value of one of the previous points and not the aggregated value.

Is this a bug or am I missing something? I'm using QGIS 2.12.1.


texts: Create a subset of a dataset with a specific grid size in ArcGIS using Data Management Tool

texts:
I am wondering how you can create a subset of a dataset that is a specific grid size (in pixels) using the Clip tool within the Data Management Tool box for raster data.

I was trying to find a hint where to set the grid size in the tool's Environment Settings, but the result, after running the tool, did not apply the new grid size.

Please does anyone handled this before?


texts: Shapely MultiPolygon Construction will not accept entire set of Polygons

texts:
I am new to Shapely, so I am using a work request to learn a bit more about it.  I need to execute a number of unions to build larger geometries from assessment neighborhoods in DC.  I am curious about why a MultiPolygon does not seem to be able to be constructed from a list unless it's sliced (or why what I am doing is incorrect).  More importantly, I cannot seem to get all of my polygons in there.  Here is the set up.

#Open assessment neighborhood shapefile
asmt=fiona.open(data_dir+'ASNBHDPLY.shp')

#Initiate iterator
nbhd=asmt.next()

#Create container for features
nbhd_list=[]

#For each feature...
for nbhd in asmt:
    #...throw the feature in the list
    nbhd_list.append(shape(nbhd['geometry']))

From here, I attempted to generate a MultiPolygon object with both of the following choices.

nbhds=MultiPolygon(nbhd_list)

and

nbhds = MultiPolygon([shape(nbhd['geometry']) for nbhd in asmt])

The second version was tried because of the success of the technique here.

In both cases, an exception was thrown:

NotImplementedError: A polygon does not itself provide the array interface. Its rings do.

I did note, however, that both the following approaches worked:

nbhds=MultiPolygon([nbhd_list[0],nbhd_list[1]])

and

nbhds=MultiPolygon(nbhd_list[:62])

The kicker, however, is that there are 71 assessment neighborhoods and I get the above error as soon as I exceed 62.


texts: raster2pgsql: "Could not allocate memory for INSERT statement"

texts:
I'm very new to raster2pgsql so please bear with me.
I'm trying to load a 60mb .tif into my postgis-enabled database via raster2pgsql using the following command:

raster2pgsql -s 5235 -C -F *.tif public.hrsl_lka | psql -h localhost
-U postgres -p 5432

However, I get the following error message:

ERROR: insert_records: Could not allocate memory for INSERT statement
ERROR: process_rasters: Could not convert raster tiles into INSERT or
COPY statements ERROR: Unable to process rasters

I'm using PostgreSQL 12.4, 64-bit. raster2pgsql has worked fine before for smaller raster files of around 3mb.
Does this have to do with the size of the image?


texts: Algorithm "qgis:difference" returns multiple heterogeneous outputs

texts:
I use the qgis:difference algorithm but it returns multi-type outputs. Sometimes it's a QgsVectorLayer but sometimes a str such as : output_....

My question is : when the algorithm returns an empty result (the input is NULL so the difference too), does QGIS create an empty layer in the output that isn't readable ? Or is it a specific problem of the difference ?


texts: How to display a control outside map container in openlayers 3.5.0?

texts:
I try to associate a control to div with OpenLayers 3.5.0.

I found this tutorial and this subject on GIS SE but the solutions were for OpenLayers 2.
I tried them with OL3 but it doesn't work.

Do you know how can I do that ?


texts: Differences in Area calculation methods Google Earth Engine?

texts:
I calculated the area of my river basin twice and got different results and do not know why/what exactly happens.

At first, I just wanted to have the whole area of my river basin, which I did by:

//calculating area whole basin
var areabasin = ee.Image.pixelArea().divide(1000000).reduceRegion({
  reducer: ee.Reducer.sum(),
  geometry: shape_basin,
  scale: 120,
  maxPixels: 1e13
});
print('Area of the whole basin: ', areabasin);

shape_basin is a shapefile I uploaded as an asset that covers the area of interest. As a result, I got 580,453, which I assume must be in km2 as ee.Image.pixelArea is calculating it in m2 and I divided it by 1.000.000
Later, I classified this area and calculated the area of the individual classes. This I did by:

// calculating the areas by class

var classes = result.select('classband');
var area = ee.Image.pixelArea().divide(1000000);

//water
var waterMask = result.updateMask(classes.eq(1));
var waterarea = waterMask.multiply(area).select('classband').rename('water');

//crops
var cropMask = result.updateMask(classes.eq(2));
var croparea = cropMask.multiply(area).select('classband').rename('cropland');

//urban
var urbanMask = result.updateMask(classes.eq(3));
var urbanarea = urbanMask.multiply(area).select('classband').rename('urban');
Map.addLayer(urbanMask, ['FFFFFF'],'urbanonly');

//bare
var bareMask = result.updateMask(classes.eq(4));
var barearea = bareMask.multiply(area).select('classband').rename('bareland');

//forest
var forestMask = result.updateMask(classes.eq(5));
var forestarea = forestMask.multiply(area).select('classband').rename('forest');

//calculating the area of the waterclasses

var area_image = waterarea.addBands(croparea)
                     .addBands(urbanarea)
                     .addBands(barearea)
                     .addBands(forestarea);

var areas = area_image.reduceRegion({
  reducer:ee.Reducer.sum(),
  geometry: shape_basin,
  scale: 120,
  maxPixels:1e13
});
print(areas);

As a result for the five classes, I get:

bareland: 769703.0088548877
cropland: 421111.16084957926
forest:   505036.9907116392
urban:    163561.49438641773
water:     21118.550734309763

Which results when added up in 1,880,531, or multiplied by 120 (scale), it results in 225,663,745.

Any ideas why the difference is so big? How could I change this, what am I doing wrong?

Here the link to the full script:
https://code.earthengine.google.com/59afd4acda14d6b99effac96bf620015


texts: Using GIS in Civil Engineering construction field?

texts:
I am thinking "How to use GIS when construction HighWays,Buildings,pipe lines etc?"

Is there any way we can utilize GIS for Civil Engineering construction industry to increase the efficiency of current construction work flows?

Any examples of using GIS in Civil Engineering construction field?


texts: I want to develop a GIS application in Java. I found some of useful software like UDIg, Openjump, KOsmo and gvsig, but I am unable to decide from which should I start as I want to create my application which I can give to my clients.


texts:
I have some big problem. I can't open GeoServer after restart the system. It shows HTTP ERROR: 503
Problem accessing /geoserver/web/. Reason:
Service Unavailable

Caused by: org.springframework.beans.factory.BeanCreationException: Error creating bean with name 'DiskQuotaMonitor' defined in URL [jar:file:/D:/GeoServer%202.9.1/webapps/geoserver/WEB-INF/lib/gs-gwc-2.9.1.jar!/geowebcache-diskquota-context.xml]: Invocation of init method failed; nested exception is java.lang.IllegalStateException: Could not locate a layer or layer group with id LayerGroupInfoImpl--39f2944:15948bbab19:-8000 within GeoServer configuration, the GWC configuration seems to be out of synch
    at org.springframework.beans.factory.support.AbstractAutowireCapableBeanFactory.initializeBean(AbstractAutowireCapableBeanFactory.java:1578)
    at org.springframework.beans.factory.support.AbstractAutowireCapableBeanFactory.doCreateBean(AbstractAutowireCapableBeanFactory.java:545)
    at org.springframework.beans.factory.support.AbstractAutowireCapableBeanFactory.createBean(AbstractAutowireCapableBeanFactory.java:482)
    at org.springframework.beans.factory.support.AbstractBeanFactory$1.getObject(AbstractBeanFactory.java:306)
    at org.springframework.beans.factory.support.DefaultSingletonBeanRegistry.getSingleton(DefaultSingletonBeanRegistry.java:230)
    at org.springframework.beans.factory.support.AbstractBeanFactory.doGetBean(AbstractBeanFactory.java:302)
    at org.springframework.beans.factory.support.AbstractBeanFactory.getBean(AbstractBeanFactory.java:197)
    at org.springframework.beans.factory.support.BeanDefinitionValueResolver.resolveReference(BeanDefinitionValueResolver.java:351)
    ... 66 more
Caused by: java.lang.IllegalStateException: Could not locate a layer or layer group with id LayerGroupInfoImpl--39f2944:15948bbab19:-8000 within GeoServer configuration, the GWC configuration seems to be out of synch
    at org.geoserver.gwc.layer.GeoServerTileLayer.getPublishedInfo(GeoServerTileLayer.java:313)
    at org.geoserver.gwc.layer.GeoServerTileLayer.getLayerInfo(GeoServerTileLayer.java:292)
    at org.geoserver.gwc.layer.GeoServerTileLayer.getLatLonBbox(GeoServerTileLayer.java:256)
    at org.geoserver.gwc.layer.GeoServerTileLayer.checkGridSubsets(GeoServerTileLayer.java:783)
    at org.geoserver.gwc.layer.GeoServerTileLayer.getGridSubsets(GeoServerTileLayer.java:771)
    at org.geowebcache.diskquota.LayerCacheInfoBuilder.findOnDiskTileSets(LayerCacheInfoBuilder.java:148)
    at org.geowebcache.diskquota.LayerCacheInfoBuilder.buildCacheInfo(LayerCacheInfoBuilder.java:110)
    at org.geowebcache.diskquota.DiskQuotaMonitor.launchCacheInfoGatheringThreads(DiskQuotaMonitor.java:413)
    at org.geowebcache.diskquota.DiskQuotaMonitor.startUpInternal(DiskQuotaMonitor.java:272)
    at org.geowebcache.diskquota.DiskQuotaMonitor.startUp(DiskQuotaMonitor.java:231)
    at org.geowebcache.diskquota.DiskQuotaMonitor.afterPropertiesSet(DiskQuotaMonitor.java:188)
    at org.springframework.beans.factory.support.AbstractAutowireCapableBeanFactory.invokeInitMethods(AbstractAutowireCapableBeanFactory.java:1637)
    at org.springframework.beans.factory.support.AbstractAutowireCapableBeanFactory.initializeBean(AbstractAutowireCapableBeanFactory.java:1574)
    ... 73 more


texts: How to enable GetFeatureinfo for WMTS Layer in ArcGIS JS v4.18?

texts: I have my server hosted locally with GeoServer and I'm using ArcGIS(v4.18) to load maps on the frontend. I don't find a way that says to enable popUp for WMTS layer.


texts: GeoServer + Geofence to secure WMTS layer

texts:
I am using the GeoServer + Geofence (embedded) plugin and hoping to use this to secure a GWC layer but I can't seem to get it working.
I have added Geofence to the Provider Chain as the first entry.
Example setup:
- Workspace: W1
- Layer: L1
- User: u1
    ROLE: view_w1
- User: u2
    NO ROLE

I am able to secure this setup without geofence for WMS/WFS (as expected based on the docs). When using geofence, I was expecting this to secure the gwc as well based on the following rule, however it doesn't seem to behave any differently to the standard security workflow.
P    ROLE     USER  SERVICE  REQUEST  WORKSPACE   LAYER   ACCESS
0 | view_w1 |   *  |   *   |    *   |     W1    |   *   | ALLOW
1 |    *    |   *  |   *   |    *   |      *    |   *   | DENY

It seems like anything I enter in the Service Security overrides Geofence, which is obviously not desired. With the above setup, I can hide the layers in W1 for WMS/WFS but they are still visible in the WMTS getCapabilities response when requested as user "u2".


texts: How to re-project shapefile to be compatible with Openlayers satellite imagery

texts:
The goal: Visualize a map using google's satellite (and projection) using my collected points.

Right now I have a shapefile of my collected points from surveying. I created it by adding my csv file. The shapefile is working and has a projection of WGS 84/ UTM 48N. My coordinates are in decimal degrees.

When I try on the fly projection to EPSG:3857 following the steps provided in this answered question: OpenLayers Plugin Error: latitude or longitude exceeded limits

I get my points in the same relative position, but can only be seen at a scale of 11:1 or greater. This also makes it so the points show up in the wrong location. (13N 112E shows up on Sumatra--the island) Why is this/how can I fix it? The scale should be more like 1:1000.

Running Mac OSX Yosemite, QGIS 2.6.1 (have tried 2.4.x and 2.6.0), OpenLayers 1.3.6


texts: Feature geometries sometimes treated as equal, sometimes not (PyQGIS)

texts:
QGIS 3.4x64 and Python 3.7 seem to sometimes treat my geometries as equal, sometimes not.

I'm debugging code using a reduced dataset with 3 points. Points A and B have the same geometry, and point C is different. I verified this using geometry.isEqual().

I need a dictionary with the format { geometry : [point_ids] }, where the key is the geometry of a point, and the value is a list of all points with that geometry.

When I run the following code, my dictionary should be { geomA : [A, B], geomC : [C] }.

reverse_geoms.setdefault(g, list()).append(v)

Instead, I get { gA : [A], gB : [B], gC : [C] }.

So this means that I have a duplicate key, because gA == gB. Or, it means the geometries aren't actually equal, but I checked and I believe they are. Is there something wrong with my code, or am I misunderstanding something about geometries?


texts: How to select all multipoint features from a point layer?

texts: I'm connecting to WFS layer by qgis. On some layer I've got points features and few multi-point features. Whole layer contains about 1000 features, how can I do some query just to know which features is multipoint geometry?


texts:
At the Moment I try to load a WMS-Map into a Java Geotools Project.
I followed every Step in the Userguide of Geotools.
(here the Link: http://docs.geotools.org/latest/userguide/extension/wms/wms.html#getmaprequest)

At the end of the Guide it shows me how to Convert my GetMapResponse to a BufferedImage.

Until now I did not found a method, to get my Image into the MapContet.
The ImageLab Tutorial from Geotools imports an Image but as a File and not a BufferedImage.
http://docs.geotools.org/latest/userguide/tutorial/raster/image.html

So Here the Code i created from the userguide

public void loadWMS() {
    URL url = null;
    BufferedImage image = null;
    try {
        url = new URL(
                "HIDDENURLHERE");
    } catch (MalformedURLException e) {
        view.setStatus("WMS-Url Fehlerhaft");
    }

    WebMapServer wms = null;
    try {
        wms = new WebMapServer(url);
    } catch (IOException e) {
        view.setStatus(
                "Fehler beim Verbinden Server evtl. Server offline");
    } catch (ServiceException e) {
        view.setStatus("Server sendet eine ServiceException");
    } catch (SAXException e) {
        view.setStatus(
                "Parsen der Server-Antwort nicht möglich oder nich valie Capabilities");
    }
    WMSCapabilities capabilities = wms.getCapabilities();

    GetMapRequest request = wms.createGetMapRequest();

    request.setFormat("image/png");
    request.setDimensions("600", "600");
    request.setTransparent(true);
    request.setSRS("EPSG:4326");
    request.setBBox("-180, -85.0511287798, 180, 85.011287798");
    for (Layer layer : WMSUtils.getNamedLayers(capabilities)) {
        request.addLayer(layer);
    }
    try {
        GetMapResponse response = wms.issueRequest(request);
        image = ImageIO.read(response.getInputStream());
    } catch (ServiceException e) {
        view.setStatus("Server sendet eine ServiceException");
    } catch (IOException e) {
        view.setStatus(
                "Fehler beim Verbinden Server evtl. Server offline");
    }


texts: Does GEE sampleRegions() include pixels along polygon perimeter?

texts:
Does anyone know if GEE sampleRegions() includes pixels that are on the edge of the polygon? I want to train a classifier with polygons, and want the classifier only to look at pixels that are completely inside by the polygon, ie not pixels that are touching the polygon perimeter.
For example, in this image, are pixels that are bisected by the polygon perimeter included in the call to sampleRegions?

Here is an excerpt from the GEE guide, second code block on this page, that shows using the sampleRegions function in context of a simple supervised classification.
// Get the values for all pixels in each polygon in the training.
var training = image.sampleRegions({
    // Get the sample from the polygons FeatureCollection.
    collection: polygons,
    // Keep this list of properties from the polygons.
    properties: ['class'],
    // Set the scale to get Landsat pixels in the polygons.
    scale: 30
});

I didn't see anything in the API docs about this.


texts: Reading/writing utf-8 encoded data using Fiona?

texts:
So I have a 'cpg' file with my 'in_file' shapefile data (definitely
set to 'utf-8'). I was using 'meta' to create the same schema/crs as
'in_file'...though it runs, my text fields are borked on anything that
isn't ascii (note the cjk - 'ZH_CN' - is coming back as codepoints).

ex:

    {'geometry': {'type': 'MultiPoint', 'coordinates':
    [(13531245.475704141, 2886003.2689278126)]}, 'type': 'Feature', 'id':
    '0', 'properties': OrderedDict([(u'EN_US', u'Taipei City'), (u'ZH_CN',
    u'\u53f0\u5317\u5e02'), (u'ID', 1668338.0), (u'DE_DE', u'Taipei
    City'), (u'ES_ES', u'Taipei City'), (u'FR_FR', u'Taipei City'),
    (u'JA_JP', u'Taipei City'), (u'KO_KR', u'Taipei City'), (u'PT_BR',
    u'Taipei City')])}

My first thought was to change:

with fiona.open(in_file, 'r') as input:

to:

with fiona.open(in_file, 'r', encoding='utf-8') as input:

...but that resulted in:

UnicodeEncodeError: 'latin-1' codec can't encode characters in
position 0-2: ordinal not in range(256)

Reading it without the encoding seems to be the start of my
problem...however I was surprised that the error message arose when I
did the above.

def process_file(self, in_file, out_file, compare_file, dist):
    # open input file and compare file
    # loop over each
    with fiona.open(in_file, 'r') as input:
        meta = input.meta
        # The outFile has the same crs, schema as inFile
        with fiona.open(out_file, 'w', **meta) as output:
            with fiona.open(compare_file, 'r') as compare:
            # Read shapely geometries from file
            # Loop through all shapely objects
            # type(item) = 'dict'
                for item in input:
                    geom = item['geometry']
                    my_shape = shape(geom)

                    # check if it is a multipoint or point
                    if geom['type'] == 'MultiPoint':
                        # break Multipoints into points
                        # i.e. <class 'shapely.geometry.multipoint.MultiPoint'>

                        for pt in my_shape:
                            single_point = self.move_point(pt, compare, dist)
                            mpoint = MultiPoint([(single_point.x, single_point.y)])
                            mpoint_for_merge = shape(mapping(mpoint))
                            # write to outfile - WHERE EVERYTHING NOW WORKS
                            output.write({'geometry':mapping(mpoint_for_merge), 'properties': item['properties']})

                    elif geom['type'] == 'Point':
                        # return of move_point is a shapely geom
                        # i.e. <class 'shapely.geometry.point.Point'>
                        my_shape = self.move_point(my_shape, compare, dist)
                        # write to outfile
                        output.write({'geometry':mapping(my_shape), 'properties': item['properties']})

                    else:
                        raise ValueError('unhandled geometry type: ' + repr(geom.type))


texts: How do I calculate a field to be equal to another?

texts:
How can I calculate the StreetNames field to be equal to the values in the TEMP field in the below feature class attribute table?

I am working in a file geodatabase in ArcMap version 10.1.


texts: Problem with Converting XY Coordinates to LatLon

texts:
Here is the function that I use to convert from pixel x,y to lat,lon coordinate:

    public static Coordinate imageToWorld(Dataset gds, int x, int y)
    {
        double[] adfGeoTransform = new double[6];
        double[] p = new double[3];
        gds.GetGeoTransform(adfGeoTransform);
        p[0] = adfGeoTransform[0] + adfGeoTransform[1] * x + adfGeoTransform[2] * y;
        p[1] = adfGeoTransform[3] + adfGeoTransform[4] * x + adfGeoTransform[5] * y;

        OSGeo.OSR.SpatialReference src = new OSGeo.OSR.SpatialReference("");
        string s = gds.GetProjectionRef();
        src.ImportFromWkt(ref s);
        //src.SetUTM(41, 1);
        OSGeo.OSR.SpatialReference wgs84 = new OSGeo.OSR.SpatialReference("");
        wgs84.SetWellKnownGeogCS("WGS84");
        OSGeo.OSR.CoordinateTransformation ct = new OSGeo.OSR.CoordinateTransformation(src, wgs84);
        ct.TransformPoint(p);
        return new Coordinate(p[0], p[1]);
    }

but this function doesn't work for utm projection and gives incorrect longitude.
for example for an image with the following WKT it gives 35.25,3.65 but it should be 35.29,61.17 @(x=0,y=0). if I uncomment line 12 it gives correct result only for this image.

  PROJCS["WGS 84 / UTM zone 41N",GEOGCS["WGS
  84",DATUM["WGS_1984",SPHEROID["WGS
  84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",0],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","32641"]]


texts: Python pyproj converting from WGS84 to Israel Grid IG 05/12

texts:
Now I see this is duplicate and is answered in
Pyproj projection 0.5 meter differnce from official results (IGD05/12,IG05/12) Israel Grid
I am trying to translate coords from WGS84 (file collected using a GPS) to ITM in Israel.
I am testing the code by comparing results with the station points in the official documentation to see if working correctly
and there are small differences.
So what is missing? Where is the mistake?
ALON
RESULT:   624087.942N 162790.452E
Expected: 624088.355N 162789.564E

My code:
from pyproj import Transformer
stril2012 = '+proj=tmerc +lat_0=31.7343936111111 +lon_0=35.2045169444445 +k=1.0000067 +x_0=219529.584 +y_0=626907.39 +ellps=GRS80 +towgs84=-24.002400,-17.103200,-17.844400,-0.33009,-1.852690,1.669690,5.424800 +units=m +no_defs'
stril = stril2012
transformer = Transformer.from_crs("epsg:4326", stril)
#ALON STATION
lat =  (28.569400 / 60 + 42) / 60 + 31
long = (23.921960 / 60 + 36) / 60 + 34
alt = 55.671
res = transformer.transform(lat, long, alt)
print("ALON")
print('RESULT:   {:010.3f}N'.format(res[1]), '{:010.3f}E'.format(res[0]))
print("Expected: 624088.355N 162789.564E")


texts: Correlating 3 data sets and extracting information from each in ArcGIS for Desktop?

texts:
I am currently trying to find a way to correlate 3 separate datasets as point shapefiles.

I would like to import all three shapefiles, and apply a 5 metre buffer to each point.

Then, for each point, I would like to select any other points within that 5 metre buffer from and add their respective point ID to the attribute table.

e.g. If point 1 from dataset A is within 5 metres of point 2 from dataset B, and point 3 from dataset C then the resultant attitubte table would be:

ID Dataset A     ID Dataset B     ID Dataset C

1                2                3

2                -                -

3                -                -

Even a long winded series of operators or scripting would work.


texts:
I'm creating maps with a specified scale for the map window in QGIS Composer. Project CRS is WGS84, layers are WGS84, onthefly-transformation enabled. The area is approx. 1° x 1° (100 x 100 km).
But when I measure and calculate the scale manually on the printed map and with the measure tool (metres) in QGIS, the scale is significantly different.
The print paper size is the same as set in Composition properties.

Is this a bug or am I doing something wrong??

I'm using QGIS 2.2.0 Valmiera 64bit Win


texts: Using GRASS plugin in QGIS to find the nearest elements in vector map

texts:
I am using GRASS processing toolbox in QGIS to find v.distance: nearest element in in vector map 'to' for element in vector map 'from.'

I am working with two shapefiles, one is a fire layer and the other is a map with admin. The fire layers are points while admin map.I wish to select fires nearest to the admin. I selected the right vector file of the fire points through the drop down  menu 'from' vector map. I also selected the second shapefile of the map admin.

I described the relation between the two nearest feature with dist.
I left the column 'from'empty and described 'to' with line.
But after hitting the run button, the command did not respond, but rather a statement that read 'Missing parameter value: 'from' type (point, line, boundary, centroid area).

What I have not done right?


texts:
I have some GeoNode instances (versions 2.4 & 2.6) running on Ubuntu virtual machines and am trying to restrict a layer to only registered users. According to the docs one should be able to edit permissions of a layer on the GeoNode GUI:

  You can now choose whether you want your layer to be viewed and
  downloaded by

  anyone
  any registered user
  a certain user (or group)

But when I go to my.geonode/layers/my_layer -> "Change Layer Permissions" (logged in as admin = superuser), all I see is this:

The problem is: I don't want to set permissions on specific users or groups but want to let ALL registered users (there are a lot, so I want to avoid adding each single user or creating a group containing all users) view this layer, i.e. deny access to anonymous users. How can I achieve this?


texts: Convert slope file to calculate the percentage of the coverage

texts: I calculated the Slope of the area that I'm working on, and I was wondering if there is a way to calculate the actual percentage of each bracket, meaning the percentage of contribution of each class. The Slope is a floating point type of file and I cant convert it to polygon. Any advice could be helpful.


texts: Polygons overlapping points (not points within the polygon)

texts:
I use QGIS 2.14. I have one layer for the points and one layer with many polygons.
I'm having trouble finding the polygons that overlapped with the points.
Is there a way to do this?
I'm not asking for how many points are within the polygon, but the vice versa.
How many polygons overlapped a certain point. Hopefully, I could have an attribute table with the points that has the overlapping corresponding polygon/s.

Here's a screenshot of what I'm working on.


texts: Purpose of "No action selected" toolbar button?

texts:
At 3.8.3, there is a button on the "Attributes Toolbar" that is always greyed out.  When I hover the cursor over it, the response is "No action selected":

The QGIS documentation makes no mention of this button.  What is its purpose?


texts: I've got a parcels layer from the state and i modified it by adding fields relating to my business as in project number and a reference plan field. The state updated their shapefile by adding new subdivisions. They're all related by oid_1 field which is a parcel id field. some of the polygons have changed and I want to make sure the polygons change and my data fields merge correctly.


texts: looking for documentation about exploiiting OSM file with c++ API

texts:
My application is written in QT4.6, i use the C++ API of QGIS 2.8.2
i had to read and display OpenStreetMap in my canvas in my application
so, i am looking for some documentation or tutorial for reading, exploring, exploiting or coding the OSM file
i must not use the QGIS application, i have a specific application for a project, and the project had to read and display OSM File


texts: I'm working on a project at work to download imagery from else where and add it to geoserver. I have no experience in GIS, nonetheless geoserver. I've been hunting through the REST api for 2.0.X of Geoserver and can create workspaces and empty coverage stores. I get that coverage stores hold raster graphics, but is it one-to-one or do they share a special relation by being in the same store? What's the URL and method (PUT, POST, etc) to add a tiff to a coverage store?


texts: How to add WMTS layer from GeoServer to a map without downloading the get capabilties xml file (in OpenLayers)?

texts:
I am new to web mapping and web in general.

I'm adding the default WMTS Italy mosaic that can be found in GeoServer (nurc:mosaic) to my map by following these examples: http://openlayers.org/en/latest/examples/wmts-capabilities.html
https://openlayers.org/en/latest/examples/wmts-layer-from-capabilities.html

I downloaded wmts-getcapabilities xml file from my GeoServer's web admin page, simply by clicking on WMTS's 1.0.0. Then saved it within the same directory with my Italy_mosaic_OpenLayers.html file.

And here is how my Italy_mosaic_OpenLayers.html looks like :

var parser = new ol.format.WMTSCapabilities();

var map;

fetch('wmts-getcapabilities.xml').then(function(response) {
    return response.text();
}).then(function(text) {
    var result = parser.read(text);

    var options = ol.source.WMTS.optionsFromCapabilities(result, {
        layer: 'nurc:mosaic',
        matrixSet: 'EPSG:4326',
        format: 'image/png'
    });

    map = new ol.Map({
        layers: [
            new ol.layer.Tile({
                source: new ol.source.OSM(),
                opacity: 1
            }),
            new ol.layer.Tile({
                opacity: 0.6,
                source: new ol.source.WMTS(/** @type {!olx.source.WMTSOptions} */ (options))
            })
        ],
        target: 'map',
        view: new ol.View({
            center: ol.proj.fromLonLat([22, 42]),
            zoom: 5,
            projection:'EPSG:900913'
        })
    });
})

This code works only when my "wmts-getcapabilities.xml" file is in the same directory with my "Italy_mosaic_OpenLayers.html" file. However if I were to renew my Italy_ mosaic raster file I would have to download "wmts-getcapabilities.xml" file again and place it to the same directory(by the "Italy_mosaic_OpenLayers.html" file) again.

To avoid replacement of the xml file again and again, I tried to put 'http://localhost:8081/geoserver/gwc/service/wmts?REQUEST=GetCapabilities' instead of 'wmts-getcapabilities-xml' in my code, however I had a blank page.
I also tried different matrixSet-projection combinations; no matter what matrixSet-projection combination I tried (I tried 9 possibilities, 3 for each (4326,3857 and 900913)), with fetch('http://localhost:8081/geoserver/gwc/service/wmts?REQUEST=GetCapabilities') I always had a blank page (not gray).

Here is my fiddle https://jsfiddle.net/qjcadvzq/15/

What am I missing ?


texts:
Is it possible to use WMTS e.g. by http://data.wien.gv.at in QGIS?

I would like to import the data published by the City of Vienna
as a layer in QGIS. Is this possible at all, and if yes, how do I do it?
ich bin kein profi und bitte daher um
unterstützung.

This is what I am talking about:

  http://data.wien.gv.at/formate/geowebservices.html

I could successfully load the WMS service, but I could not do the same with WMTS.


texts: How to create a branch network from several input points to an unique outlet in minimizing the distances?

texts:
I am searching a way to "hydrologically" connect upstream points to a unique downstream outlet by creating an optimized stream network as pictured below...

The cummulated lenght of the stream lines should be minimized.

Otherwise the stream network should remains confined within a defined polygon...

I only know the input points (red) and the output (green) and the polygon. There is no other constraints

Does anyone know how to do it ?...

I searched on the web all the morning but without succeed

No I rightly try to do it without considering the topography... Just a flat (or approx flat) surface.
That's why I am searching a way for so long time :-)

It is an hyopthetical area.

Imagine that you are in you bath.. the water level rises up and overflow on a downstream threshold - at a corner of you bath..

Then you add a few particles at different corners of your bath (upstream) and you try now to sketch what will be their circulation to reach the threshold by forming an optimized network.. The stream network (as I draw on the right) should be minimum regarding the particles input.

Am I understood ?
:-)


texts: How to save results of postgresql to csv/excel file using psycopg2?


texts: OS Translator plugin replaces the unique feature references (TOIDs) with ogc-fid which is an integer

texts:
I've successfully loaded OS MasterMap Topography GZs files into my PostgreSQL (9.5) database, using OS Translator II plugin within QGIS (2.18.18).

But somewhere during the translation process the unique feature references (TOIDs) gets replaced with the ogc_fid ([PK] integer) which is an integer instead of the TOIDs in the GZs files. Here is a screenshot of data in PostgreSQL table:

Does anybody know how I can get the TOID into the PostGIS?


texts: Problems filtering MOD09A1 in GEE using the QA layer

texts:
I'm trying to mask each modis image using the QA layer. As far as I understand, the part  .select("QA").eq(0) should remove all pixels with no ideal quality for all bands. However, when I make a plot with Map.addLayer(ls7.first()) it seems like there is no pixel in any part of the world with the best quality.
It seems strange, is something wrong with my code? When I use .select("StateQA").eq(0) it seems that only some sparse pixels (most of them in the sea) are ok. Any suggestions?
This is my code:
//select collection

var ls7SR = ee.ImageCollection('MODIS/006/MOD09A1');

// calculate normalized burn ration

var ls4_7_Indices = function(lsImage){
  var nbr = lsImage.normalizedDifference(['sur_refl_b02', 'sur_refl_b07']).toFloat();
  nbr = nbr.updateMask(lsImage.select("QA").eq(0));
  return nbr.copyProperties(lsImage, ['system:time_start']);
  };

// Map functions across modis Collection

var ls7 = ls7SR.map(ls4_7_Indices)

Map.addLayer(ls7.first())


texts: Exporting vector from raster classification gives deformed vector in Google Earth Engine?

texts:
I am very new with Google Earth Engine and I try to make a landcover classification using Sentinel-2 imagery to get the urban area over a very large area (country wide ideally).

I was able to realize unsupervised classification and convert it into vector ("small" area for test) and the result displays well. But when I want to export on my drive and use it in ArcMap, the results is very different from what I can see in GEE.

Below is the code I used:

var AOI = table
function maskS2clouds(image) {
  var qa = image.select('QA60');

  // Bits 10 and 11 are clouds and cirrus, respectively.
  var cloudBitMask = 1 << 10;
  var cirrusBitMask = 1 << 11;

  // Both flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudBitMask).eq(0)
      .and(qa.bitwiseAnd(cirrusBitMask).eq(0));

  return image.updateMask(mask).divide(10000);
}

// Map the function over one year of data and take the median.
// Load Sentinel-2 TOA reflectance data.
var dataset = ee.ImageCollection('COPERNICUS/S2')
                  .filterDate('2018-12-01', '2018-12-30')
                  .filterBounds(AOI)
                  // Pre-filter to get less cloudy granules.
                  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20))
                  .map(maskS2clouds);

var median = dataset.median();

// Set the training AOI
var training = median.sample({
  region: AOI,
  scale: 10,
  numPixels: 5000
});

// Train the cpu to make classification
var clusterer = ee.Clusterer.wekaKMeans(15).train(training);
var result = median.cluster(clusterer);

// Only select Nola area.
var Nola = result.clipToCollection(AOI);

// Create a mask to only keep urban area.
var maskUrban = result.updateMask(result.eq(10));

// Conversion To vector.
var clusterValue = maskUrban.select('cluster');

var classes = clusterValue.reduceToVectors({
  reducer: ee.Reducer.countEvery(),
  geometry: AOI,
  scale: 10,
  maxPixels: 1e8
});
var result2 = ee.FeatureCollection(classes);
Map.addLayer(result2, {},'Urban vector');

// Export Vector data.
Export.table.toDrive({
  collection: result2,
  description: 'CAR_2018_VectorUrban_Nola_2',
  fileFormat: 'SHP'
});

Map.addLayer(median, {bands: ['B8','B4','B3'], min: 0.0, max: 0.3}, 'S2 - 842');
Map.addLayer(result.randomVisualizer(), {}, 'clusters');
Map.addLayer(maskUrban, {},'Mask urban');
Map.centerObject(AOI,9);

The result I got in GEE:

The exported result as seen in ArcMap:

If anyone of you has an idea of what is going on and/or how it is possible to fix it, let me know.


texts:
On a fairly regular basis at my place of work we get inquiries of the form "where is Lake Potato?" or "where is Poppingstone Creek?" or... The questions come from co-workers, the general public and other organizations. We're a government agency and sell maps, so a natural target for that kind of query.

My method for answering generally resembles:

search the official names db (in Canada we have the Canadian Geographical Names DB),
Scan our local data stores and maps
A general internet search
Send an email 'round to field workers (biologists, conservation and parks officers, etc.)
...

...but what to do if you get all the way to #5 and still need an answer? I could post a story on my website, and one of the 7 people who read it this year might know. GIS Stack Exchange is peopled by the kind of folk who might know, but it's off topic. Indeed, even this question about the question is close to being off topic (read as: it's ok with me if the general consensus is to move this to meta).

So my question is, is there any place where "where is such-n-such?" is on topic?


texts:
We have a number of mapservices in ArcGIS Server 10.0 that reads data from FGDBs.
We are now moving that data in to ArcSDE and the plan is to use this data from AGS.

Is it better to connect to SDE using Direct Connect or going through the ArcSDE services performance wise?


texts: Display 8-Bit Signed Rasters with Negative Values in QGIS / MapServer

texts: When loading 8-bit signed rasters with values from -128 to 127 into QGIS or MapServer, negative numbers are not recognized. Is there a way to change the default display settings in QGIS and/or configure the map file for MapServer to recognize the negative values? For MapServer, if classes were setup for the display properties would that work (may just try this out). FYI - These rasters display correctly in ArcGIS Desktop.


texts:
I tried to define area of NDWI from 2017-2019 but only show list of date and Area (F can not change).

Can you fix my code for me?

https://code.earthengine.google.com/f0955df742de0f0b533f7778fd6c9349

or

// Data
var Sentinel2 = ee.ImageCollection("COPERNICUS/S2");

var Lake = ee.FeatureCollection("ft:1DZkHTe3rWHS7DDfBcxMQ-DUTtkrDgcYmw1TVO6Gt")
Map.centerObject(Lake,15);
Map.addLayer(Lake,{},"Ha Thuong")
// filter data and location
var collection = Sentinel2.filterDate
('2017-12-12', '2019-9-5').filterBounds(Lake);
print(collection);
print(collection.toList(100).get(0));

var a = collection.toList(collection.size()).map(function(img){
img = ee.Image(img)
// Display true color image
//var Sentinel2_ = Map.addLayer(collection.min().clip(Lake), {bands: ['B4', 'B3', 'B2'], max: 2048}, 'Sentinel2_');

// calculate ndwi from sentinel
function s2ndwi(img){
  var ndwi = img.normalizedDifference(['B3', 'B8']).rename('NDWI');
  return img.addBands(ndwi);
}
// run ndwi on all images
collection = collection.map(s2ndwi);

// select maximum NDWI
s2ndwi = collection.select("NDWI").max();

// set the threshold

var THRESHOLD = 0.1;

// select pixels greater than threshold
s2ndwi = s2ndwi.gt(THRESHOLD);

// set visualization
//var ndwi_viz = {bands:"NDWI", min:-0.3, max:0.3, palette:"000000, 0000ff"};
// add the map as a layer
//var ndwi_2017 = Map.addLayer(s2ndwi.updateMask(s2ndwi).clip(Lake), ndwi_viz, 'ndwi_2017' );

var area = ee.Image.pixelArea().divide(10000)
var area_image = s2ndwi.multiply(area).select([0],['F ha']);

var areas = area_image.reduceRegion({
  reducer:ee.Reducer.sum(),
  geometry: Lake, //.geometry(),
  scale: 10,
  maxPixels:5e13,
  // return: collection
})

var idate = ee.Image(img).date();

return [idate, areas]
})
var l = a.map(function(b){
  var c = ee.Feature(b)
  return c
})
print (l)


texts: Identifying differences in polygons across shape files

texts:
I am trying to compare the 2000 Census county boundaries with the 2010 Census county boundaries. Each county is uniquely identified in a given year by a FIPS code and I want to compare if the polygon represented by given FIPS code z in 2000 is different from the polygon represented by the same code (z) in 2010. By different I mean any change to the geometry, regardless of whether it modifies the area of the polygon.

My ideal output is a table/spreadsheet with two columns - one containing the 2000 Census FIPs codes and the other indicating some way of identifying if a that code references different geometries/polygons in the two years. This indicator could be a boole (TRUE/FALSE), a binary indicator (1/0), some continuous variable, or any other indicator.

I am looking for a solution in either R or QGIS.

Notes:

This post is similar but isn't concerned with geometry while I am concerned with geometry: Polygon overlay, calculating changes, in QGIS?
It looks like this would provide the solution: http://postgis.net/docs/ST_Equals.html, but I am preferably looking for a solution in QGIS or R.
Something somewhat related is this example in ArcGIS: Finding differences between shapefiles using ArcGIS for Desktop?
Also similar but using ArcGIS: Finding differences between shapefiles using ArcGIS for Desktop?


texts: Link to google maps ignores url parameters for streetview

texts:
I have a web site where I need to generate link to Google Maps in Street View mode. I have available the coordinates to generate the URL.

When I try this URL in a desktop browser, it does what I want, and takes me to a street in Norway, and activating street view:
https://maps.google.com/maps?spn=0.008443,0.022767&ctz=-60&t=h&z=16&layer=c&cbll=58.819625,5.75462&cbp=12,44.85,,0,18.41

When I try exaclty the same URL in a mobile browser, it just ignores the parameters, and redirects to http://maps.google.com .
When I try another link that is not streetview, it works.

I noticed it first in an iPad, and then tried in a Samsung Galaxy using Firefox, and same there.

I even tried using Firefox on desktop using the mobile device emulator plugin emulating an iPad, and same there. So it is Google that does this if the user agent reports itself as a mobile device.

Then I tried to navigate to a specific location in Google Maps on the mobile device, and clicked on the street view icon. Then a new window was opened with this URL:
http://maps.gstatic.com/m/streetview/?panoid=vZ_OD_yZLEJYC5wQQ9q5zg&cbp=0,247.65935203825757,,0,0

This works for this location, but I need to generate the URL based on coordinates, so I don't have a valid panoid to use. Is panoid the only way to use url parameters for displaying street view on mobile device?
If so, is there any API available for generating a panoid based on coordinates?


texts: Are there limits on the number of unique polygons created by the QGIS Union tool?

texts:
Is there a limit on the number of polygons that the Union tool will compute? I have overlapping radii from a variety of sources pictured below, but the tool continues to give me an "GEOS geoprocessing error: intersection failed" error every time I try with a larger dataset.
I've tried to fix geometries and checked validities for topology, are my overlaps just too numerous so the tool fails? Is there anything I can try, other than splitting my data into smaller and smaller peices and Unioning them individually until I get my final solution?


texts: How to keep (complex) color table when exporting from GRASS

texts:
I understand how to work with color tables within GRASS - but is there any way to preserve color tables when opening the raster in a different program?

For example, I created a color table in GRASS that is very detailed (certain number ranges get certain colors), but then when I open the file in QGIS, all that work is gone. It does not seem to matter what type of raster file I export as from GRASS (I have tried geotif, ArcGIS Binary, and ArcGIS ASCII). Is there any trick I am missing?

Thanks if you can help.


texts: Getting layer in WMS and editing in WFS-T with OpenLayers

texts:
I have huge datasets on PostGIS. For example building table includes more than 100k features in POLYGON geometry. When I want to get all data in QGIS (2.1) or other GIS software by using direct DB connection I wait about 30 sec. But if use WFS connection (GeoServer) it takes more time (don't try)... In a web GIS project (using OpenLayers) people don't want to wait in front of the screen. And I think OpenLayers can't handle with that amount of data... So I use to serve data by using WMS. It's really too fast and maybe later I try tilecaching etc... OK everything is so far so good. But I want to edit features (only insert maybe delete, but no update) in the browser.

I have checked OpenLayers examples for WFS-T. These examples load all data to by WFS. I don't want loading data just inserting new features by drawing and than clicking save button. So new features sending to GeoServer...

I tried with OpenLayers WFS-transaction example link. Firstly modified data map.addlayers([wfs,gpyh]); to map.addlayers([gpyh]); to load only googlemaps baselayer. There is no error during drawing but when saving I get error. I says

this.layer.map is null
var local = this.layer.map.getProjectionObject();

What do you advise me when working with big layers.Is my way correct. How can I handle with this issue? My JavaScript knowledge is limited.

EDIT

Some info about the system which I want to build:

I want to serve buildings, roads, boundary etc. layers to people. Not only showing but also editing capabilities it should have.

MY PROBLEM:

I don't want to load features to OpenLayers, I just want to draw new features and put them into database. And I use WMS for serving layers.

I use OpenLayers lib in client side (apprentice in JavaScript).

The data is stored in PostGIS., but I can use Spatialite.

GeoServer for WMS and WFS service.

I am not bad in Python and PHP. Maybe I can use FeatureServer. Also I can serve feature data with PHP or Python.


texts: Reading files from grass into R with loop

texts:
I have a connection between GrassGIS and R. I want to read thousands files from my mapset into R. I have lines with names "line1","line2"...
I know I have to use a loop for this so I did

list_lines=list()
count=1
for (i in 1:3000){
    list_lines[count]=readRAST(c(line[i]),cat=c(FALSE))
    count=count+1
}

But it doesn't work. What should I do?


texts: Determining GPS Position from Elevation, Azimuth, SNR

texts:
I'd like to calculate the receiver GPS location myself, without relying on high-level API. A few questions there: I am looking at an Android app that can show GPS satellites, on a diagram like this

However Android (until API 25) does not return satellite pseudorange, which every source I read is required to compute a location yourself. Android API does expose signal-to-noise ratio, elevation, azimuth, and PNR (pseudorandom number) for each satellite in view. Now I believe the elevation and azimuth are given in relation to the receiver (hence the image above), so I was wondering, given this information and precise location of each satellite (the PRN uniquely identifies a satellite AFAIK, I could look that up from an almanac, or ephemerisis file), couldn't we compute a rough location?

This gentleman here can compute the location of a GPS satellite given receivers location, satellite elevation and azimuth. I want to merely do the reverse, given satellite location (from almanac), azimuth and elevation, compute receiver's location.

Another link on azimuth / elevation and Android.

And this person seems to be talking about the distance computation, and even computing the receiver location..?


texts: How to merge adjacent polygons with no common attributes in GRASS?

texts:
Many topics are talking about merging polygons based on their common attributes, but few (or none) give a real answer to location-based merging problem.

Explanation.

Here is a fraction of the base map:

Some polygons are adjacents but belong to different entities. I want them to merge when they touch and stay as they are if they don't.
Then, I want to extract the area of the new entities out of it.

Due to the large size of the maps, I'm working on Grass.

To solve the problem, I first made a tiny buffer:

grass_util.g('v.buffer', input='map', output='map_buffer', type=area distance=0.001)

My output is this one:

Everything looks pretty good, except that making a buffer brings the entities to merge into the one same entity (as displayed by the attr table that I can display here due to my small reputation).
More precisely, the number of entities is almost the same, a bit less, the ones that didn't merge plus the new merged ones, but they are all affected "1" in the only category: "cat".

So we still have our different polygons, different entities, but they are merged.

Following this, when I try to add a new column and fill it with geometrics (my area in ha):

grass_util.g('v.db.addcolumn' map='map_buffer' col="area int")
grass_util.g('v.to.db', map='map_buffer', option=area, column=area, unit=h)

It lamely fails !
The area isn't calculated for each entity and attributed to it. It is calculated for the whole package, and every entity gets the area of the whole.

I have to say that I don't really see how to get this thing sorted.
I tried dissolve, v.clean, v.build, multiple to one and one to multiple, I tried to deconstruct the polygones to lines and to reconstruct them again.

An option which need the processing module that I don't have yet is to use this command:

processing.runalg('qgis:exportaddgeometrycolumns', input, calc_method, output)

which works when I try on a small area in Qgis Gui.

The question is twofold:

Find an other way to merge the adjacent polygons based on their location (to avoid the destructive buffer)
Find a way to get back my entities and retrieve their personal areas.


texts: Calculating slope and aspect from DTM lidar elevation data in ArcGIS Desktop?

texts:
I have an ASCII 2m lidar DTM file which I converted into a raster in ArcGIS with the From ASCII to Raster tool. From this DTM file I would like to calculate the slope and aspect. I have used the spatial analyst -> surface toolbox for this. The DTM has a BNG projection and all values are calculated in degrees. I followed several different tutorials and they all show a same result, a 'equal and smooth' raster with slope or aspect values (ArcGIS and GIScourse).

However, when I try to do this for my DTM raster, I get both 'smooth' parts (left side of pictures) but also what seems like contour lines (right side of pictures). It seems like they follow the lines of the DTM where values change, not really changes in slope etc.

Had anybody had the same problem and/or can anybody tell me what I am doing wrong?


texts:
In QGIS python console, I can run the following to extract the week number from date column. How can I use the same in qgsExpression as a where clause? Any other option is welcome as long as it uses PyQGIS. I intend to filter a vector layer in QGIS to return records that match a particular week.

for x in layers:
if x.name() == "test":
    idx = x.fieldNameIndex("date_from")
    values = x.uniqueValues( idx )
    b =values[0]
    print b.weekNumber()
    b.weekNumber()[0]

This expression gives me an indication of the week number. Can I use QgsExpression like
 expr = QgsExpression( "\"date_from.weekNumber()\" = 99"  )?


texts: How to create a fishnet/grid along a polyline or polygon

texts: I am trying to generate a fishnet/grid along a line/polyline or polygon boundary in ArcGis 10.3 or QGIS. While it is totally intuitive and easy to generate a box-standard grid, what I am trying to do is to basically generate a series of cells along the boudary of the Italian peninsula. One grid should have the center of each cell positioned along the boundary itself, and a second grid shoud have cells along the boundary but tangent to the boundary itself.


texts: QGIS raster calculator using second raster to set value in first

texts:
I have two rasters.  ras1@1 has values from 1-50, ras2@1 has values 0 and 1.  I want to convert values from ras1@1 >= 4 to 1, except where ras2@1 = 1, which I want to set to 0.  I cannot work out how to do this.  I was thinking something like:

(ras@1 >= 4 AND ras2@1 = 0) * 1

or

((ras1@1 >= 4) * 1) - ras2@1

but neither seem to return what I want.  Everything I have tried seems to return just ras2@1 = 1.


texts: Selecting records by SQL Query and saving SUM or COUNT result into differnt table

texts:
I am trying to create a python toolbox with the following function:
I want to select points or lines within a polygon and save the selected feautres as a temporarily Feature. (That works)

Then I want to use SQL-Queries for example to count all Features of the Point-Feature-Class or the total length of a Line-Feature-Class. The result of the queries should be written into an existing Table. For example the result of the total length of a Line-Feature-Class, should be written into the Field: TOTAL_LENGTH.

At the end, the Table should be exported as an excel-file (Should be no problem)

I don't want to get the complete finished code. Just an Idea how to do the query and saving the result into a different existing table.

Does somebody has an idea?


texts: Pansharpening Landsat 8 16 bit raster images

texts:
Situation:
I have merged six Landsat 8 image tiles to one big tile for the bands 4,3 and 2. These are Landsat level 2 products, so already corrected to surface reflectance and need no further processing.
I have separately downloaded six band 8 image tiles for the same region. They are of a different date.merged those to one tile too.
All the raster images are 16 bit.
What options do I have to perform pansharpening on the 432 band combination?
Grass Gis apparently only takes 8 bit.
Semiautomatic classifications plugin in QGIS requires a metadata file. Of course no one single metadata file applies once the tiles are merged to one. Also, the band 8 is of a different date than the other bands.
please advise.


texts:
I need to set up loadMask for OpenLayers.Control.GetFeature. I know I should hide the mask on featureselected event or on protocol callback but when should I show the mask?

Here is code (assume I have _mask global variable with show and hide methods):

var selectionControl = new OpenLayers.Control.GetFeature({
            protocol: this.map.layerUtils.createWFSProtocol({featureType:featureType, srsName:'EPSG:'+ this.map.projectionCode}),
            clickout:true,
            box: true,
            toggleKey: "ctrlKey"
        });

selectionControl.events.register("featureselected", this, function(e) {
            _mask.hide();
        });


texts: Export mosaic to Google Drive, generating more than just one file and can't download this

texts:
I'm exporting a mosaic via google earth engine to google drive, however it has generated many files inside a folder.

// Create a geometry representing an export region.
var roi = ee.FeatureCollection('ft:15Y2IjO7bMc-COe8xOm3x1tfZdgb7GO4hK1TDWxGg');

// Load the Sentinel-2 ImageCollection.
var sentinel = ee.ImageCollection('COPERNICUS/S2').filterBounds(roi);
var s2 = sentinel.filterDate('2017-01-01', '2017-12-31');
var s3 = s2.filterMetadata('CLOUDY_PIXEL_PERCENTAGE', 'less_than', 0.1).sort("ascending");

// Mosaic

var mosaic = s3.median()
var bands = mosaic.select('B2', "B3", "B4", "B8")
Map.addLayer(bands,{}, "mosaic")

// Export
var task = Export.image.toDrive({
  image: bands,
  region: roi.getInfo()["coordinates"],
  scale: 10,
  description:'MOSAIC',
  folder: 'MOSAIC',
  fileNamePrefix: 'MOSAIC',
  maxPixels: 1e13
})

How to export to generate fewer files. How to download large files as mosaics to your computer?

Result from export


texts: Merge a list of spatial polygon objects in R

texts:
I have a list of spatial buffers (30000 buffers) that I built with the function lapply:

buff.pts <- lapply(1:nrow(pts.prj), FUN=function(l){
  buff <- gBuffer(pts.prj[l,], width=1000) ## 1km
  return(buff)
}))

> head(buff.pts)
[[1]]
class       : SpatialPolygons
features    : 1
extent      : 307941.8, 311941.8, 4994518, 4998518  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=tmerc +lat_0=0 +lon_0=-73.5 +k=0.9999 +x_0=304800 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs +towgs84=0,0,0

[[2]]
class       : SpatialPolygons
features    : 1
extent      : 307226, 311226, 4991153, 4995153  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=tmerc +lat_0=0 +lon_0=-73.5 +k=0.9999 +x_0=304800 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs +towgs84=0,0,0

From this list, how can I merge all spatial buffers to obtain a shapefile with the 30000 buffers (or features) ? (This shapefile will then be used in the function aggregate to aggregate spatial polygons by attributes.)

I tested this code but I obtain this error message:

test <- as.data.frame(do.call("rbind", buff.pts))
Error in as.data.frame(do.call("rbind", buff.pts)) :
  error in evaluating the argument 'x' in selecting a method for function 'as.data.frame': Error in validObject(res) :
  invalid class “SpatialPolygons” object: non-unique Polygons ID slot values


texts: Django app can't upload multi-band-image

texts:
I try to upload a TIFF raster image with 5 bands in Django app and I get error message in form : Upload a valid image. The file you uploaded was either not an image or a corrupted image. I am sure that image is valid .

Here the code :

def dataset(request):
    uploadimages = UploadImagesForm(request.POST or None, request.FILES or None)
    if uploadimages.is_valid():
        if request.FILES.get('multipleimages', None) is not None:
            images = request.FILES.getlist('multipleimages')
            for image in images:
                name = str(image)
                .........
    return render(request, 'mylist.html', {'uploadimages': uploadimages})

forms.py

class UploadImagesForm(forms.Form):
    multipleimages = forms.ImageField(label="Upload Multiple Images:",widget=forms.ClearableFileInput(attrs={'multiple': True}), required=False)

Any idea how to fix it ?


texts: Delete custom toolbar grayed out in ArcGIS 10.2

texts: Recently I worked on adding a custom toolbar to ArcMap. It took a few tries, and now there are my initial 'drafts' of the toolbar still added to my ArcMap list of toolbars. According to the help, I should just be able to delete them. However, the Delete button is grayed out when I try to select my custom toolbar, as if it were one of the default toolbars. Has anyone else had this problem and come up with a solution?


texts: Retrieve average Z value of a 3D geometry in postgis

texts:
I need to compute the average Z value of a 3D geometry in PostGIS.
Two functions ST_ZMin and ST_ZMax are available to return the min and the max, but the function ST_ZAverage returning the average value does not exist unfortunately.

Is there a way to compute this average value?

Maybe there is an existing solution based on ST_NPoints, ST_PointN and ST_Z.


texts:
I have two different data sources that I would want to use for my application

  1) Time series data on a road network

  2) Spatial data of the road network

In essence I would want to combine both the information and build visualizations. For example to show variation of speed with time on a particular road using a time slider.

To do this I created a view joining the two tables using sdetable tool and built a layer using the view. But for some reason I am not able to host it as a service.

Please do suggest if there are better ways to do it?


texts: OD Cost Matrix does not work from every origin to every destination?

texts: I'm trying to run OD Cost Matrix function for my term project and i run into a problem. I have 2566 origins and 3 destinations and I am trying to calculate distance base origin cost matrix from every origin point to every destination point but as you can see in figure, only 43 lines appeared. The problem is not about my network. I used the default settings for building the network. The error message is : Warning: No "Destinations" found for "Location 1" in "Origins" for every origin. I could not figure what the problem is.


texts: Is it possible de create a custom line style in QGIS, as mentioned in the attached file?


texts: Cost Path in ArcGIS Desktop drawing no path?

texts:
I'm using the Cost Path tool in ArcGIS Desktop 10.3 to draw thalwegs (line of lowest elevation) on a channel network DEM, but running into some problems. In the most simplified case, the raw DEM is used as the cost distance, and the results of Flow Direction are used as the back link (as suggested by ESRI here), which works. The example below shows the output Cost Path over the Flow Direction, with the destinations drawn as polygons.

Unfortunately, this is a bit over-simplified, mainly because flow direction requires a Fill operation, and so discards pools in the raster and allows the cost path to "cut corners". I created a weighted raster to account for both slope and proximity to deeper areas, and Cost Distance works fine with this and produces a back link raster which looks good. But running Cost Path with the new cost distance and cost back links simply fills the destination polygons, with no connections drawn between them. For example, here is the weighted raster with the output cost path:

As you can see, the destinations themselves are weighted (as usual), but no path is drawn between them. The strange thing is, any Cost Path I run with Flow Direction as the back link always works, but it also gives me the same output path regardless of the cost raster used. But any time I use the back link produced by other tools, it fails to draw a path.

Some things that I have tried that did not solve this:

ensuring geographic projection, processing extent, snap raster, and cell size are identical
reclassing the Back Link to the same scale as Flow Direction (1, 2, 4, 8 .. 128)
running Cost Back Link separately from Cost Direction
using Path Distance instead of Cost Distance or Path Distance Back Link instead of Cost Back Link
trying the process on a smaller raster (the original is already only 1421 x 1619)
trying Cost Path between two destinations only
redrawing the destination polygons so they are fully within the raster bounds

None of these things have worked.

Any ideas?

I'm also aware of a similar function in the River Bathymetry Toolkit Add-In, but that only works in 10.2.

===========================================

Here's some clarification on why using flow direction as the back link fails for my needs. This is based on a high-resolution (1 m) bathymetry DEM from sonar data. The top panel in the image below show the results of simple flow direction in a cost path overlaid on the DEM. Several corners are cut to produce the straightest path by excluding areas over pools. In order for flow direction to work, Fill must first be run, which effectively sets all pools and large swaths of the bottom to a constant elevation. The areas filled by this operation are illustrated in the second panel. The third panel demonstrates the final method I settled on, use of slope as the cost surface, and the line is clearly a better representation of the thalweg here.

Here's the final output:


texts: "Node tool: could not snap to a segment on the current layer" message

texts:
I have a strange issue. Editing a layer with node tool I get always the following message:
"Node tool: could not snap to a segment on the current layer."
This error happens to me even with a fresh install of new 2.6 version only on my xp system (erased .qgis directory in profile and with no plugin at all)
I get this error in all the currently installed release (2.0, 2.2, 2.4, 2.5 weekly, 2.6) and appears to be system dependant as other xp systems near to me are not affected.
What can I do?

(source: nabble.com)


texts: Exporting Zonal Statistics As Table output in ModelBuilder

texts:
I have an extremely simple model that takes as input a shapefile for zones along with a value raster. All the model does is produce a table of my desired zonal statistic for each zone (the mean). Next, I hope to get Modelbuilder to export those results into a txt or other file I can easily manipulate later in Excel.

I have about 48 rasters for which I would like to do this, and I can copy and paste the desired columns side-by-side into the Excel file myself, so although it will be a bit tedious I'm not worried about that. But it would be nice if there were a way to automatically export the table. Is this possible in Modelbuilder?

Other background info, if helpful:
I found the following exchange intriguing (Iterate through multiple Raster files and compute statistics using Zonal Statistics), and had hoped to be able to iterate my rasters and have them automatically append all in one table which I could later export myself in one easy click. I was unsuccessful in getting this to work, perhaps due to my misunderstanding of properly calling the target files, column names, etc. (I got the models described here to work perfectly well on my data... except that I could not successfully automate adding and naming new fields... my result had one field with the final raster/statistic having overwritten the others.)

I am using ArcGIS Desktop 10.2.1 and Spatial Analyst.


texts: How centroid is getting linked to one Parcel if centroid is lying on boundary of adjacent parcel

texts: In my case there Parcel is who is supplier and building centroid is requestor in SpatialRelator of FME Script. You can see black round in one below picture. So even if centroid is falling on adjacent parcel boundary it got link of one parcel. How come ? why it's not gone for one to many matches ?


texts:
Following the directions to install geodjango on linux ubuntu 18.04, so far here are my commands:

sudo apt-get install binutils libproj-dev gdal-bin

wget http://download.osgeo.org/geos/geos-3.4.2.tar.bz2

tar xjf geos-3.4.2.tar.bz2

sudo apt-get install build-essential

cd geos-3.4.2

./configure

make

I get this error

In file included from ../../../include/geos/geom/Geometry.h:25:0,
                 from IndexedPointInAreaLocator.cpp:17:
../../../include/geos/platform.h:112:2: error: #error "Can not compile without isnan function or macro"
 #error "Can not compile without isnan function or macro"
  ^~~~~
In file included from ../../../include/geos/geom/Coordinate.h:158:0,
                 from ../../../include/geos/geom/Envelope.h:25,
                 from ../../../include/geos/geom/Geometry.h:27,
                 from IndexedPointInAreaLocator.cpp:17:
../../../include/geos/geom/Coordinate.inl: In member function ‘bool geos::geom::Coordinate::isNull() const’:
../../../include/geos/geom/Coordinate.inl:38:10: error: ‘ISNAN’ was not declared in this scope
  return (ISNAN(x) && ISNAN(y) && ISNAN(z));
          ^~~~~
../../../include/geos/geom/Coordinate.inl:38:10: note: suggested alternative: ‘SNAN’
  return (ISNAN(x) && ISNAN(y) && ISNAN(z));
          ^~~~~
          SNAN
../../../include/geos/geom/Coordinate.inl: In member function ‘bool geos::geom::Coordinate::equals3D(const geos::geom::Coordinate&) const’:
../../../include/geos/geom/Coordinate.inl:77:21: error: ‘ISNAN’ was not declared in this scope
   ((z == other.z)||(ISNAN(z) && ISNAN(other.z)));
                     ^~~~~
../../../include/geos/geom/Coordinate.inl:77:21: note: suggested alternative: ‘SNAN’
   ((z == other.z)||(ISNAN(z) && ISNAN(other.z)));
                     ^~~~~
                     SNAN
Makefile:372: recipe for target 'IndexedPointInAreaLocator.lo' failed
make[5]: *** [IndexedPointInAreaLocator.lo] Error 1
make[5]: Leaving directory '/home/ralph/geos-3.4.2/src/algorithm/locate'
Makefile:391: recipe for target 'install-recursive' failed
make[4]: *** [install-recursive] Error 1
make[4]: Leaving directory '/home/ralph/geos-3.4.2/src/algorithm/locate'
Makefile:445: recipe for target 'install-recursive' failed
make[3]: *** [install-recursive] Error 1
make[3]: Leaving directory '/home/ralph/geos-3.4.2/src/algorithm'
Makefile:476: recipe for target 'install-recursive' failed
make[2]: *** [install-recursive] Error 1
make[2]: Leaving directory '/home/ralph/geos-3.4.2/src'
Makefile:367: recipe for target 'install-recursive' failed
make[1]: *** [install-recursive] Error 1
make[1]: Leaving directory '/home/ralph/geos-3.4.2'
Makefile:677: recipe for target 'install' failed
make: *** [install] Error 2

so following the docs does not seem to be working...

I do not know how to work around this error, any suggestions

are there any good write ups on how to install geos on ubuntu?


texts: How to split polygons in one layer into polygons of a different layer in ArcGIS Pro?

texts: I am using ArcGIS Pro and I have 2 layers that are both polygons. The first layer is called ELC, and the second layer is called SPU. My ELC polygons overlap my SPU polygons, sometimes the ELC layer overlaps multiple SPU polygons. I want to split my ELC polygons into the SPU polygon they overlap with. The reason I want to do this is so that I can get the total area of ELC polygons within each SPU polygon. I have tried using the split tool but it gives me errors when I try to do this for multiple polygons. Does anyone know how to do this?


texts: What is meaning of 'accuracy' and 'confidence' in Nominatim search results?

texts:
I have the search result for the string "12 street" from Nominatim in GeoJSON (Python geocoder library):

'geometry'  : {'coordinates': [-110.0121359, 53.2489107],
               'type': 'Point'},
'properties': {'accuracy': 0.225,
               'address': u'12 Street, Lloydminster, City of Lloydminster, Alberta, T9V 1V4, Canada',
               'city': u'Lloydminster',
               'confidence': 9,
               'country': u'Canada',
               ...
               } # Truncated for simplicity.

What do the the terms 'accuracy' and 'confidence' mean exactly, in terms of geocoding?


texts: Adding box to toolbox created with ArcPy (only visuals)?

texts:
I was wondering if there is a way I could easily include three parameters into a box (red box in the example) and so when there are a lot of inputs it's easier for a user to spot what data belongs to what sets of inputs.

Is there a way to give a title to the box ("Summary text" in red in the example) as well? In the example below, for example, I want to include param0, param1 and param2 into a box.

I'm hoping there is something built into arcpy rather than using third-party libraries.

import arcpy, os

class Toolbox(object):
    def __init__(self):
        self.label = "Toolbox"
        self.alias = ""
        self.tools = [Tool]
class Tool(object):
    def __init__(self):
        self.label = ""
        self.description = ""
        self.canRunInBackground = False

    def getParameterInfo(self):
        param0 = arcpy.Parameter(
            datatype = "Folder", direction = "Input",
            parameterType = "Required", name = "Summary_1",
            displayName = "Field 1")

        param1 = arcpy.Parameter(
            datatype = "String", direction = "Input",
            parameterType = "Required", name = "Summary_2",
            displayName = "Field 2")
        param1.type = "ValueList"
        param1.filter.list = ["Item 1", "Item 2", "Item 3"]

        param2 = arcpy.Parameter(
            datatype = "File", direction = "Output",
            parameterType = "Required", name = "Summary_3",
            displayName = "Field 3")

        param3 = arcpy.Parameter(
            datatype = "Boolean", direction = "Input",
            parameterType = "Optional", name = "Date_1",
            displayName = "Replace Date?")

        param4 = arcpy.Parameter(
            datatype = "String", direction = "Input",
            parameterType = "Optional", name = "Date_2",
            displayName = "Enter Date to replace")

        params = [param0, param1, param2, param3, param4]
        return params

    def execute(self, parameters, messages):
        return


texts: Connecting points (bus stops), which don't lie on lines (LINESTRING), to network?

texts:
I need to connect bus stops (points) to a network layer (OSM data). These bus stops do not lie directly on the lines (see screenshot) nor should their location be moved. I use PostGIS, pgrouting and QGIS and the network is already routable with source and target columns etc.

Mainly I want to do two things afterwards:

Getting the distances between the bus stops using shortest path analysis.
Creating isochrones with walking distances from the bus stop using the OSM network.

To get exact values it is necessary, that the routing 'starts' and 'stops' most closely to the bus stops. In many cases the closest existing node will be too far away to get exact values. But there should not be a routing to the actual point location of the bus stop. In my example on the picture you can see how the routing between stops should look like.

Is there a possibility to insert new nodes automatically into the network (LINESTRING) which are closest to the bus stops or is it possible to start routing on a kind of 'dummy point' which is set just for the query (similar to what the road graph plugin in QGIS does)?


texts: Decimal degree coordinates: how many decimal places to use for ...?

texts:
When using coordinates in decimal degrees, how many decimal points actually make sense?

For example, say one is producing a table of locations of base camps in a report, there's no point in displaying a number like 62.946916, -135.018639 when 62.95, -135.02 is sufficient. On the other hand when giving someone a waypoint for GPS navigation 62.947, -135.019 is more warranted (manually keyed in, why type 6 digits when 3 will do?), while for locating a geocache 4 positions is desired, and 6+ for a survey stake.

In other words, what is the easiest way to determine for any given purpose how many decimal places should be displayed?

There are technical answers at Algorithm for offsetting a latitude/longitude by some amount of meters, but they're not readily adaptable to a quick reference handbook style of usage. I'm looking for a handy chart to refer to like "the nth decimal place is X distance in ground units (at Y latitude)", or "the nth decimal place is X distance in ground units (at ___ major city/landmark)".


texts: Using use_spheroid=false in ST_DWithin of PostGIS

texts:
I have created 2 geography columns and indexed them using GIST. I am looking to use ST_DWithin for a larger dataset and to make it faster trying to use the use_spheroid=false.
SELECT pois.objectid, ST_Multi(ST_Union(polygon.coordinates)) as coords
              into table3
              FROM table1 as polygon
              INNER JOIN table12 as pois
              ON ST_DWithin(polygon.geog,pois.geog, 10* 1609.34, use_spheroid=false)
              group by pois.objectid
            limit 1;

but I get the following error,
ERROR:  column "use_spheroid" does not exist

Am I using the use_spheroid=false argument correctly


texts: How can I get Google Fusion Tables to work in QGIS?

texts:
I am using GDAL/OGR 1.10 and QGIS 2.0.1.

I am trying to load a google fusion table into QGIS using the method mentioned in this blog post. The only difference is that I am not declaring a email and password but using the OAuth refresh token I have generated. My VRT code is shown below.

<OGRVRTDataSource>
    <OGRVRTLayer name="ps">
       <SrcDataSource>GFT:refresh=<MY_VERY_LONG_TOKEN_IN_HERE> tables=1kQSkw-2cM_PU_jnu_T_PIksBkbyXm1feYaIooq8</SrcDataSource>
    </OGRVRTLayer>
</OGRVRTDataSource>

But if I then add this layer to QGIS I get an OGR error that said "OGR[3] error 1: HTTP error code : 400", followed by another "OGR[3] error 1: Fetching OAuth2 access code from auth code failed."

I thought that this would have worked because via the commandline I can use OGR to open the data source as shown in the image below.

I have also tried structuring my VRT layer like in the blog post also with no luck. I'm a bit stuck on how to make this work. Any suggestions?

Thanks
Ando


texts: How can I set multiple values to null in one raster and replace all other values with another raster?


texts: Raster TIFF file does not project onto QGIS basemap correctly

texts:
Forgive me if this is a simple fix, I am new to QGIS 3.8. I have a .tif file that when opened into QGIS is projected on the (0,0) lat/long coordinate of the base map and is very small. I have confirmed the CRS for both is the same. My hunch is the lat and long on the .tif file may not match or is inconsistent with the base map. I tried to change the Metadata>Extent to match that of the base map with no luck. I have attached the .tif file below. Does anyone have any thoughts on how to get this data to project onto the base map correctly? The data is located over the Antarctic Peninsula. The .tif range is between Latitude(-60, -75) and Longitude(285,310)

https://www.dropbox.com/s/b6ri8c4inu3lzlv/2016_foehn_melt.tif?dl=0


texts: OpenLayers + GEOExt + extjs. VS OpenLayers+jquery

texts:
Could some one maybe tell me if there are any real benefits of using OpenLayers + GEOExt + extjs. as oppose to using only OpenLayers and maybe jquery for client side webgis apps?

Edit: thanks for replying guys. I'm doing my graduation project which will be a web app for finding the best bus route to do a certain trip and route the user to the departure stop from the departure point and from the destination stop to the destination point. There are 3 layers one for the good looking map from osm, one from my roads containing the topology and one from the stops, i'll be using pgrouting and postgres+ postgis for the database.


texts: Exporting ESRI image Tiles designed for Export?

texts:
I have been trying to directly export ESRI tiles designed for offline use (http://www.arcgis.com/home/item.html?id=226d23f076da478bba4589e7eae95952) to tile package but all the methods I know (Create Map Tile Package georprocessing tool and ArcMap->File->Sharing->Tile Package) fail complaining about layer type is not supported. So I am wondering how I am supposed to "Export" these for offline use.

The end goal is to have some base map in a Mobile Map Package (MMPK) file that includes other data so that crew working in remote areas have some reference about locations.

I have seen the other posts, none really answers this question (or even the original question). Closest questions are:

Exporting ESRI World Imagery from World Imagery (for Export) map service? and
Using ESRI World Imagery (for Export)?


texts: Simplifying polygons in ArcGIS considering contains-one-point topology rule?

texts:
I want to simplify polygons using the "Simplify Polygon (Cartography)" tool in ArcGIS 10.2. Additionaly there is one point in each polygon and this condition should still be fulfilled after simplification. Therefore i created a topology, added the two feature classes and created a rule "Contains One Point".

When executing the simplify tool I selected "RESOLVE_ERRORS" (http://resources.arcgis.com/en/help/main/10.2/index.html#//007000000011000000), but the result isn't correct, there are polygons containing not exactly one point (as defined by topology).

Before:

After:

Have I done something wrong or is it a bug?


texts: Image Collection Statistics in Google Earth Engine?

texts:
How can I extarct mean of each image collection in google earth engine?
for example for this code:

var modis = ee.ImageCollection("MODIS/006/MOD13Q1")
.filterBounds(geometry)
.filterDate("2000-01-01","2001-01-01")
.select("NDVI");

print(modis)

var mod13 = modis.map(function(img){
  return img.multiply(0.0001)
  .copyProperties(img,['system:time_start','system:time_end']);
});


texts:
In version 2.18, when using the Reshape tool on multilines and polygons in 3D, the Z coordinate of new points are interpolated based on Z coordinates of existing points. For example, a line with "z" : 100,150 if I reshape Z stay : 100,115,130,135,150.

Now in QGIS version 3.4 all new points I create with the Reshape tool have z=0. Last example stays : 100, 55, 0, 75, 150, because in the New point z=0.

Is there an equivalent to the old Reshape tool in QGIS 3?


texts: Using modeled surfaces together to make predictions?

texts:
My question is about stacking multiple rasters of predicted surfaces, that have some of the same predictors.

For example,  if I make a predicted surface using all the terrain attributes as predictors, and then add that interpolated surface to a new stack with the original predictors (eg slope, aspect, TPI, and the new surface) will this violate any modeling principles? Is it commonly done?

It seems problematic because I have the original predictors in the stack, and then the new interpolated surface based on those predictors in the stack too. It is like a geospatial double jeopardy.


texts: Fully symbolize a vector field marker in PyQGIS

texts:
I'm designing a symbol for an external application that needs to be a vector field marker combined with a simple point marker, displaying a line at a certain angle around the center point. In PyQGIS, it's easy to instantiate the vector field class, as well as the field type (Polar), and angle units (Degrees). However, since it's polar, I need to set the length and angle attributes from the table of the layer I'm trying to display. I'm seeing in the API documentation for the QgsVectorFieldSymbolLayer that there are functions to set the X and Y attributes, which is necessary for a cartesian vector field marker, but irrelevant for a polar vector field marker. I need something like "setLengthAttribute()" and "setAngleAttribute()", much like what is possible in the QGIS GUI:

My relevant snippet of code currently looks like this:

    symbol = QgsMarkerSymbolV2.createSimple({})
    symbol.deleteSymbolLayer(0)

    symbol_layer = QgsSimpleMarkerSymbolLayerV2()
    symbol_layer.setSize(.5)
    symbol_layer.setColor(QColor("yellow"))
    symbol_layer.setOutputUnit(QgsSymbolV2.MapUnit)
    symbol.appendSymbolLayer(symbol_layer)

    symbol_layer = QgsVectorFieldSymbolLayer()
    symbol_layer.setColor(QColor("yellow"))
    symbol_layer.VectorFieldType(QgsVectorFieldSymbolLayer.Polar)
    symbol_layer.AngleUnits(QgsVectorFieldSymbolLayer.Degrees)
    symbol_layer. # lengthAttribute?
    symbol_layer. # angleAttribute?
    symbol.appendSymbolLayer(symbol_layer)

    renderer = QgsSingleSymbolRendererV2(symbol)

    self.mylayer.setRendererV2(renderer)

Any idea how to set the angle and length attribute?


texts:
I tried to use Tilecache with my mapserver/openlayers/postgis/cygwin configuration, but no success so far, I just got tired of keeping digging in the whole system to get that thing working.

So, could anyone please suggest some cache tools that instead can be used with mapserver/openlayers/postgis/cygwin? (...yeah, I have to use cygwin)

Any hints are welcomed,


texts: How to get the latitude & longitude of the 4 corners of the map being displayed

texts:
I am working on OpenLayers (ol3) and I want to get the Latitude, Longitude of the corners of the map being shown at a point. I tried using map.getView().getProjection().getExtent() and converting the points from 'EPSG:3857' to 'EPSG:4326' but I get -180,-85.05112877980659,180,85.05112877980659 irrespective of the zoom level of the map I have set.

I would like to know the map view corner points dynamically.


texts: Loading data in a Spatialite DB with custom CRS from QGIS

texts:
I created a custom CRS in QGIS and proceeded to create a Spatialite database with layers using this custom CRS. In order to do this, I edited spatial_ref_sys and spatial_ref_sys_aux tables in the database and added my custom CRS.
When I try to load data into this database using QSpatialite I get a Destination SRID isn't valid.
Loading through the DB Manager succeeds but the table doesn't get an SRID. My only option to load data is by using the standalone Spatialite GUI which works perfectly fine.

So my question is: is there a workaround in order to be able to successfully load data from QGIS into a Spatialite database with a custom CRS?

Note: I haven't messed around with QGIS srs.db file.


texts: I was wondering what is the easiest way to get XY geocaching coordinates (by a particular Area of Interest (AOI)) data in a spread sheet format that I could covert to a GIS layer?  On the Geocaching.com website you can search by different AOI, but I do not see a way to mass export XY coordinates.


texts: Cannot open gdb folder from external drive using geopandas

texts:
I am trying to open a .gdb folder from an external drive with geopandas:

import geopandas as gpd
data = gpd.read_file('/Volumes/SSD/data.gdb')

and I am getting the following error:

DriverError: '/Volumes/SSD/data.gdb' not recognized as a supported file format.

However, when I am trying the same on my local storage, geopandas can strangely read the folder just fine. Any ideas where this discrepancy in behavior might originate? I have observed this error with two different datasets so far.

EDIT:

I am on Mac OS 10.14.6.

I had a thought the files may have gotten corrupted when transferred to my external drive, but when copying back to my local drive I am able to open them just fine once again.


texts: Removing NaN from raster using ArcGIS Pro 2.0.1?

texts:
I have an image (RGB) with a lot of NaN pixels (black) and I want them to have no color for visualization on a map.

My raster is derived surface reflectance values from Sentinel 2. After atmospheric correction I stacked the blue, green and red band in R and exported the rasterStack (.tif) to ArcGIS Pro for making a map of the RGB raster on a basemap.

I have added a basemap that cannot be seen due to the black color and when I click the pixels I also see that they have No Data (nan)

But I cant seem to get rid of the black color. When I look in Symbology, background and NoData already is set to No Color. When I try to change the color of No Data, it has no effect.

How can I make my black NaN pixels no color/transparent using ArcGIS Pro?

Edit: I have added a picture showing the edge of my raster against the basemap, and the information I get when I click a black pixel.


texts: How do I open a shapefile in QGIS for upload to Garmin GPSmap 62?

texts:
I now have a polygon and have successfully raised 25 random sample points and saved as a shapefile. I have checked that GPS utility and the shapefile are set at EPSG:27700 BNG but when I try and open the sample points shapefile in GPS utility it says cannot open shapefile.

Can anybody help and is there an easy solution?


texts: Sample method returns empty featureCollection (GEE)

texts:
I have some troubles with the sample method from Google Earth Engine. On some images, it returns an empty feature collection (0 element).

Why this is happening and how do I deal with that ?

It's problematic to train a model on an empty dataset...

Here one illustration:

// Example on two images : one working, not the other
var image_ok = ee.Image("COPERNICUS/S2/20150627T102531_20160606T223605_T31RFN");
var image_empty = ee.Image("COPERNICUS/S2/20151228T002843_20151228T085259_T54HYD");
var VisImg = {"opacity":1,"bands":["B4","B3","B2"],"min":293.6805388296989,"max":2220.5270083401124,"gamma":1};

// Sampling
var training_ok = image_ok.sample({numPixels:1000});
var training_empty = image_empty.sample({numPixels:1000});

// Display
print("Ok: ", training_ok);
print("Empty: ", training_empty);

Map.addLayer(image_empty, VisImg, "Image empty");
Map.addLayer(image_ok, VisImg, "Image ok");
Map.centerObject(image_empty);


texts:
I did this 3D map with ThreeJS in QGIS:3D population density map

I almost done, but I need help to remove or fill perforations/hole.

Can anyone help me to do this?


texts: How do I get the coordinates of a click on vector feature/layer in OpenLayers2?

texts: I need to get the coordinate of the click when the user clicks on a vector feature on the OpenLayers map. The SelectControl only gives the feature that was clicked and not the coordinates of the click. Anyway to get at the coordinates of the click on a vector? I need to show AnchoredBubble at the point of the user click.