texts: Differences between Arc 10.0 and Arc 10.1 Layer files

texts: We have a program that cannot import Layerfiles from 10.1, but instead changes them to 10.0 layer files.  I'm wondering if there is any Layer functionality that is lost in this transition.  Right now only a couple of people are using 10.1 in our organization and they have been saving their layer files as 10.0, no issues.  But going forward if we all change to 10.1 will this automatic saving to a 10.0 version of a layer file lose anything that we would notice?


texts: Shows previous Atlas features grayed out

texts:
I try to make a set of maps using QGIS Atlas. My maps will consist of selected feature with it's attributes (using $id = @atlas_featureid, and that's working). Basically - districts in particular city with population, from the lowest to highest. But I want also to show the previous shown features grayed out/in different color/etc. Like here GeoGIF.

I've found that the solution could be a formula $id <= @atlas_featureid. But it's not working properly in my case. The problem I've found is that ID of features in my project does not correspond to the order of features I want to show.

In my project - the order of maps generated with atlas is based on population value, from lowest to highest. How to handle this problem?


texts: How do you convert to degrees and minutes from 8/9-digit lat/lon DMS code?

texts:
I'm trying to map the national bridge inventory, but the coordinates are stored in this 9-digit code for longitude (and 8-digits for latitude)

EXAMPLE:

Longitude is 81°5.8'
Code is 081054800

Page 9 of document found here explains:
http://www.fhwa.dot.gov/bridge/mtguide.pdf

The data can be found here:
http://www.fhwa.dot.gov/bridge/nbi/ascii.cfm

I've seen there are tools like Mimee, but I can't figure out how to run it on Ubuntu (tried to install correct libraries, etc.)

Is there maybe another way to do this?


texts: QGIS point array in JSON format to polygon


texts:
I received a layer package (lpkx) from a client. Finding out I couldn't open it in ArcMap led me to discover that this is a newer version of the software (ArcGIS Pro).

I managed to install that and license it. But once I opened it in ArcGIS Pro I quickly found out that the functionality to add proportional piecharts is not (yet) implemented in Pro, so I can't do the things I need to do.

Now it appears that the entire file from the client is useless, as I can't open it in ArcMap either.

Are there any possible fixes?


texts:
I'v got a detail map for an event as PDF an made an openlayers slippy map with Maptiler. I now need to place markers at the position of gps/gsm-trackers, which store their coordinate on a server.

The map's north isn't at top and i need some kind of formula (projection?) to translate the gps coordinates to the coordinates of my map.

Is this somehow possible, e.g. via. 2 points where I know the coordinates in both systems?

I don't expect a finished solution but maybe you can give a hint where to start, I've no idea how to do this.

Example:
GPS 48.26771, 16.36882 (somewhere in Vienna, Austria) in my projection is
0, 0

GPS 48.26673, 16.3696 is
100, 0


texts: CartoDB – How can I place labels from one layer on top of an overlapping layer?

texts: I have two overlapping polygon layers and I'd like to display the labels from the bottom layer on top of both layers so they're easily visible. Is there a way to do this with CartoCSS? Some sort of z-index for the labels?


texts: PostGIS difference between polygons - Ring Self-intersection?

texts:
I want to take ST_Difference or ST_Intersection to take the difference between two polygons.
SELECT a.id, a.geom, st_difference(a.geom, b.geom)
FROM liechtenstein.flaechen_norden a, liechtenstein.naturschutzgebiet b

From geometry a I want only these areas which are not in geometry b. But my code - whether ST_Difference or ST_Intersection - doesnt work, I get this error:
FEHLER:  lwgeom_difference: GEOS Error: TopologyException: Input geom 0 is invalid: Ring Self-intersection at or near point 505930 5252800 at 505930 5252800

I think the problem is that there are many polgons which have a hole in themselfes. How do I include these polygons so that I get the difference?


texts: Update field based on selected attribute from dropdown list

texts:
Is it possible to assign an attribute from a list (table) based on a chosen attribute (code) from a dropdown list field?

For the selected code, field "Sifra", update field "Opis" with something from the codelist table.


texts: Solving QGIS error Python Plugin is Broken?

texts:
I tried following the steps to build a python plugin for QGIS. (step 3 is where I am stuck at)

1) I used the Plugin Builder plugin to build the initial set of files. I call my plugin "PointWithinRadius". After using the plugin builder I have the files: init.py, pointwithinradius.py, resources.qrc, etc..(the usual files)

2) Next I ran pyrcc4.exe and was able to get resource.py from resource.qrc

3) Tried running pyuic4.bat as follows: (but get the errors as seen below)

![C:\Program Files\Quantum GIS Copiapo\bin>pyuic4.bat -d -o "C:\Users\stalmaki\.qg
is\python\plugins\PointRadius\ui_pointradius.py" -x "C:\Users\stalmaki\.qgis\pyt
hon\plugins\PointRadius\ui_pointradius.ui"
Traceback (most recent call last):
  File "C:\Python25\lib\site-packages\PyQt4\uic\port_v2\invoke.py", line 44, in
invoke
    exit_status = driver.invoke()
  File "C:\Python25\lib\site-packages\PyQt4\uic\driver.py", line 69, in invoke
    self._generate()
  File "C:\Python25\lib\site-packages\PyQt4\uic\driver.py", line 103, in _genera
te
    self._opts.pyqt3_wrapper, self._opts.from_imports)
AttributeError: Values instance has no attribute 'from_imports'][1]

I am unable to figure out the reason for this and searching for this error didn't give me much success either.

Finally when I try opening my plugin in QGIS, I get the following (see attached screenshot jpeg)

I want to build a sort of hello world plugin at first but these errors are not allowing me to get out of the starting blocks. Any help on solving would be great.

My development environment:
Windows 7 32 bit
QGIS 1.6.0
Python 2.5
Qt 4.7


texts: Building pyramids of virtual raster not permanent when file is moved

texts:
I have a folder saved under:

C:\QGIS\Active Layers\

In the "Active Layers" folder there our two folders; one called "Geotiffs" and the other "Virtual Raster".

In QGIS 2.18.20 I created a virtual raster with the .tif files contained in the "Geotiffs" folder and saved the virtual raster in my "Virtual Raster" folder. I created pyramids using the internal if possible option which created a .ovr file alongside my .vrt file. Everthing worked as expected.

Now if I move the "Active Layers" folder onto someone else's computer and save it in the exact same file path the pyramids don't work properly causing the virtual raster to disappear when I zoom in so far.

What have I done wrong and is there a way around it?


texts: Can QGIS automatically snap points to a polygon in a layer?  From that can the points then be put into a text delimited file in sequencial order?


texts: Dissolving Features from different Layers

texts:
I'm trying to dissolve polygon features from different layers by their FID. If two polygons have the same FID the result should be one single polygon. Maybe this drawing will help to make clear what I'm trying to get:

Here's how I was going to proceed:
First I performed a Union of the two layers joining only the FIDs. Thus I got a new layer with an attribute table containing two fields (called FID_Layer1 and FID_Layer2) with the initial FID of each layer and the value -1 if the polygon is from the respective other layer. Then I tried to add a new field combining those two fields but erasing the negative values and replacing them with the actual value. Based on this field I am going to perform a Dissolve hoping that it will give me the results I want.
In order to get the new field on which to perform the Dissolve I wrote a script in Python:

def FID_Union(inValue):
if !FID_Layer1! < 0:
    return !FID_Layer2!
else:
    return !FID_Layer1!

However, this doesn't work. The error message tells me that there's a syntax error in Line 2 but I can't figure out what's wrong.

So can anybody tell me what I have to change in the script to make it work?
Or is there another (easier) way to dissolve the polygons from different layers? I'd be happy about any suggestions.


texts: Calculating nearest distance to landscape element taking direction into account (i.e. slope, prevailing winds) using ArcGIS Spatial Analyst?

texts:
I am trying to build a layer of distance to relevant landscape elements in a watershed. To be relevant, a landscape element must be downstream with respect to the runoff direction.

Clearly, I want to assign to each point belonging to a regularly spaced points layer, the distance to the nearest landscape element, which is in the direction of the slope.

To do so, I have :

A polygon shapefile of « LANDSCAPE ELEMENTS » which are polygons of hedgerows, and forests (that are not distinguished).
A Digital Elevation Model
A shapefile of a regular point grid (50 m x 50 m, a total of 500 000 points )

Note : I'm using ArcGIS 10.1 with Spatial Analyst


texts:
I was playing around with the maputnik style editor and I wondered if it is possible to have styles that depend on the proximity of nearby features.
E.g. buildings within a certain proximity of water bodies having a different styling compared to other builds.


texts: Finding center point of MULTILINESTRING using PostGIS?

texts:
How to find the center point of MULTILINESTRING, which lies on the multiline?

I am using SQL in PostGIS


texts: ArcPy Calculate Null Values to Zero in Multiple Fields

texts:
I have a table that has Null values and I would like to insert a portion in my script to calculate these Null values to 0 in only a few of my columns in my table.

Here is my script so far:

    with arcpy.da.UpdateCursor(target, "ELEM_STUD") as cursor:
        for row in cursor:
            row[0] = None
            cursor.updateRow(row)

    del row

Target is my table, which is a parameter that is defined prior to this portion of my code.

The error that I am getting is in the first line, stating that 'ELEM_STUD' is not defined. That is the first of multiple fields I would like this script to run through.

I have also tried to add removeNull ("ELEM_STUD") after my cursor and haven't had any luck there either.


texts:
Is there a way to get the mymap even it was restricted on sharing.
var  src = "https://www.google.com/maps/d/kml?mid={id}&forcekml=1"

geoXml = new geoXML3.parser({
      map: map,
      infoWindow: infowindow,
      singleInfoWindow: true,
      afterParse: useTheData
     });
     geoXml.parse(src);

     function useTheData(doc) {
      var geoXmlDoc = doc[0];
      console.log(geoXmlDoc)
    }

Access to XMLHttpRequest at 'https://www.google.com/maps/d/kml?{id}&forcekml=1' from origin 'http://127.0.0.1:5501' has been blocked by CORS policy: No 'Access-Control-Allow-Origin' header is present on the requested resource.


texts: Zoom to Selected Feature for point layers in QGIS 2.0

texts:
I find the QuickFinder plugin for QGIS 2.0 very useful to quickly search, find, and select data. However, I'm under the impression that checking the "scale" option doesn't seem to do anything. I thought that perhaps it would change the scale of the map and zoom it to a certain level, but my results are the same whether or not this option is checked.

It seems as though Zoom to Selected Feature doesn't change the scale of the map either. I'm wondering if this might be why the scale feature for the Quickfinder doesn't do so. I'm trying to confirm if these are the general effects of these tools or if it is simply my implementation.

After finding this very old ticket: https://trac.osgeo.org/qgis/ticket/154 I still feel confused about the current status of this topic. Does anyone have any ideas of what I may be missing? Thanks in advance.


texts: Change color of polygon (country) in QGIS

texts:
I have a txt file added like a layer to QGIS 3.4.2. File number 1. It is an attribute table.

I have a shapefile (a map) added to QGIS 3.4.2. This shapefile has all the countries of the world.

Both files are joined by a key field.

In the file number 1 I have a list of countries:

Key file - Name country - Color

1 - Canada - #000000
2 - US - #FFCC00

How can I change the background of all the countries of the map with the field "color" of every country?


texts: Tiles from lower resolutions still visible in the map

texts: I'm using OpenLayers to display information on a tile grid with multiple resolution levels, with different tile sizes on each level. However, when zooming IN, the lower level (larger) tiles are still visible, but when zooming OUT, the higher level tiles (smaller) disappear. Is there anyway to only show the tiles for the particular resolution interval you're in?


texts:
Is there another use for "Calculate Geometry" aside from calculating area, length and points? I already read the help section ArcGIS Help (10.2, 10.2.1, and 10.2.2)  and it's very informative.

However I was told today that you can use calculate geometry to check that the coordinates of a feature are correct. If you have features with coordinates and they're displayed in the correct location, why would you need to use calculate geometry to check that. So although this makes not sense to me, I went ahead and used calculate geometry to check that the feature (points) was accurately displayed.

CALX & CALY are the new x, y coordinates using the calculate geometry option.

So my questions are:

1) Is this a valid/correct alternative use of "calculate geometry"?

2) If so, why are my new coordinates different from the original x & y?


texts: Extracting non-spatial data from imported KML file using R?

texts:
I have a KML file which was created using Google's My Maps.

The original file can be downloaded here:
Google My Maps

Using R, I can import this using the "readOGR" function of the rgdal library
This brings the KML file in as a SpatialPointsDataFrame (SPDF) - which i am calling asf52

In this SPDF, the spatial data is contained under @coords and is readily extracted into a dataframe using code like

df  <- data.frame(asf52@coords[,1:2])

However, I am struggling to come up with a way to neatly extract the the non-spatial data - contained under @data$Description - and turn it into a dataframe with a column for each variable.


texts:
I have the problem that the WMS layers are hidden from the basemap, I try to solve the problem with the zIndex but it does not work. Any other ideas?

WMS Layer

var two = L.tileLayer.wms('http://localhost:8080/geoserver/Routen/wms', {
        layers: 'Routen:engelbecken',
        type: 'FeatureCollection',
        format: 'image/png',
        zIndex: 1,
        transparent: true
        });

Basemap

var OpenStreetMap_Mapnik = L.tileLayer('https://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png', {
maxZoom: 19,
zIndex: 2,
attribution: '&copy; OpenStreetMap contributors',
});


texts: Reprojecting a large raster using Python and rasterio

texts:
I'm trying to reproject a raster that describes land cover over time in the Brazilian Amazon.  The raster can be downloaded at http://storage.googleapis.com/mapbiomas-public/COLECAO/2_3/CONSOLIDACAO/AMAZONIA.tif. It is roughly 2.1 gigabytes on disk. The metadata looks like this (there is one band per year):

{'count': 17, 'crs': CRS({'init': u'epsg:4326'}), 'dtype': 'uint8', 'affine': Affine(0.0002694945852358565, 0.0, -75.00007357655362,
       0.0, -0.0002694945852358565, 6.000027445691108), 'driver': u'GTiff', 'transform': (-75.00007357655362, 0.0002694945852358565, 0.0, 6.000027445691108, 0.0, -0.0002694945852358565), 'height': 90112, 'width': 123904, 'nodata': None}

You can visualize the data at http://mapbiomas.org/map.

The raster is originally in http://spatialreference.org/ref/epsg/wgs-84/; I'd like to reproject it to an equal area projection in meters, e.g. http://spatialreference.org/ref/sr-org/brazil-albers-equal-area-conic-wgs84/proj4/, at a resolution similar to the original (I think 30 meters would do).

My personal laptop has only 6 gigabytes of RAM. Can I reproject this raster using rasterio and Python?

I tried the following code, but it requires me to load the original raster values into RAM, and I get a memory error:

import numpy as np
import pyproj
import rasterio
from rasterio.warp import reproject, Resampling
from rasterio import crs, transform

def main(infile_path='./mapbiomass/AMAZONIA.tif',
         outfile_path='./mapbiomass/AMAZONIA_reprojected.tif',
         origin_lat=-11.188404,
         origin_lon= -58.338657,
         resolution=30.0):
    destination_crs_string = ' '.join(['+proj=aea +lat_1=10 +lat_2=-40 +lat_0=-25 +lon_0=-50',
                                       '+x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_def'])
    pyproj.transform(pyproj.Proj(init='EPSG:4326'), pyproj.Proj(destination_crs_string), origin_lon, origin_lat)
    origin_x, origin_y = pyproj.transform(pyproj.Proj(init='EPSG:4326'),
                                          pyproj.Proj(destination_crs_string),
                                          origin_lon,
                                          origin_lat)
    # See https://github.com/mapbox/rasterio/blob/master/examples/reproject.py
    dst_shape = (64, 64)  # Test
    dst_transform = transform.from_origin(origin_x, origin_y, resolution, resolution)
    dst_crs = crs.CRS.from_string(destination_crs_string)
    reprojected_mapbiomass_values = np.zeros(dst_shape, np.uint8)

    original_mapbiomass_raster = rasterio.open(infile_path)
    original_mapbiomass_values = original_mapbiomass_raster.read()  # Memory error
    reproject(original_mapbiomass_values,
              reprojected_mapbiomass_values,
              src_transform=original_mapbiomass_raster.transform,
              src_crs=original_mapbiomass_raster.crs,
              dst_transform=dst_transform,
              dst_crs=dst_crs,
              resampling=Resampling.nearest)
    pdb.set_trace()

    with rasterio.open(outfile_path,
                       'w',
                       driver='GTiff',
                       width=dst_shape[1],
                       height=dst_shape[0],
                       count=reprojected_mapbiomass_values.shape[0],
                       dtype=np.uint8,
                       nodata=None,
                       transform=dst_transform,
                       crs=dst_crs) as dst:
        dst.write(reprojected_mapbiomass_values)

if __name__ == '__main__':
    main()

Edit: I have

$ gdalwarp --version
GDAL 1.11.3, released 2015/09/16

and based on the suggestions in the comments I am trying

gdalwarp -tr 30.0 30.0 -s_srs EPSG:4326 -t_srs "+proj=aea +lat_1=10 +lat_2=-40 +lat_0=-25 +lon_0=-50 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs" AMAZONIA.tif AMAZONIA_reprojected.tif


texts: Error occured while running script in QGIS?

texts:
I am using QGIS version 2.18.9. I first create a graphical modeler and it is working. Then I export it to python script. While I am running the script I got following error massage.

  "NoneType' object has no attribute 'getitem See log for more detail".

Can you explain the reason?

My script is shown below:

##Testing=name
##drawingdxf=file
##cleanf=output vector
outputs_GRASS7V_IN_DXF_1=processing.runalg('grass7:v.in.dxf',drawingdxf,'1',True,False,True,False,False,True,['200000.000000,700000.000000,200000.000000,700000.000000'],0,None)
outputs_GRASS7V_CLEAN_1=processing.runalg('grass7:v.clean', outputs_GRASS7V_IN_DXF_1['output'],1,0.1,['200000.000000,700000.000000,200000.000000,700000.000000'],-1.0,0.0001,cleanf,None)


texts: Batch processing DEM rasters in QGIS/GDAL

texts:
I am attempting to add relief shading to map using .asc raster files, following this tutorial. It works well, but I have only tried it on single tile of raster data which is 20x20km. I want to generate the shading for a much larger area, which means multiple tiles.

I understand I can input GDAL commands in the Hillshade/DEM dialogue in QGIS, and the tutorial suggests this means I can perform a batch operation, but I don't know enough (read any) GDAL to do this.

I am using QGIS 1.9 on Ubuntu (though I do have a Windows installation as well).

Any ideas?

EDIT: I have written the following script, demraster.sh:

!/bin/sh
for f in *.asc
do
echo "Processing %f"
gdaldem hillshade %f -z 1.0 -s 1.0 -az 315.0 -alt 45.0
done

I have set it to be executable with chmod, and run it with sh demraster.sh but get a syntax error:

Syntax error: word unexpected (expecting "do")


texts: How to create a routable road network from Ordnance Survey data?

texts:
I've been following the "Beginners Guide to pgRouting" in an attempt to create an in-house application that provides drive-time isochrones.
I have a server, postgreSQL and PostGIS and pgRouting - no problem. I have drivers to I can talk to the database from ASP.Net pages - easy again.
I've got the data from OS too - imported shapefiles into the database - all easy.

Now I've run the bits of code in that article and the first two are easy again (and quick) so I've got my Node table and the road_ext view OK. However the third statement that creates the "Network" table is killing me, or more specifically the server. I have been trying with just the OS Grid tile "SO" which is not the most remote part of the UK by any means but still not going to have anything like as many roads as the London / Birmingham / Manchester area tiles surely. I'm running it on a quad-core box with a decent chunk of RAM, it's been running for 15 hours and not done yet. Box is not doing anything else, it's dedicated to the task. Even a little tile like the tip of Scotland (John O'Groats) took about 3 hours.
At this rate I'm looking at over a week to build a network which is ridiculous. Am I doing anything obviously wrong here or is it just going to take that long and I have to deal with it.

Ultimately I want to produce these polygons and overlay them on Bing Maps. If I have gone off-course above, could anyone please suggest an alternate dataset like OSM or something that's ready made I can use?

Thanks,
Michael


texts: PyQGIS code takes a lot of time to execute

texts:
I've written the Python code below in PyQGIS which loops through a ~30 MB shapefile. It adds an attribute with features which are calculated based on another attribute.
Is it normal that it takes around 20 minutes to fully execute the code?
Sometimes it doesn't execute it at all - QGIS just loads forever (waited over 30 minutes) until I close it.
The shapefile has ~130.000 features and 27 attributes.
My laptop has 16 GB RAM and an Intel Core i5-1135G7.
Code:
from qgis.PyQt.QtCore import QVariant
from qgis.core import QgsProject

layer = QgsProject.instance().mapLayersByName("677_5334")[0]
iface.setActiveLayer(layer)

caps = layer.dataProvider().capabilities()

if caps & QgsVectorDataProvider.AddAttributes:
    res = layer.dataProvider().addAttributes([QgsField("Stockwerke", QVariant.Double)])
layer.updateFields()

exp1 = QgsExpression('"HOEHEGEB"/3.5')

context = QgsExpressionContext()
context.appendScopes(\
QgsExpressionContextUtils.globalProjectLayerScopes(layer))

with edit(layer):
    for f in layer.getFeatures():
        context.setFeature(f)
        f["Stockwerke"] = exp1.evaluate(context)
        layer.updateFeature(f)

Updated the code (but doesn't work any faster unfortunately):
from qgis.PyQt.QtCore import QVariant
from qgis.core import QgsProject

layer = QgsProject.instance().mapLayersByName("675_5331")[0]
iface.setActiveLayer(layer)

caps = layer.dataProvider().capabilities()

if caps & QgsVectorDataProvider.AddAttributes:
    res = layer.dataProvider().addAttributes([QgsField("Stockwerke", QVariant.Double)])
layer.updateFields()

visited_index = layer.fields().indexFromName("Stockwerke")
attr_map = {}

exp1 = QgsExpression('"HOEHEGEB"/3.5')

context = QgsExpressionContext()
context.appendScopes(\
QgsExpressionContextUtils.globalProjectLayerScopes(layer))

for f in layer.getFeatures():
    context.setFeature(f)
    attr_map[f.id()] = {visited_index: exp1.evaluate(context)}
layer.dataProvider().changeAttributeValues(attr_map)


texts:
The navigation of a 3D qgis2threejs map for an Android tablet is not optimal. Test on 7 and 10 zoll.
Problem: If I want to zoom in or out with two fingers the zoomlevels are to big. It should be small zoom levels. In this time the map is very difficult to handle.
One finger for rotation and three fingers to move the whole map are OK.
Any idea to change this?


texts:
I was using the pgrouting plugin to display the shortestpath. I used the plugin as suggested by underdark to display the shortest path in qgis.

My source node is 283 and my target node is 172. When I used the pgrouting plugin to give me the shortest path I got one as attached in the screenshot.

However, when I queried the database of the graph which was generated by the tool osm2po. I found only one edge which was connected to node 283.

The result of the query is given in the figure below

My question is if you look at the results of the query, the edge is from node 279-283 and is one way as well with 283 being the target. Now since I used qgis to show the shortest path from 283-172. There is no starting from 283 as source. So my assumption was it shouldn't have shown any path. So I was wondering if there are some issues with the pgrouting plugin.

How come it gave me the shortest path when there is no edge with 283 as source? As I have explained there is only one edge connected to 283, but 283 is target and also the edge is one way from source to target.


texts: Geostatistical Analyst toolbar options not accessible

texts:
I am trying to perform cross-validation for a Kriging interpolation I did on my dataset. For this I wanted to use the Geostatistical Analyst Wizard in ArcGIS Desktop.

I have the Geostatistical Analyst extension installed and activated in ArcCatalog.

However, when trying to access any options of the Geostatistical Analyst toolbar, everything is grayed out and inaccessible.

I've tried to select a layer (both the original points layer as well as the Kriging-generated raster layer), and I've also tried to close and re-open the program multiple times - still the same story.

Would be glad if someone could help me find out what I'm missing here!


texts:
I am trying to produce a set of monthly solar radiation maps based on a monthly parameterization of the diffuse proportion and transmissivity (i.e. new values for both parameters each month). The PointsSolarRadiation function in arcpy produces a full year of monthly radiation values for each function call (i.e. 12 months all based on the same diffuse proportion and transmissivity values). I am trying to extract only the specific month from each of the twelve output files (i.e. extract the "April" data from the file run using the "April" parameters).

To do this I have a subsequent for loop which calls the individual files, identifies the appropriate data column in the file, extracts the data (using SearchCursor) then appends it to the placeholder. The append method fails on the second iteration, which does not happen if I run this with a simpler example of appending defined lists.

My code is as follows:

tableMonthHolder = []
for (j, value) in enumerate(diffuse):
    b = str(j)
    outFeatures = "c:/.../output" + b + ".shp"
    outTableFinder = "T"+b
    tableMonth = [row[0] for row in arcpy.da.SearchCursor(outFeatures,(outTableFinder))]
    tableMonthHolder = tableMonthHolder.append(tableMonth)

Any advice is greatly appreciated.


texts:
I have done an application with Geoserver for WMS Map rendering. In this application I have made a Dynamic Layer Group concept in which all the layers which I have checked on will be made as a layer group and request will be sent only once instead of single request for each layer. There are totally 89 layers. Another concept I have used in Load balancing. Totally 4 apache tomcat, in which application is in 1st tomcat and other 3 tomcat has geoserver. The WMS requests will flow in 3 geoserver
The following screenshot shows balancing of geoserver map rendering using "apache balance manager"

Eventhough I have made this many configurations, the getmap requests is very slow and it is not upto expected level. Since this application is going to used by more than 20 users, this speed is not sufficient. Geowebcache is not possible since map data is dynamic. We don't know what layer will be selected by the user.
My Queries :

I need to tune my application to work for 30 user simultaneously presently i having an HP server with XEON Processor + 4 GB RAM
ANy other peculiar configurations to be made for faster rendering of maps ? If so please suggest me !!


texts: Enabling public viewers and private editors of feature service on ArcGIS Online For Organizations?

texts:
In ArcGIS Online For Organizations, how does one publish a hosted feature service that is shared with Everyone for viewing, while restricting editing to specific users or groups, and different kinds of editing for different groups?

In the help at the end of Editor permissions for feature services under Configuring multiple tiers of access it says

  If you need to support multiple tiers of users, each requiring
  different operations, the recommended approach is to create one
  service for each level of user. For example, you could create a
  Professors service that might have all operations allowed. Along with
  that, you could create a less-privileged Students service with only
  Create, Query, and Upload operations allowed.

Which is great, that's what I need to do. However what is left out is how to reconcile the multiple Feature Services so that when the Professors create new features the Students see the new records, and vice versa.

The web browser management Create Service wizard allows you to use one feature service as a source for a new one, but only the data model is used (symbol definitions, field names & types, etc.) while the actual features are left behind, the new service is blank. This is the intended behaviour.


texts: Access PostGIS on Heroku instance from GeoServer

texts:
In Heroku, I have a PostgreSQL DB with all the credential to be accessed by outside (Host, Port, Database name, User, Password, URI). Then I connected pgAdmin with this database using the above credentials. Finally, I created a PostGIS extension for this database and I added a .shp file through the PostGIS Import/Export tool.

Now, I am trying to access this database as GeoServer Vector Data Source. When I enter all the db parameters (e.g Host, port and so on) in the "Connection Parameters" section I get an error:

Error creating data store, check the parameters. Error message: Unable to obtain connection: Cannot create PoolableConnectionFactory (FATAL: no pg_hba.conf entry for host "130.233.84.159", user "myuser", database "dbname", SSL off)

How can I solve this issue? Any hints?


texts: Stats on individual grid points nearest neighbours

texts:
I have shp of points where each point is the centroid of the individual cells of a 1*1km grid (~ 500 000 points). For each of these points, I have their coordinates and elevation.

What I want: for each of these points, I want to compute the mean, std dev, min, max and median elevation of the nearest 8 surrounding points.

What I did so far: I wrote an R script to identify each of the nearest points for each of the points of the grid. But considering the size of the dbf and the general slowness of R, it takes forever.

My tools: QGIS, ArcGIS, SAS, R or other free software if needed.

Does any one have a suggestion to do this (1- identify the direct neighbors, 2- compute summary stats from the elevation of these neighbors) in a more efficient way, using GIS?


texts: Failed to create a raster-brick from a netCDF file (.nc)

texts:
I am trying to build a raster-brick from a netCDF file (network Common Data Form)

leave this link for more information () but to make a brick (which is part of the package raster) throws me the following error:

  Loading required namespace: ncdf4
  Failed with error: ‘there is no package called ‘ncdf4’’
  Error: requireNamespace("ncdf4") is not TRUE

after that try installing the package ncdf4 but throws the following error when trying to install:

  Error, nc-config not found or not executable. This is a script that comes with the
  netcdf library, version 4.1-beta2 or later, and must be present for configuration
  to succeed.
  If you installed the netcdf library (and nc-config) in a standard location, nc-config
  should be found automatically. Otherwise, you can specify the full path and name of
  the nc-config script by passing the --with-nc-config=/full/path/nc-config argument
  flag to the configure script. For example:
  ./configure --with-nc-config=/sw/dist/netcdf4/bin/nc-config

  Special note for R users:

  To pass the configure flag to R, use something like this:
  R CMD INSTALL --configure-args="--with-nc-config=/home/joe/bin/nc-config" ncdf4
  where you should replace /home/joe/bin etc. with the location where you have
  installed the nc-config script that came with the netcdf 4 distribution.

How I can resolve this error?

I am new user in Linux, my Linux distro is Linux mint 17.1 and my version of R is 3.2.2


texts: How to access metadata from feature in a geodatabase with ogrinfo from command line?

texts:
Before loading a feature class from a File Geodatabase (.gdb) to PostgreSQL, I would like to check some general information such as the number of polygons, the bounding box or the attribute table.

The following GDAL versions is installed on my computer: GDAL 2.2.1, released 2017/06/23

Is there a way to do this in the command line? I would try it with ogrinfo, but did not find information about that in the internet, but some hints on this Cheat Sheet only provided information for ESRI Shapefiles (.shp)


texts: Using ArcGIS workspace defined as version in ArcPy?

texts:
I am trying to define my workspace as the current version that I am connected to.  The Python tool that I am running utilizes the arcpy.da.Editor() function, and I used a workspace variable that is currently hardcoded to reference the database connection file that I am operating out of.  Like this:

# Start an edit session
edit = arcpy.da.Editor(workspace)
edit.startEditing(False, True)
edit.startOperation()

*SCRIPT TOOL HERE*

# Stop editing session
edit.stopOperation()
edit.stopEditing(True)

This works just fine, but the problem is that since the workspace variable is hardcoded to reference the database connection file, it is referencing the sde.default version and does not account for the version that I am currently editing within, thus I get an error that states:

  workspace already in transaction mode

If I reconnect to the default version, I can rerun the tool without issue, but that is not ideal for obvious reasons.

Instead of the workspace variable being a reference to the database connection, I need to also add in there for the workspace to acknowledge the version that I am connected to.


texts: Seeking program similar to ArcCatalog in QGIS?

texts:
i'm qgis user.
There is one inconvenience with using qgis, it's file management.

when i used ArcGIS Desktop, ArcCatalog made me convenience for file management(especially file name change).

QGIS, provide Brower Panel, but it's function so simple.
Can not change file name....

so, when i need to change file name
i must go to folder, and change file name 4 or 5 times(shp, shx, dbf, prj....and so on)

Is there a way to avoid this inconvenience?


texts:
My goal is it to write a plugin, that automates the print composition. The idea is that you choose the scale and the paper size and based on that a map with legend, north arrow and block title is produced. The legend should only show the active layers. The finished map should not be automatically saved. Instead the print composer window should open so that the user can change things manually afterwards. I read the PyGIS cookbook and the PyGIS Programmer's Guide however there wasn't much information about the automation of the print composer.

What I want to know is, if it is generally possible to write a plugin that does this stuff? If it is possible, where do I find more information about automating and programming the print composer? Or does anyone know a tool that does something similar so I don't have to write my own tool?


texts: Counting number of points in 1 km grid covering surface of Earth?

texts:
I want to cover surface of earth with equal sized shapes and count number of points within each cell. shape of cell doesn't matters (triangle, hexagon ...), neither area, but equality of them is essential.

The problem is that I can not do this in a software like ArcGIS and need to calculate it using formula.

Here is what I am thinking:

convert geographic coordinates into Cartesian coordinates,
divide surface of the earth to equally distributed points
round all points to nearest such point.

Problem is step 2. How to distribute points equally?
Have any ideas?


texts:
I use QGIS print layout to create atlas with many features.
I am tring to find a way to show page number left, right, up and down of atlas feature.
Is there any way to insert these page numbers to print layout?
Example is on a picture below. I need to show page number 4, 6, 8 at atlas with page number 5 as is seen at atlas overview.
Unfortunately, script described at post below doesnt work.https://gis.stackexchange.com/questions/214300/how-to-determine-neighbouring-tile-ids-in-qgis


texts: Arcpy code only works when run in two chunks, not as one piece

texts:
I have a weird thing happening with my script. When I copy/paste it into the python cmd line within ArcMap, I get error messages. However, if I run the first part, then copy in the very last line by and run it by itself everything goes fine.
Here's my code:

import os
import arcpy
from arcpy import env
from arcpy.sa import *

arcpy.env.workspace = r'H:\Arc'
#dem = r'C:\Users\Lucy\Documents\GISthings\DCEW_VegStatistics\miscrasters\bare1m_fl_int'
dem = r'H:\Arc\dem_coarseres'
# create and store aspect
outAspect = Aspect(dem)

'''set up vegetation rasters'''
pftVeg = r'H:\Arc\veg_4classes'
# raster of n and s aspects
veg_n = Con((outAspect < 45) | (outAspect > 315), pftVeg)
veg_s = Con((outAspect < 225) & (outAspect > 135), pftVeg)

''' buffer around weather stations '''
featureToBuffer = r'H:\Arc\DCEW_WS.shp'
#outBuffer = r'C:\Users\Lucy\Documents\GISthings\DCEW_VegStatistics\bufferedWSveg.gdb\buff50m'
distanceField = "1000 meters"
# creates shapefile of buffers
# if you get an error here, try restarting computer to release schema lock
arcpy.Buffer_analysis(featureToBuffer, 'outBuffer', distanceField)

# Execute PolygonToRaster of buffers
valField = 'Field1'
assignmentType = "MAXIMUM_AREA"
priorityField = 'Field4'
cellSize = 0.00001  # 0.000001
arcpy.PolygonToRaster_conversion('outBuffer', valField, 'outRaster', assignmentType, priorityField, cellSize)
#outRaster = Lookup(Raster(outBuffer), 'Field1') #superfluous?

i = 4

# Process: Raster Calculator
arcpy.gp.RasterCalculator_sa("Con(\"outRaster\" == i,\"veg_s\")", 'temp_s')

The Raster calculator line is where everything gets fuzzy. My error message is this:

  Runtime error  Traceback (most recent call last):   File "",
  line 39, in    File "c:\program files
  (x86)\arcgis\desktop10.3\arcpy\arcpy\geoprocessing_base.py", line
  504, in      return lambda *args: val(*gp_fixargs(args, True))
  ExecuteError: ERROR 000539: Error running expression: rcexec()
  Traceback (most recent call last):   File "", line 1, in
     File "", line 5, in rcexec   File "c:\program files
  (x86)\arcgis\desktop10.3\arcpy\arcpy\sa\Functions.py", line 255, in
  Con     where_clause)   File "c:\program files
  (x86)\arcgis\desktop10.3\arcpy\arcpy\sa\Utils.py", line 53, in swapper
  result = wrapper(*args, **kwargs)   File "c:\program files
  (x86)\arcgis\desktop10.3\arcpy\arcpy\sa\Functions.py", line 239, in
  Wrapper     ["IfThen", in_conditional_raster,
  in_true_raster_or_constant]) RuntimeError: ERROR 000732: Input Raster:
  Dataset veg_s does not exist or is not supported  Failed to execute
  (RasterCalculator).

Does anyone know what could possibly be going on?

Edit in response to @FelixIP:

I'm using that syntax for my raster calculator because that's what model builder spit out for me. I'm still new to python, so I've been using it for training wheels. I tried your code with expr on one line and Rastercalculator on the other but I got the same error message. When I ran them by themselves everything worked perfectly. Does that make any sense? Ideally I'd like to put this within a for loop, so it would be good to not have to run it in two separate sections if possible.


texts:
I'm using QGIS 2.18. I need to create sector lights for navigational purposes on a map.

I have the light sector data as fields giving becon_id, start degree, end degree and color in a shapefile with over 500 buoys, beacons and lighthouses that need to be shown on the map. For each beacon, there can be many rows, each describing one light sector (for instance a white sector)

The end result should look about like below: The light sectors in the right colors, with the color marked as the character in the color field (R G W), and dotted lines from 100m to 1000m from the buoy/beacon/lighthouse.

This should most likely be created as a ruled-based symbol, but needs some python I guess?

Here below is an example of the shapefile data for a lighthouse (unfortunately not the above one) which has a green sector between 114 and 154 degrees, a white sector between 154 and 168 degrees, a red sector between 168 and 237 degrees, a green sector between 237 and 314 degrees, a white sector between 314 and 320 degrees, a red sector between 320 and 337 degrees (for some reason, 0 is not north but south):


texts: Table linked to dropbox file won't update

texts:
I have a table in my cartodb account linked to a file in my public dropbox folder. I have the settings in Cartodb such that the table syncs with every hour. When I look at the table in Cartodb, it says that it is syncing fine. However, it's not actually updating the table. When I make a change in the dropbox file, Cartodb doesn't register it.

Example: right now, my table says it synced 14 minutes ago and that my table was last updated 3 hours ago. This is incorrect. I have continually been updating my table over the last three hours, the last time being about 30 min ago.

Does anyone know why Cartodb may not be recognizing when my table is updated? Has anyone else had this problem?


texts: QGIS/Python script to generate random points many times

texts:
Can anybody provide a simple code to be run in the Python Console of QGIS (I am using 1.8 Lisboa) to generate X number of random points (say e.g. 5000) over a certain global shape of polygons (all Vector data). I then need to count the number of points per polygon and calculate the %. All that can be done via GUI but the catch is that I would need to repeat the procedure 100 times, i.e. generate 5000 random points 100 different times and then calculate the mean and standard deviation.
Does anybody have a simple code for this (repeating the random generation many times)?

I´m a total beginner with Python.


texts: Feature class imported table as points and now want to display

texts:
In ArcCatalog and ArcMap 10.

I have set up a feature class with all the fields I want to populate. I have data in a table format which includes wasting and birthing columns. I used the load data to get the data into the feature class. All the data came in successfully including the easting and northing. Now I want to display the points. I cannot find anywhere to display the positions and use the easting and northing data.

How can i display this data? Or set the easting and northing columns to be the position of the points?

The ultimate aim is to take a table and in one step out that data into a displayable feature class with set fields.


texts: Can FieldMappings be used to drop all attributes?

texts:
We have a script that regularly copies a bunch of feature classes out to shapefiles, keeping only the most important attributes (e.g. customer category) and dropping useless/private information (e.g. customer name).

Right now, it repeats the same set of steps without much variation. I am attempting to significantly simplify it by defining a function to do all those steps (the entire purpose of functions, of course.) Also, by introducing the arcpy.FieldMappings objects instead of verbose field mapping strings, it's possible to quickly understand what fields are coming from which data sources.

This works great for most of my feature classes.

def generate_field_mapping(inTable, inFieldList):
    fm = arcpy.FieldMappings()
    for field in inFieldList:
        vars()[field] = arcpy.FieldMap()
        vars()[field].addInputField(inTable, field)
        fm.addFieldMap(vars()[field])
    return fm

def subset_attributes(origFC, fieldList, outFC):
    fStart = time.time()
    name = origFC.split('.')[-1]
    # make feature class with only desired fields
    stringFieldMappings = generate_field_mapping(origFC, fieldList)
    tempfc = os.path.join(scratchWS, "temp_{}".format(name))
    arcpy.FeatureClassToFeatureClass_conversion(origFC, scratchWS, outFC, "", stringFieldMappings)
    print "  > done with {} ({} seconds)".format(name, round((time.time() - fStart), 2))

But, I have one feature class for which we want all of the fields dropped and I am not sure if I can accomplish this in the existing function workflow.

If I pass an empty list of fields, it breaks: FieldMap: Error in adding input field to field map
If I pass an empty ("") value, and put in a conditional to skip fieldmapping if there's no list, then all the original fields are passed.
If I pass OBJECTID (logic being, that would definitely be in the output), it breaks: Field mapping error: merge rule not valid for output field OBJECTID

I'm happy to just bypass the fieldmapping part of the function with a conditional, but I am not sure how to drop all the resulting attributes. (I don't want to DeleteField_management each one, that would be clumsy and I switched to FieldMappings to get away from that.) Any ideas for a workaround?


texts: Seeking GIS Solutions for land management?

texts:
I am working at real estate company.

We need to use a GIS solution for land management.

Which is better Open source or Esri Solution?

A big problem with Esri for me is its Licence. I need more licences for more employees and we need more Servers ,hardware.

Do you have a solution for this problems or Integrated between Esri and open source?


texts: Create Database View with ST_Within Query

texts:
I see several other questions similar but none that apply directly to mine.

I am using ArcGIS Desktop 10.6.1 with SQL Server running my Enterprise Geodatabase.

I have a very large point feature and a smaller polygon feature class.
My goal is to create a view to only show points that fall within the polygon.

The documentation I found for ST_Within on ESRIs page does not exactly help me.
In addition, when I tried to run a basic ST-Within query I get an error stating ST-Within is not a recognized function.


texts: How to extract coordinate from GML File in Java using geotools with STAX

texts:
I would like to extract all coordinates from my GML file and convert them in geometry type:
<geometrie>
                <gml:MultiCurve srsName="EPSG:3944" srsDimension="3">
                    <gml:curveMember>
                        <gml:LineString>
                            <gml:posList>1787810.8386000022 3214988.2920000004 0 1787810.7756000012 3214988.686999998 0</gml:posList>
                        </gml:LineString>
                    </gml:curveMember>
                </gml:MultiCurve>
            </geometrie>

....
<geometrie>
                <gml:Point srsName="EPSG:3944" srsDimension="3">
                    <gml:pos>1787697.0650999993 3214943.1966999993 0</gml:pos>
                </gml:Point>
            </geometrie>

I found the class GMLFilterGeometry in geotools but it's implemented in SAX and not STAX. It exists any over class like GMLFilterGeometry in STAX? Or another library can do that?


texts:
Do you know any tool that provides that? I need to manipulated GTFS file but I would like to have data available in any well-known format like XML or JSON.

EDIT

GTFS are Google Transit Feed Specification


texts: Counting Polylines on top of one another in ArcGIS Desktop?

texts:
I have a walking network, split into links.. Each link has more than one feature on top of each other which represents the number of times it was used within my network when plotting routes from destinations.

Is there a way I can merge each link back to one feature but with a field for the count on how many features were merged into one ?

Obviously this is one dataset, which over a 1,000 links so can I carry this out over the entire dataset... End result beings the network with a count field of the the number of time each link was used.


texts:
I've got an ST_Difference query runnning error free in PostGIS, and in PGAdmin I can see the output table but it's not viewable in QGIS - The table name appears in DB Manager but says "There is no entry in Geometry Columns".

What do i need to do to view the results?

Here's my query:

DROP TABLE IF EXISTS processing.newdifference;

CREATE TABLE processing.newdifference AS
SELECT 1 AS gid, ST_Difference(tdhab.geom, tdbuf.geom) AS "diff_geom"
FROM processing."tdhab", processing."tdbuf";

Any pointers?


texts: ArcGIS Pro Auto save file format (.save)

texts:
My ArcGIS Pro 1.4 has crashed and I noticed it saves an Auto-save file (.save) but I can't open it. The only file formats I can open are .aprx and .ppkx. Arc Pro doesn't suggest any auto-save version (like Microsoft Word for example) as I reopen the project neither.

Any methods to restore the auto-save?


texts: Labeling points with field value except those with a value of 0?

texts:
I am trying to label spot elevation points. Here is what I currently have:

Label expression:

Function FindLabel ( [OBJECTID], [CVC_1_Z], [CVC_2_Z] )
  FindLabel = [OBJECTID] & " ELEV(CVC_07)=" & [CVC_1_Z] & " " & " ELEV(CVC_18)=" & [CVC_2_Z]
End Function

Result:

How can I only label those points with non-zero Z values?


texts: Arcpy Create Polygon Geometry or Shapefile/FeatureLayer from WKT

texts:
I have a polygon feature class that i am looping through row by row. I want to take each row/record and create it into its own Polygon Geometry or Shapefile/FeatureLayer using WKT

PolygonGeometry Attempt:

summary = 'fc'
for row in arcpy.da.SearchCursor(summary, ["SHAPE@"]):
    wellknown = row[0].WKT
    arcpy.Polygon(wellknown)

I have tried

arcpy.Polygon(*[wellknown])

also

arcpy.Polygon(arcpy.Array([arcpy.Point(*wellknown)]))

Traceback (most recent call last):
  File "C:\Users\rzagha\Desktop\RTK\python\FreighMunicipality.py", line 121, in <module>
    hm = arcpy.Polygon(arcpy.Array([arcpy.Point(*wkt)]))
TypeError: __init__() takes at most 6 arguments (36467 given)

whats the correct syntax to create a polygon geometry from WKT and is there an easier way to create a featurelayer or shapefile from current row/record within the for loop like

arcpy.MakeFeatureLayer_management(row[0], "zone")

or

arcpy.MakeFeatureLayer_management(row[0].wkt, "zone")


texts: Google Docs CSV layer in ArcGis doesn't update automatically

texts:
I have added a layer to ArcGis Online from a Google Docs spreadsheet, following the instruction on this blog post

https://blogs.esri.com/esri/arcgis/2012/02/21/using-google-docs-in-your-arcgis-online-maps/

But the layer doesn't update automatically. Is this expected behaviour?


texts: Batch ENVI IDL PointCloudFeatureExtraction results in unprocessed areas

texts:
I have an IDL-ENVI code to batch PointCloudFeatureExtraction for lidar tiles across a large area, however, when run, the output results in lower areas of each tile as "unprocessed".

So my question is whether people have had this issue before or if you can see an issue in my code? The Point Cloud Feature Extraction works fine within ENVI+Lidar.

  e = ENVI(/HEADLESS)

  search_dir = 'E:\FilteredLidar4ppm\'
  filelist = file_search(search_dir + '*.las')
  task = ENVITask('PointCloudFeatureExtraction')

  TASK.DEM_GENERATE = 1
  TASK.DEM_FILTER_LOWER_POINTS = 10
  TASK.DSM_GENERATE = 1
  TASK.ORTHOPHOTO_GENERATE = 1
  TASK.TREES_GENERATE = 1
  TASK.TREES_MAX_HEIGHT = 6500
  TASK.TREES_MIN_HEIGHT = 300
  TASK.TREES_MAX_RADIUS = 2000
  TASK.TREES_MIN_RADIUS = 200

  foreach file, filelist, index do begin
  ; Open a las file
  pointcloud = e.OpenPointCloud(file, /CLOSE_PREVIOUS)
  ; Define inputs
  Task.INPUT_POINT_CLOUD = pointcloud

  Result = task.Validate(VALIDATION_EXCEPTION=msg)
  print, 'Executing Point Cloud Feature Extraction Task'
  Task.Execute
  pointcloud.Close

  endforeach

  END

I have tried:

Adding TASK.GENERAL_CLIP_MIN_HEIGHT = 0 to try and stop it being cut off - doesn't work.
the code on single file and it works fine, so perhaps there is something wrong with my looping foreach? Could settings from the previous loop not be removed?


texts: Using Field Calculator to auto-increment in ArcMap?

texts:
I am trying to edit the basic sequential numbers in a field script using Python in the Field Calculator.

I set the Parser to Python, and for my Pre-logic Script:

rec=0
name="manhole"
def autoIncrement():
 global rec
 pStart = 1
 pInterval = 1
 if (rec == 0):
  rec = pStart
 else:
  rec += pInterval
 return name, rec

And for the box below:

autoIncrement()

Would this not be how to go about this?

Essentially I just want "manhole1" then "manhole2" then "manhole3", etc.


texts: Installing QGIS for Viewing purpose only and not edit?

texts: I would like to find out if it is possible to install QGIS to be able to view and query data but not be able to make edits,is this possible.Can someone please advice how to do this.


texts:
I'd like to insert a row with a point geometry with pyspatialite. The sql looks like this:

insert into eqdata(origin_time, url, mag, lat, lon, depth, m_a, region, mt, the_geom)
                 values ('2013-11-10 07:39:59', 'http://geofon.gfz-potsdam.de/event.php?id=gfz2013walq', 4.300000, -19.840000, 70.050000, 45.000000, 'M',
                 'Near Coast of Northern Chile', 'None', geomfromtext('POINT(-19.840000 70.050000)', 4326));

This does not work and produces this error:

spatialite.dbapi2.IntegrityError: eqdata.the_geom violates Geometry constraint [geom-type or SRID not allowed]

If I execute the same SQL in the same database with spatialite_gui it works without any problem.

How can I handle this problem with pyspatialite?


texts: How can I create polygons for a 3D map from a DEM in grass/qgis

texts:
I have a SRTM 30m DEM of my area. I want to create a physical 3D terrain map from it (from felt, 3D printer, you name it). I have two questions about this process:

With grass r.contour I created contour lines from the DEM, and with v.generalize I smoothed the contour lines. Some of those contour lines still have gaps in them though. How can I make those contour lines smoother and interpolate in between so that two isolines that have the same elevation and have a 10m gap between them, are actually connected.

And more importantly, I would like to know how I can make polygons that correspond to the whole surface area with that height, so I can print those polygons (and cut them from felt, or input them in a 3D printer) to create polygons that I can stack on top of each other and create a 3D terrain map. How can I do this? I already have the isolines, the land outline and the bounding box outline for my area. I just need to somehow connect the isolines to the bounding box and land outline to create a polygon, for each isoline.

Any suggestions? I mainly use grass and qgis, but any other FOSS package is fine as well.


texts: MapXtreme compatibility with with .NET CORE

texts: In my office we are looking to use MapXtreme 9.0 for a new project we have. Is there a way to have MapXtreme support using .NET Core?


texts:
Is there a way to compare two attributes which can not be related spatially?
For example: I have .CSV file that has names of roads only, no coordinate, no geometry just names, and .SHP that contains the spatial info of those roads as well as their names. Can I compare the two files somehow to make sure the names in the .SHP are the same as the ones in the .CSV?
I'm looking for a method to do this comparison using QGIS or FME.


texts: Can you change the processing extent in environment settings using python to a shapefile?

texts:
I have managed to do this manually at Geoprocessing>Environments>Processing Extent and setting the shapefile in the drop down bar.

I was wonder if this could be automated with a python script?

My goal is to create a loop to change the processing extent prior to using spatial analysis tools on a various rasters that represent the areas of the shapefiles. This will keep their extent size minimal.

Edit (Current Code): (Not working)

for i in range (31):
... arcpy.env.extent = (r"C:\Folder\Shapefile"+str(i)+".shp")

Receiving a runtime error stating "Error in accessing environment "


texts:
I am trying to load geojson data with timedimension plugin as shown here with difference I am loading line objects. This is what I did:
   var timeSeriesGeoJSON =
  {
  "features": [{ "type": "Feature","id": -1549271008,"properties": {"id": 827793,"time": "2021-02-25T00:30:00","loading_percent": 4.584,"from_txt": "VTMI3_11","to_txt": "VVARNA11" },"geometry": {"type": "LineString","coordinates": [[27.338,41.483],[27.53,43.322]]}},
                {"type": "Feature","id": -1549271008,"properties": {"id": 827794,"time": "2021-02-25T01:30:00","loading_percent": 15.741,"from_txt": "VTMI3_11","to_txt": "VMA_IZ11"},"geometry": {"type": "LineString","coordinates": [[27.338,41.483],[25.863,42.122]]}}
    ],
    "type": "FeatureCollection"
    };

var map = L.map('map').fitBounds([[44,28],[42,27.5]]);

var timeDimension = new L.TimeDimension({
        period: "PT1H",
    });
map.timeDimension = timeDimension;

var player = new L.TimeDimension.Player({
    transitionTime: 100,
    loop: false,
    startOver:true
}, timeDimension);
var timeDimensionControlOptions = {
    player:        player,
    timeDimension: timeDimension,
    position:      'bottomleft',
    autoPlay:      true,
    minSpeed:      1,
    speedStep:     1,
    maxSpeed:      15,
    timeSliderDragUpdate: true
};
var timeDimensionControl = new L.Control.TimeDimension(timeDimensionControlOptions);

map.addControl(timeDimensionControl);

function getColor(d) {
    return d > 120 ? 'rgba(215,25,28,1.0)' :
        d > 100  ? 'rgba(237,110,67,1.0)' :
        d > 80  ? 'rgba(254,186,110,1.0)' :
        d > 60  ? 'rgba(255,232,164,1.0)' :
        d > 40   ? 'rgba(231,245,203,1.0)' :
        d > 20   ? 'rgba(183,223,227,1.0)' :
        d > 0   ? 'rgba(117,177,211,1.0)' :
              'rgba(44,123,182,1.0)';
  }

function style(feature) {
    return {
      weight: 10,
      opacity: 1,
      color: 'rgba(35,35,35,1.0)',
      dashArray: '',
      fillOpacity: 1,
      fillColor: getColor(feature.properties.loading_percent)
    };
  }

var timeSeriesLayer = L.geoJSON(timeSeriesGeoJSON
,{style: style}
);

map.fitBounds(timeSeriesLayer.getBounds())

var geojson = L.timeDimension.layer.geoJson(timeSeriesLayer);

geojson.addTo(map);

On display I get nothing. I tried console.log with geojson object. Object is populated and created, but it can't be displayed.


texts: Importing custom coordinate system into ArcGIS Desktop?

texts:
I have custom coordinate system (Korea).

If I want to import it, what should I do?

I'd like to know another way instead of (Define projection -> import -> one by one files). I'd like to know how to import a whole folder in C drive using ArcGIS Desktop.


texts: Limit to number of items using 'IN' SQL syntax in arcmap/arcpy?

texts:
I'm writing a python script that uses 'NOT IN ('a', 'b', ...)' statements on very large lists (like comparing 500 items in one feature class against >10k items in another) and am getting weird unexpected results (sometimes it works, sometimes it doesn't and it corrupts the geodatabase...?!?).

I've looked and can't find any ArcMap/arcpy documentation on how large a list can be using IN/NOT IN syntax and was wondering if anyone knew. Using ArcGIS 10.3.1 with 64-bit python.

Update: All of the data in question is stored in file geodatabases.


texts: Applying OTB mean shift smoothing to a folder of rasters using Python 2

texts:
I am very new to Python, and am having trouble coming up with a loop script that will allow me to apply a smoothing process to a folder of rasters. This is the script I am using for a single file right now:

import otbApplication

MeanShiftSmoothing = otbApplication.Registry.CreateApplication("MeanShiftSmoothing")

MeanShiftSmoothing.SetParameterString("in", "C:/Users/riddella/Desktop/Raster_Clip.tif")

MeanShiftSmoothing.SetParameterString("fout", "C:/Users/riddella/Desktop/Raster_Clip_Smooth.tif")

MeanShiftSmoothing.SetParameterInt("spatialr", 4)

MeanShiftSmoothing.SetParameterFloat("ranger", 256)

MeanShiftSmoothing.SetParameterFloat("thres", 0.1)

MeanShiftSmoothing.SetParameterInt("maxiter", 100)

MeanShiftSmoothing.ExecuteAndWriteOutput()

How can I reconstruct this script to read files from a folder, rather than just giving it a single file input and output? Apologies for my lack of base knowledge of Python, I only started working with it this week.


texts:
this is my first post to this forum. I hope to start participating actively.

I have a table with several records (aprox 200.000), each with its own coordinates. And I need to identify the "Polygon" where each point belongs (I got another table with the vertex for earch polygon) . I am working with Excel and Access at this point. Thanks for the help!

(the polygons doesn’t have holes and doesn’t overlap with other polygon)

I know there is an algorithm that checks if the number of times a line cross another line from the desired point to another point is even, the point is inside polygon, if the number is odd, the point is outside the polygon. But I don’t know how to implement this. Here are some reference links with the algorithm:

http://alienryderflex.com/polygon/

Here is a sample:

I got a table similar to this:

Point_ID / Coordinates
1 / -70.559349,-33.411913
2 / -70.579004,-33.407542
etc. (aprox. 200.000 points)

And I got the following table defining the boundaries of each polygon:

Polygon_ID / Polygon_Sequence / Coordinates
1 / 1 / -70.567331,-33.408803
1 / 2 / -70.584755,-33.414249
1 / 3 / -70.586987,-33.409692
1 / 4 / -70.585527,-33.404175
1 / 5 / -70.571537,-33.398084
1 / 6 / -70.567331,-33.408803
2 / 1 / -70.567245,-33.409047
2 / 2 / -70.564756,-33.416570
2 / 3 / -70.552998,-33.418576
2 / 4 / -70.551796,-33.408474
2 / 5 / -70.560122,-33.407542
2 / 6 / -70.562525,-33.407327
2 / 7 / -70.567245,-33.409047
etc...

So, this is the output I need (in bold):

Point_ID / Coordinates / Polygon_ID
1 / -70.559349,-33.411913 / 2
2 / -70.579004,-33.407542 / 1
etc...

Thanks a lot!
Alejandro


texts: Can client-side caching for GeoServer / GeoWebCache be customized?

texts:
Right now client-side caching appears to be working with its default configuration, using the following headers:

Response
HTTP/1.1 200 OK
Content-Type: image/png
Content-Disposition: inline; filename=States-US_HistCounties.png
Transfer-Encoding: Identity
Cache-Control: no-cache
ETag: 29e79cc2e4fe182660936eb2ff1f5aab
X-Frame-Options: SAMEORIGIN
Last-Modified: Thu, 20 Dec 2018 17:30:02 GMT
geowebcache-tile-bounds: -7924991.091503907,5117000.420810547,-7915207.151884766,5126784.360429689
geowebcache-crs: EPSG:900913
Server: Jetty(9.2.13.v20150730)
geowebcache-gridset: EPSG:900913
geowebcache-cache-result: HIT
geowebcache-tile-index: [1238, 2571, 12]
geowebcache-layer: xxx:xxx

I have found that the downside to this caching mechanism is that an extra request gets made back to the server to check if anything has changed since the last request(?) for this layer.

My goal is to not have that extra HTTP request, and instead immediately serve the cached content. This particular layer is static and will never change once cached.

Is there any way to configure that on a layer-by-layer level? I was thinking something along the lines of Cache-Control: public might be useful.

I have browsed the GeoWebCache documentation and tried to google articles, but I have been unable to find anything that talks about this.


texts: Control Order of Layers in QGIS Merge Vector Tool?

texts:
It may be that the answer to my question is the same as this one: How to sort the layer list in Multuple selection window, i.e. that there is way to do so in QGIS 2.18 but there will be in the 3.0 release. If so, I apologize for the redundant question. Still, I'd like to be absolutely sure there's not another solution to my problem.

I need to merge a few vector layers in QGIS to make, essentially, a "bullseye" around a specific point which I can then use to visualize the intensity of events around that point within each concentric perimeter.

Merge vector layer seems simple enough except for the fact that I can't control the order in which these layers are rendered/merged, so some of the smaller features are ending up underneath the bigger ones.

Is there any way to control the rendering/merge order of layers within the "Merge Vector Layer" tool in QGIS?

If any of this helps:

I am using QGIS 2.18.14
I am running Windows 10
"Control feature rendering order" under the style tab of the properties window doesn't affect the order in which layers appear in the list of features available for merging in the Merge Vector Layer tool (this was my first attempt at a fix)


texts: Extract time-series rainfall from multiple point with date information

texts:
I would like to have a CSV file which contain a time-series of 30 minutes rainfall from last 10-days before Landslide occurs from certain (landslide) location. I used NASA GPM IMERG data for the rainfall, and try to write a GEE script (compiling from others too) as below:

// Define a point of interest as a landslide location.
var POI = ee.Geometry.Point(111.19165, -8.08875); // Adjust the coordinate
Map.centerObject(POI, 9);

// Date when landslide occurs
var lsevent = new Date('2017-11-29'); // Adjust date period with landslide event
var start = new Date(lsevent.getTime() - 10*24*60*60*1000); // 10-days before
var end = new Date(lsevent.getTime() + 1*24*60*60*1000); // 1-day after
print(start);
print(end);

// Import NASA GPM IMERG 30 minute data.
var imerg = ee.ImageCollection('NASA/GPM_L3/IMERG_V06');
var imergHH = imerg.filterBounds(POI)
                   .filterDate(start, end)
                   .select('precipitationCal');

// Add rainfall accumulation into map
var precip = imergHH.select('precipitationCal').sum();

// Rainfall vis parameter
var palette = [
  '000096','0064ff', '00b4ff', '33db80', '9beb4a',
  'ffeb00', 'ffb300', 'ff6400', 'eb1e00', 'af0000'
];
var precipVis = {min: 0.0, max: 1000.0, palette: palette, opacity:0.5};

Map.addLayer(precip, precipVis, "10-days rainfall", false);
Map.addLayer(POI, {palette:"#ff0000"}, "Landlside location", true);

// Create a function that takes an image, calculates the mean over a geometry and returns
// the value and the corresponding date as a feature.
var timeSeries = imergHH.map(function (image) {
  var imergdate1 = image.date().format('yyyy-MM-dd-hh-mm');
  var value = image
    .clip(POI)
    .reduceRegion({
      reducer: ee.Reducer.mean(),
      scale: 30
    }).get('precipitationCal');
  return ee.Feature(null, {value: value, date: imergdate1});
});

var accumulation = imergHH.map(function (image) {
  var imergdate2 = image.date().format('yyyy-MM-dd-hh-mm');
  var value = image
    .clip(POI)
    .reduceRegion({
      reducer: ee.Reducer.sum(),
      scale: 30
    }).get('precipitationCal');
  return ee.Feature(null, {value: value, date: imergdate2});
});

// Create a graph of the time-series.
var graphTS = ui.Chart.feature.byFeature(timeSeries,'date', ['value']);
print(graphTS.setChartType("LineChart")
           .setOptions({title: 'NASA GPM IMERG 30-minute rainfall time-series',
                        vAxis: {title: 'Rainfall estimates (mm)'},
                        hAxis: {title: 'Date'}}));

// Create a graph of the accumulation.
var graphAcc = ui.Chart.feature.byFeature(accumulation,'date', ['value']);
print(graphAcc.setChartType("LineChart")
           .setOptions({title: 'NASA GPM IMERG 30-minute rainfall accumulation',
                        vAxis: {title: 'Rainfall estimates (mm)'},
                        hAxis: {title: 'Date'}}));

// Export the result to Google Drive as a CSV.
Export.table.toDrive({
  collection: timeSeries,
  description:'LS_20191129',
  folder:'GEE',
  selectors: 'date, value',
  fileFormat: 'CSV'
});

// End of script

So far, I got the time-series Chart and CSV, and I happy for that as my first goal.
Next step:
(1) I would like to have a Chart and CSV for the accumulation data from start to end of the period. The graphAcc (2nd chart) doesn't seem to give me the results I want. Did I missed something in the script?

(2) If I have a CSV file with information on ID,Lon,Lat,Date
How to integrate this data as input? So the script above can directly read the list and do rainfall extraction and export to multiple CSVs or append to the list. Any reference or example?

For reference, attached is the data from Chart 1, and I calculate the accumulation manually, and total rainfall for start to end is 558.557 milimeters.


texts:
I have an application built in Angular which will invoke a RESTful service in back-end application which in turn will call the ArcGIS REST API.

For: a feature will be created in the Angular application which will call addFeature (RESTful Service) in the back, which in turn will call the ArcGIS REST API which is for now:

https://services7.arcgis.com/02tfVoY2YvL3RJ26/ArcGIS/rest/services/iconlayer/FeatureServer/0/addFeatures

The back-end service would receive the data in GeoJSON format and the REST API is not accepting it. I guess it's expecting an ArcGIS format which looks like this:

[
  {
    "geometry" : {"x" : -118.15, "y" : 33.80},
    "attributes" : {
      "OWNER" : "Joe Smith",
      "VALUE" : 94820.37,
      "APPROVED" : true,
      "LASTUPDATE" : 1227663551096
    }
  },
  {
    "geometry" : { "x" : -118.37, "y" : 34.086 },
    "attributes" : {
      "OWNER" : "John Doe",
      "VALUE" : 17325.90,
      "APPROVED" : false,
      "LASTUPDATE" : 1227628579430
    }
  }
]

I'm now thinking of building a mapper which maps from GeoJSON to ArcGIS JSON but need help with how to map field wise.  Any guidance on this please?

I've also seen a thread (Converting ArcGIS Server JSON to GeoJSON?) saying the ArcGIS REST API accepts data in GeoJSON format but it doesn't seem to be the case with the API I'm calling.  Does the ArcGIS API accept GeoJSON?  If so, I don't have to do the mapping.


texts: How to use Route Task in Network analysis tools

texts:
I am using esri network analysis tool and I just want to draw the path between two clicked point on the map (by using route task)

I use the example code block in the following link
http://help.arcgis.com/en/webapi/silverlight/samples/start.htm#Routing
which is included in the ESRI ARCGIS silverlight API sample page and I only change the route task and Dynamic layer URLs... In the end I always receive the error message as "Routing error: 'stops' parameter is invalid"

Can anyone help me?


texts: Which QGIS processing tool to use to select non-shared features in two layers?

texts:
I have two layers each of which have hundreds of polygons. I want to create a new layer, which is reflective of the features which the two layers don't share.

I've been using the QGIS difference and the symmetrical difference tools, and even though the map displays right, the feature count is totally off.

For example, the first layer has 3,000 features and the second has 1,000. The new layer should have 2,000. And, again, even though it displays just fine with the difference tool, the feature count is like 2999 on the new layer. I don't get what's happening. I need the feature count to be accurate because I need the CSV file.


texts: How to mark a particular point in MapPanel?

texts:
There is a component called MapPanel, which displays OpenStreetMaps.

I want to mark some places on the displayed map, without changing the underlying map. Those marked places are identified by latitude and longitude.

In order to implement this, I wrote following class, which should convert latitude and longitude to pixel coordinates.

@Override
public Point convert(double aLatitude, double aLongitude, Rectangle bounds) {

    final Point pointInMapUnits = mapPanel.computePosition(new Point.Double(aLatitude, aLongitude));

    Point mapCenter = mapPanel.getCenterPosition();
    Dimension mapSize = mapPanel.getMapSize();

    final Rectangle panelSize = bounds;
    final Point2D.Double panelCenter = new Point.Double(panelSize.x + panelSize.getWidth()/2,
            panelSize.y + panelSize.getHeight()/2);

    final double ourObjectXFromMapCenterInMapUnits = pointInMapUnits.getX() - mapCenter.getX();
    final double ourObjectYFromMapCenterInMapUnits = pointInMapUnits.getY() - mapCenter.getY();

    final double ourObjectXFromMapCenterInPixels = convertToPixelsX(ourObjectXFromMapCenterInMapUnits, panelSize,
            mapSize, mapPanel.getZoom());
    final double ourObjectYFromMapCenterInPixels = convertToPixelsY(ourObjectYFromMapCenterInMapUnits, panelSize,
            mapSize,mapPanel.getZoom());

    final Point result = new Point((int)(panelCenter.getX() + ourObjectXFromMapCenterInPixels),
            (int)(panelCenter.getY() + ourObjectYFromMapCenterInPixels));

    return result;
}

But it doesn't work as following screenshots illustrate.

Step 1: Map is centered at the point in question (132.18381, 43.35514). That place is marked by a red dot.

Step 2: Now we scroll the map so that the marked place is located to the left of the position in step 1.

The blue dot shows where the red dot should be.

Step 3: Now we scroll the map so that the marked place is at the top.

Again the blue dot shows the expected position of the red dot.

Here's a video of the process.

You can find the code here.

How can I change it, in order for the dot to always be displayed at the point (132.18381, 43.35514) ?


texts: An effective way to count points / cluster using hex bins in QGIS?

texts:
Is there a best practice in QGIS for summarizing point data in small cells over a large area?
Background/Problem:
In QGIS, it is "possible" to create a grid of hex cells over an area using the extent of a layer.  Then one can summarize the points which fall in that polygon grid.  However this is very computationally wasteful to construct for points that cover a large extent, but which need to be analyzed in relatively small cells.  For example, if points that cover an entire region need to be binned at a 100-500M resolution - this will leave the user waiting for hours while QGIS constructs the grid, then a while longer for it to do the summary stats.
In ArcGIS Pro, I can [mis]use the Optimized Outlier Analysis tool, which allows me to determine the grid type, size and distance band, ignoring other stats and returning only polygons drawn and populated with the counts of the points.  This gives me the same sort of output (summarized hex cells) but at a tiny fraction of the time.
Is there a workflow/plugin in QGIS that will get me to a similar output?
I have tried the Hotspot Analysis plugin in QGIS, but A.) it appears abandoned and B.) it requires statistics fields as forced, not optional, and C.) pysal is tricky to get working correctly with QGIS on Windows.
UPDATE: Description of the Outlier Analysis Tool as requested, from the documentation:
Optimized Outlier Analysis executes the Cluster and Outlier Analysis (Anselin Local Moran's I) tool using parameters derived from characteristics of your input data. Similar to the way that the automatic setting on a digital camera will use lighting and subject versus ground readings to determine an appropriate aperture, shutter speed, and focus, the Optimized Outlier Analysis tool interrogates your data to obtain the settings that will yield optimal analysis results. If, for example, the Input Features dataset contains incident point data, the tool will aggregate the incidents into weighted features. Using the distribution of the weighted features, the tool will identify an appropriate scale of analysis. The classification type reported in the Output Features will be automatically adjusted for multiple testing and spatial dependence using the False Discovery Rate (FDR) correction method.


texts: Open the object attribute dialog for a selected feature with the Python QGIS API

texts:
I want to open the attribute dialog programmatically for a selected feature. I need to achieve the same effect as if the user clicks on the "object identify" action ( CTRL + SHIFT + I).
I searched the QGIS API docs for QGIS 3.2 but I was unable to find a solution.

I found the classes QgsAttributeFormInterface and QgsMapToolIdentifyFeature, but I was not able to use them in a useful way.
I also tried to trigger the action from the toolbar, but I guess I need some kind of click event with x,y coordinates for this to work.

Does anyone know if you can open the attribute dialog from a plugin programmatically?

Andi


texts:
I'm a newbie with GeoData Processing and Map Making. I'm doing a city map for someone and I've got the vector layers in with the polygons and the lines, but I need the street names. How do i make a vector layer that shows the street names?

Some additional Specs: I'm using QGIS with OSM data. I also run a Mac with OS X 10.11.6 aka "El Capitan". I do not have access to a PC.


texts: Accessing the spatial index of a layer (shapefile) with PyQGIS

texts:
I built a manual index on a shapefile, through the process "Create spatial index", which can be found in the process toolbox > general vector tools.
This process creates a persistent spatial index, a .qix file.
I want to optimize the query: "How many polygons of a layer are contained within one or more larger polygons that I select" using PyQGIS.
Both layers are based on files in polygon shapefile format. The layer of the polygons that will be evaluated has more than 300.000 entities.
I can build and use a custom index with this code:
pry=QgsProject.instance()
parcelas=pry.mapLayersByName('parcelas')[0]
radios=pry.mapLayersByName('radios_censales')[0]

tiempo_inicial = time.time()
index = QgsSpatialIndex() # Spatial index
index = QgsSpatialIndex(parcelas.getFeatures())
duracionSI=time.time()-tiempo_inicial  #  6.816403388977051 sec

csi=0
lisIds=[]
tiempo_inicial = time.time()
for i in radios.getSelectedFeatures():
   bbx= i.geometry().boundingBox()
   listIds=index.intersects(bbx)
   for j in listIds:
        if i.geometry().contains(parcelas.getFeature(j).geometry()):
            csi=csi+1
duracionSI=time.time()-tiempo_inicial

But what I really want to do use leverage the existing .qix file directly, so I don't have to spend six seconds compiling the index every time.  How can I do this from PyQGIS?
Note: I am aware that the spatial index file must be updated if geometries are modified, removed or added.


texts: Using Merge GeoProcess task in ArcGIS Engine with C#?

texts:
Having the hardest time getting the correct syntax when adding layers as gp parameters to call the Merge_management GeoProcessing Task.  Below is just one of MANY, MANY ways I have tried with no luck.  Can anyone help?

I've successfully called many GeoProcessing tasks via code in ArcEngine, but this is the first one to take an array ( ie. [layer1,layer2,etc]) as one of the parameters.

Intial Attemp:

    public void RunMergeGeoProcess(ILayer SourceLayer1, ILayer SourceLayer2, string ResultFeatureClassName, string GeoDBAbsFileName)
    {

        IFeatureLayer featureLayerSource1 = (IFeatureLayer)SourceLayer1;
        IFeatureLayer featureLayerSource2 = (IFeatureLayer)SourceLayer2;

        _geoProcessor.SetEnvironmentValue("Workspace", GeoDBAbsFileName);
        //_gpParams.Add(featureLayerSource1);
        ////_gpParams.Add(featureLayerSource2);

        IWorkspace refWS = ((IDataset)featureLayerSource1.FeatureClass).Workspace;
        string dbPath = refWS.PathName;

        IDataset myDS1 = featureLayerSource1.FeatureClass as IDataset;
        string mysDS1 = dbPath + "\\" + myDS1.Name;
        IDataset myDS2 = featureLayerSource2.FeatureClass as IDataset;
        string mysDS2 = dbPath + "\\"+ myDS2.Name;

        string myResult = dbPath + "\\" + ResultFeatureClassName;

        string param1 = "[" + mysDS1 + ", " + mysDS2 + "]";

        _gpParams.Add(param1);
        _gpParams.Add(myResult);

        try
        {
            _geoProcessor.OverwriteOutput = true;
            _geoProcessor.Execute("Merge_management", _gpParams, null);
            //returnMessages(_geoProcessor);
        }
        catch (Exception ex)
        {
            MessageBox.Show("There was a GeoProcessing Error." + ex.ToString());
            returnMessages(_geoProcessor);
        }

        _gpParams.RemoveAll();

    }

Working Version:

    public void RunMergeGeoProcess(ILayer SourceLayer1, ILayer SourceLayer2, string ResultFeatureClassName, string GeoDBAbsFileName)
    {

        IFeatureLayer featureLayerSource1 = (IFeatureLayer)SourceLayer1;
        IFeatureLayer featureLayerSource2 = (IFeatureLayer)SourceLayer2;

        _geoProcessor.SetEnvironmentValue("Workspace", GeoDBAbsFileName);

        IWorkspace refWS = ((IDataset)featureLayerSource1.FeatureClass).Workspace;
        string dbPath = refWS.PathName;

        IDataset myDS1 = featureLayerSource1.FeatureClass as IDataset;
        string mysDS1 = dbPath + "\\" + myDS1.Name;
        IDataset myDS2 = featureLayerSource2.FeatureClass as IDataset;
        string mysDS2 = dbPath + "\\" + myDS2.Name;

        string myResult = dbPath + "\\" + ResultFeatureClassName;
        string param1 =  mysDS2 + ";" + mysDS1  ;
        _gpParams.Add(param1);
        _gpParams.Add(myResult);

        try
        {
            _geoProcessor.OverwriteOutput = true;
            _geoProcessor.Execute("Merge_management", _gpParams, null);
            //returnMessages(_geoProcessor);
        }
        catch (Exception ex)
        {
            MessageBox.Show("There was a GeoProcessing Error during the Merge Operation." + ex.ToString());
            //returnMessages(_geoProcessor);
        }

        _gpParams.RemoveAll();

    }


texts: MapProxy GetFeatureInfo: passing additional parameter to QGIS Server

texts:
I have MapProxy Version 1.12 configured and working very well with QGIS Server 3.10.

I need to use cached layers suporting GetFeatureInfo requests. But I also need to pass to additional parameters to QGIS Server:

WITH_GEOMETRY=TRUE&
FI_POINT_TOLERANCE=10&

These additional parameters should only be added to the request for GetFeatureInfo requests.

But when I use a MapProxy cached layer, the additional GetFeatureInfo parameters are not passed back to the QGIS Server.

If I issue a GetFeatureInfo request against Mapproxy:

http://mapproxy.qgis.demo/mapproxy/service?
SERVICE=WMS&
VERSION=1.1.1&
REQUEST=GetFeatureInfo&
FORMAT=image%2Fpng&
TRANSPARENT=true&
QUERY_LAYERS=geotuga_costanova_short_house_numbers&
LAYERS=geotuga_costanova_short_house_numbers&
TILED=true&
INFO_FORMAT=application%2Fjson&
WITH_GEOMETRY=TRUE&
FI_POINT_TOLERANCE=10&
X=235&
Y=91&
WIDTH=256&
HEIGHT=256&
SRS=EPSG%3A3763&
STYLES=&
BBOX=-36975%2C120827%2C-36900%2C120902

The request if forwarded to QGIS Server as:

http://qgis.demo/postgresql/geotuga/public/costanova/cgi-bin/qgis_mapserv.fcgi?
transparent=True&
layers=short_house_numbers&
bbox=-36975.0,120827.0,-36900.0,120902.0&
width=256&
height=256&
x=235&
y=91&
query_layers=short_house_numbers&
info_format=application%2Fjson&
format=image%2Fpng&
srs=EPSG%3A3763&
request=GetFeatureInfo&
version=1.1.1&
service=WMS&
styles=

The additional parameters do not reach the server.

My source is configured with featureinfo: true in mapproxy.yaml:

  geotuga_costanova:
    type: wms
    wms_opts:
      featureinfo: true
      legendgraphic: true
    req:
      url: http://qgis.demo/postgresql/geotuga/public/costanova/cgi-bin/qgis_mapserv.fcgi
      transparent: true

How can I have all GetFeatureInfo requests going back to the server with additional parameters?

This is related to this Custom Vendor Parameters for WMS and Caches discussion. GetFeatureInfo requests are not cached, so it maybe easier to support.


texts: Convert grib (grb2) precipitation band to visible PNG

texts:
I'have been trying to convert band 7 (precipitation) of a grib file (you can download it if you like) to a visible PNG image which I can load via KML into a google maps viewer I have.

The precipitation of the grib displays just as I'd like to see it inside quantum gis if in Layer properties I choose to render band 7 as single band gray with a pseudocolor color map.  Furthermore specifying 0 for nodata values.

I have tried using this command with gdal (version 1.9.2-2~precise4):

gdal_translate -a_nodata 0 -of PNG -b 7 grib.grb2 grib.png

but the resulting grib.png 's size is just 295 Bytes and it is just transparent when opened in a viewer.

I need someone to explain to me how to do this and why.  And if you know other command-line tools which work to convert this file, that'll help, too.


texts: Widget Expand is not functioning correctly

texts:
I have placed 2 widgets (baseMapGallery, LayerList) inside a small map in a rather large popup. To keep it from crowding out the view, I am trying to use an Expand to launch it from a small icon button in the top-right corner of the map view. As you can see, there are 2 maps in this popup; each with its own separate div/view, independent of the main map view in the foreground.

Currently, when I try to launch the widgets from the icon button (Expand), it attempts to open in the main map view to the right of the popup. Also, there is nothing showing (no basmaps, no layer list). Here's my current code for the left map panel with the widgets, expand, and a handler:
 //Left panel map
    var popupDiv1 = document.createElement("div");
    popupDiv1.classList.add("mapViewLeft");
    popupDiv1.style.marginBottom = "10px";

    var popupView1 = new MapView({
        container: popupDiv1,
        map: map_left,
        actions: [],
        //center: view.center,
        center: [lon, lat],
        zoom: view.zoom,
        //scale:
        //    view.scale *
        //    2 *
        //    Math.max(view.width / 250, view.height / 250),
        constraints: {
            rotationEnabled: false
        },
        ui: {
            components: ["zoom"]
        }
    });

    var layerList_left = new LayerList({
        view: popupView1,
        listItemCreatedFunction: function (event) {
            const item = event.item;
            if (item.layer.type !== "group") {
                // don't show legend twice
                item.panel = {
                    content: "legend",
                    open: false
                };
            }
        }
    });

    bgExpandLayerList_2 = new Expand({
        view: popupView1,
        content: layerList_left.container,
        expandIconClass: "esri-icon-layer-list"
    });

    var basemapGallery_2 = new BasemapGallery({
        view: popupView1,
        container: document.createElement("div")
    });

    bgExpandBaseMapGallery_2 = new Expand({
        view: popupView1,
        content: basemapGallery.container,
        expandIconClass: "esri-icon-basemap"
    });

    expandHandle_a = watchUtils.pausable(bgExpandLayerList_2, "expanded", function (newValue, oldValue) {
        if (newValue === true) {
            expandHandle_a.pause();
            setTimeout(function () {
                expandHandle_b.resume();
            }, 100);
        } else {
            expandHandle_a.resume();
        }
        if (bgExpandLayerList_2.expanded) {
            bgExpandLayerList_2.collapse();
        }
    });

    expandHandle_b = watchUtils.pausable(bgExpandBaseMapGallery_2, "expanded", function (newValue, oldValue) {
        if (newValue === true) {
            expandHandle_b.pause();
            setTimeout(function () {
                expandHandle_a.resume();
            }, 100);
        } else {
            expandHandle_b.resume();
        }
        if (bgExpandBaseMapGallery_2.expanded) {
            bgExpandBaseMapGallery_2.collapse();
        }
    });

    popupView1.ui.add(bgExpandLayerList_2, "top-right");
    popupView1.ui.add(bgExpandBaseMapGallery_2, "top-right");

I've tried removing the handler and just adding a single widget, but I am still seeing the same problem (opening in the wrong view). Any ideas as to what I am doing incorrectly here? The widget works fine by itself (see right panel map), just not with the Expand. Note: I have a baseMapGallery and LayerList widget in the main view as well that are both working. I have named the variable differently for each. Is there a restriction on adding more than one of the same widget Expand in a web app?


texts: How to show Two Labels of One feature In Geoserver…? (Part 2)

texts:
This is how my SLD 1.0.0 rule looks like

  <sld:Rule>
      <sld:Title>Reser Type</sld:Title>
      <sld:PolygonSymbolizer/>
      <sld:TextSymbolizer>

        <sld:Label>
          <ogc:PropertyName>fp_no</ogc:PropertyName>
          <ogc:PropertyName>reser_use</ogc:PropertyName>
        </sld:Label>

        <sld:Font>
          <sld:CssParameter name="font-family">Arial</sld:CssParameter>
          <sld:CssParameter name="font-size">12</sld:CssParameter>
          <sld:CssParameter name="font-style">normal</sld:CssParameter>
          <sld:CssParameter name="font-weight">normal</sld:CssParameter>
        </sld:Font>
        <sld:LabelPlacement>
          <sld:PointPlacement>
            <sld:AnchorPoint>
              <sld:AnchorPointX>0.0</sld:AnchorPointX>
              <sld:AnchorPointY>0.5</sld:AnchorPointY>
            </sld:AnchorPoint>
          </sld:PointPlacement>
        </sld:LabelPlacement>
        <sld:Fill>
          <sld:CssParameter name="fill">#f42f08</sld:CssParameter>
        </sld:Fill>
      </sld:TextSymbolizer>
    </sld:Rule>

My question is Is there any way we can display these 2 labels above and below and NOT in one line ?

Please refer this question for more details.

Note 1: this is not a duplicate question read before reacting !!!!

Note 2: I tried \n it does not work


texts: Export Thematic Raster Files algorithm in eCognition not resulting in Raster output

texts:
I am trying to export a classification from eCognition in raster format and have tried both the Export Thematic Raster Files and Export Classification View algorithms and get a generic File as an output.
I specified the .img as the output format and tried all three options for the Export Type and all resulted in the same generic "File".


texts: Replacing NULL values, loop skips records using Python, OGR, Rtree?

texts:
I have a point layer with several NULL values and based on a near analysis with another layer, I would like to replace those NULLs. However, with the code I use  it only does a couple and then stops. Is there something wrong with my for or if loops? I am sure that more values fall within the distance band of 1000 that I apply.

inFile0 = r'C:\data\areas.shp'
inData0 = driver.Open(inFile0, 0)
areaLayer = inData0.GetLayer()
inFile1 = r'C:\data\points.shp'
inData1 = driver.Open(inFile1, 1)
pointLayer = inData1.GetLayer()

index = r.tree.index.Index(interleaved = False)
for i in range(0, areaLayer.GetFeatureCount()):
    areaFeat = areaLayer.GetFeature(i)
    geometry = areaFeat.GetGeometryRef()
    xmin, xmax, ymin, ymax = geometry.GetEnvelope()
    index.insert(i, (xmin, xmax, ymin, ymax))

distance = 1000
for point in range(0, pointLayer.GetFeatureCount()):
    pointFeat = pointLayer.GetFeature(point)
    if pointFeat.GetField('CLASS') is None:
        geometryP = pointFeat.GetGeometryRef()
        xmin, xmax, ymin, ymax = geometry.GetEnvelope()
        searchEnvelope = (xmin - distance, xmax + distance, ymin - distance, ymax
        + distance)
        for area in list(index.intersection(searchEnvelope)):
            areaFeat = areaLayer.GetFeature(area)
            geometryA = areaFeat.GetGeometryRef()
            if geometryP.Distance(geometryA) <= distance:
                classType = areaFeat.GetField('FEAT_TYPE')
                distance = geometryP.Distance(geometryA)
                pointFeat.SetField('CLASS', classType)
                pointFeat.SetField('DISTANCE', distance)
                pointLayer.SetFeature(pointFeat)
    else:
        continue

FIXED: the loop reassigns a value to the variable 'distance' which should not happen, over time this would reduce the distance band breaking the loop.


texts: Projections when using google maps and rest service on openlayers

texts:
I am trying to use a rest service on an openlayers map that has 2 google maps layers.  The google map layer is set to it's own custom projection (which I had to code into the library by hand, but I believe is the standard EPSG of 900913)..  However the layer that I would like to sit on the google map is in Oklahoma State Plane North.

var arcrest = new OpenLayers.Layer.ArcGIS93Rest("Rest Layer", arcserver_rest, {
                            LAYERS: 'show:18,19,20,21,29,30,31',
                            IMAGESR: '2267',
                            BBOXSR: '2267',
                            FORMAT: 'png',
                            TRANSPARENT: true
                        },
                        {
                        isBaseLayer: false,
                        units: 'ft',
                        maxExtent: new OpenLayers.Bounds(2034036.00510709,  136689.343284666, 2069119.14841628, 148575.857284666)
                        }
                        );

Why is this projecting into the Congo, and not Oklahoma?


texts: How do I reach a layer inside a GroupLayer?

texts: I am having problems in ArcObjects. I am struggling to get to the individual (feature-)layers inside a GroupLayer. Does anyone know how to do that ?


texts: MapProxy - Is it possible to get uniform WMS GetfeatureInfo info_format?

texts:
This is my very first mapproxy tagged question, so I hope it makes sense.

I am developing a WebGIS client in my spare time, and now I am mainly focusing on my first type of OGC service: WMS.

Everything worked quite fine so far except one big issue: I could not request a uniform format for any GetFeatureinfo. This is due to the specification (version 1.3.0) not having a default MIME type for the info_format. Everytime I receive is either a WMS exception or HTML data that display very differently from service to service.

I have to fix this problem and a quick search led me to MapProxy for the first time and I am very excited about it. Especially, this part of the documentation where it reads:

  MapProxy supports XSL transformations for more control over feature information. (...) You can use XSLT for sources to convert all incoming documents to a single, uniform format and then use outgoing XSLT scripts to transform this format to either HTML or XML/GML output.

Please can anybody confirm this is one (the?) way to go for my problem, before I start to invest more effort during my spare time to learn how it works?

I must admit I still have to figure out what MapProxy is, but I guess I could even use it to speed up my WMS requests which are super slow (I am currently using CORS-ANYWHERE in my development env), but this is another argument...


texts: Formatting OGR VRT file to convert CSV to KML?

texts:
I am attempting to get a nice KML for simple review of POI/geocoding from a CSV. I a using ogr2ogr and a VRT. You can see the VRT file below, with two attempts to get this right.

Basically, the issue is using ogr2ogr -f KML ... with the VRT w/o renaming results in a lot of "SimpleData" tags in the KML.

I'd like to get some CSV fields as the KML Placemark "Name" and "Description", concatenated fields, actually. -dsco NameField=RegionName option in ogr2ogr, but I don't see a way to concatenate fields into name, e.g Street, City, State..

In the .vrt below, the SrcSQL field results in no records. It does open, but with layer info. The Field tag (when not commented out) does result all records with in the new alias, but the "name" is the only field that comes through, it drops all other fields.

Is there a way to concatenate fields using ?

Of course, perhaps ogr is not capable of concatenation ...

<OGRVRTDataSource>
    <OGRVRTLayer name="roofing-poi2">
        <SrcDataSource>roofing-poi2.csv</SrcDataSource>
        <GeometryType>wkbPoint</GeometryType>
        <LayerSRS>WGS84</LayerSRS>
        <GeometryField encoding="PointFromColumns" x="Longitude" y="Latitude"/>
<!--         <Field src="StreetAddress" name="name" /> -->
        <SrcSQL>SELECT StreetAddress as name</SrcSQL>
    </OGRVRTLayer>
</OGRVRTDataSource>


texts:
I might have tomatoes on my eyes, but how can I generally activate Python macros in QGIS 3.12.

Before, I used QGIS 3.8 and had no problems at all with that...


texts: I convert file to GTFS（using transitfeed）and I haven't create shape.txt

texts:
I using transitfeed by python to convert my file to GTFS

but I haven't connect shape.txt and the trip.txt

and it show

  [the shapes with the following shape_ids aren't used by any trips：1]


texts: Python code for generating random points not working!

texts:
I am trying to create some code to generate 10 random points within a vector layer's districts:

layer=iface.activeLayer()
import ogr
import sys
import random
from shapely import wkb, geometry

Moz_no=layer.featureCount()

for i in xrange(Moz_no):
  district_tmp=layer.getFeatures(i)
  district=wkb.loads(district_tmp.GetGeometryRef().ExportToWkb())

  random_point_count=0

  sys.stdout.write("\r")
  sys.stdout.flush()

  while random_point_count<11:
    maxx, maxy, minx, miny = district.bounds
    random_point=geometry.Point(minx_+(random.random()*(maxx-minx)), (miny+ (random.random()*(maxy-miny))))

However, when I run the code, it gives me this error:

TypeError: QgsVectorLayer.getFeatures(QgsFeatureRequest request=QgsFeatureRequest()): argument 1 has unexpected type 'int'

Could someone please help me, and explain what is going on?


texts:
I search the config (.ini) file for QGIS 2.18 on my pc.
For QGIS 3.4.8 the file is stored under this path
C:\Users\AppData\Roaming\QGIS\QGIS3\profiles\default\QGIS\QGIS3

I can´t find this file for QGIS 2.18.

Where is this file normaly stored?