texts:
I am trying to run mosaicking using the script below:
saga_cmd grid_tools "Mosaicking" -TARGET_DEFINITION 0 -GRIDS "/media/geodev/Data2/test_data2/2019-06/2019-06-19/Flight_Plan_B_2/clip/IMG_0036_3_gcp.sdat;/media/geodev/Data2/test_data2/2019-06/2019-06-19/Flight_Plan_B_2/clip/IMG_0045_3_gcp.sdat;/media/geodev/Data2/test_data2/2019-06/2019-06-19/Flight_Plan_B_2/clip/IMG_0047_3_gcp.sdat" -NAME "Mosaic" -TYPE 8 -RESAMPLING 0 -OVERLAP 0 -BLEND_DIST 1.0 -MATCH 0 -TARGET_USER_SIZE 0.111078 -TARGET_USER_FITS 0 -TARGET_OUT_GRID "/media/geodev/Data2/test_data2/2019-06/2019-06-19/2019_06_19_Flight_Plan_B_2_1_3_v6.sdat"

But I am getting the error below:
Error: grid: memory allocation failed [1879.20mb]
Error: executing tool [Mosaicking]

I am using Ubuntu 16 with 64GB RAM, SAGA 2.3.1 from a repository shared here.


texts:
I am trying to replicate my erosion model in ArcGIS for Desktop in QGIS 2.4.

I tried to calculate the upstream distance along the flow path for my agricultural field using QGIS. I did the equivalent in ArcGIS with the tool Flow length and it is simple and effective.

I see two potential tools that can deliver a similar result in QGIS, r.flow (GRASS) or SAGA. The first tool, r.flow gives me a flow length path but the result make no sense and some output cannot be saved. The second tool flow path length with SAGA, I don't understand what is the seed and I don't see the option to calculate the upstream distance.

Arcgis Flow length (Value 0 to 840 meters, black to white):

QGIS r.flow (GRASS) (Value 0 to 6.5, black to white):

Here my input for the tool r.flow (GRASS) in QGIS 2.4

LiDAR elevation, I did the treatment Fill Depression.
Aspect raster.
Barrier, I made a raster contains all 0 value. (I'm not very sure here, I didn't find so much information for this input).
Calculate for upstream length.
Others parameters are the default ones.

Suggestions for the appropriate tool to use are welcome.


texts: I want to know what happen when we switch between the baseLayers of an openLayers Map, I mean when we switch the baseLayer, does it use the layer.setVisibility() to false for all the remaining layers or it uses map.removeLayer() for all others


texts: Previously I had Windows 7 and I recently installed Windows 10 (64-bit) on my computer. Upon installation of QGIS ( newer /older versions) the installation completes taking more time than it usually did and when I open it, an error dialog box pops up. The error message is attached. It says re-installing may fix the problem but reinstalling didn't help.


texts: How to build an artificial 'dam' across a river?

texts:
I have a DEM.
Created all files related to watershed including accumulation, direction, streams etc...
I'm stuck creating a dam.
I request you to please let me know how to build an artificial 'dam' across a stream/river.

Please suggest any software/methods to do.


texts:
I'm attempting to upgrade to a new version of GeoWebCache from 1.3.  I'm trying to upgrade to version 1.7.2.  I blew away the geowebcache directory under the Tomcat webapps and also deleted the geowebcache directory under Tomcat work.  I put the new WAR file in the Tomcat webapps directory and can not get it to start correctly.  I've been looking all over trying to see if anyone ran into this issue with no luck.  The server has Java 1.6 installed on it and according to the GWC documentation it supports 1.5 or later.

Here are the snippets from the stack trace.

org.springframework.beans.factory.BeanCreationException: Error creating bean with name 'restFilterInvocationInterceptor' defined in ServletContext resource [/WEB-INF/acegi-config.xml]: Cannot resolve reference to bean 'restFilterDefinitionMap' while setting bean property 'objectDefinitionSource'; nested exception is org.springframework.beans.factory.CannotLoadBeanClassException: Error loading class [org.geowebcache.security.RESTfulDefinitionSource] for bean with name 'restFilterDefinitionMap' defined in ServletContext resource [/WEB-INF/acegi-config.xml]: problem with class file or dependent class; nested exception is java.lang.UnsupportedClassVersionError: org/geowebcache/security/RESTfulDefinitionSource : Unsupported major.minor version 51.0 (unable to load class org.geowebcache.security.RESTfulDefinitionSource)
    at org.springframework.beans.factory.support.BeanDefinitionValueResolver.resolveReference(BeanDefinitionValueResolver.java:328)
    at org.springframework.beans.factory.support.BeanDefinitionValueResolver.resolveValueIfNecessary(BeanDefinitionValueResolver.java:106)

...

Caused by: org.springframework.beans.factory.CannotLoadBeanClassException: Error loading class [org.geowebcache.security.RESTfulDefinitionSource] for bean with name 'restFilterDefinitionMap' defined in ServletContext resource [/WEB-INF/acegi-config.xml]: problem with class file or dependent class; nested exception is java.lang.UnsupportedClassVersionError: org/geowebcache/security/RESTfulDefinitionSource : Unsupported major.minor version 51.0 (unable to load class org.geowebcache.security.RESTfulDefinitionSource)
    at org.springframework.beans.factory.support.AbstractBeanFactory.resolveBeanClass(AbstractBeanFactory.java:1265)
    at org.springframework.beans.factory.support.AbstractAutowireCapableBeanFactory.createBean(AbstractAutowireCapableBeanFactory.java:433)

...

Caused by: java.lang.UnsupportedClassVersionError: org/geowebcache/security/RESTfulDefinitionSource : Unsupported major.minor version 51.0 (unable to load class org.geowebcache.security.RESTfulDefinitionSource)
    at org.apache.catalina.loader.WebappClassLoader.findClassInternal(WebappClassLoader.java:2822)
    at org.apache.catalina.loader.WebappClassLoader.findClass(WebappClassLoader.java:1159)


texts:
I am trying to get a list of NDVI values as calculated for each image in an ImageCollection for a given point. This list will then be matched to a list of dates for the point, to produce an array of values. However, when I use aggregate_array it only collects 54 values, and not the 100 in the ImageCollection. Could anyone shed some light on this? The code is below:

var aoi = ee.Geometry.Point([-96.17,18.17])
var start = '2017-01-01'; //start date Sentinel 2
var end = '2018-12-31'; //final date Sentinel 2

function maskOutClouds(img) {
  var mask = ee.Image(0).where(img.select('QA60').gte(1024), 1).not();
  return img.updateMask(mask);
}

function addNDVI(img) {
  return img.addBands(img.normalizedDifference([('B8'), ('B4')]).select([0],['NDVI']));
}

var imgs = ee.ImageCollection('COPERNICUS/S2').filterDate(start, end)
                                              .filterBounds(aoi)
    .map(maskOutClouds)
    .map(addNDVI)

print(imgs)

var imgsDate = imgs.map(function(image){
  return image.set('Date', ee.Date(image.get('system:time_start')).format('yyyy-MM-dd'));
})

print(imgsDate,'imgsDate')
var list = ee.List.sequence(1,imgs.size())

var dateList = ee.List(imgsDate.aggregate_array('Date'))
print(dateList)

var dateNDVI = (imgsDate.select(['NDVI']))
print(dateNDVI,'dateNDVI')

var NDVIprop = dateNDVI.map(function(image){
  return image.set('NDVI', (image.reduceRegion(ee.Reducer.first(), aoi, 30)).get('NDVI'));
})

print(NDVIprop,'NDVIprop')

var NDVIvalue = function(ndvi){
  return ndvi.get('NDVI')
}

var NDVIdateSel = NDVIprop.filter(ee.Filter.inList('Date',dateList))
                          .filterBounds(aoi)
print(NDVIdateSel,'NDVIdateSel')

var NDVIdateSelAdd = NDVIdateSel.map(function(ndvi){

  return ndvi.set('NDVI',(ndvi.get('NDVI')))

})

print(NDVIdateSelAdd,'NDVIdateSelAdd')

var NDVIlist = NDVIdateSel.aggregate_array('NDVI')
print(NDVIlist)


texts:
I want to connect from QGIS to a PostGIS cluster via SSL.

I'm able to connect via SSL from PgAdmin4 and psql using the same certificates but not from QGIS.

QGIS throws the error:

  SSL error: certificate verify failed

PSQL (works fine):

psql "sslmode=verify-ca sslrootcert=server-ca.pem \
      sslcert=client-cert.pem sslkey=client-key.pem \
      hostaddr=XXXXXX \
      port=5432 \
      user=postgres dbname=postgres"

PGAdmin 4 (works fine):

QGIS configuration (doesn't work):

This is the connection form. I use a PKI stored identify certificate

I create the PKI stored identify certificate as a PKI PEM/DER Certificate Paths. I provide the path of the client's cert (client-cert.pem) and the client's key (client-key.pem).

The root certificate is also there:

Server error log:

"[214]: [1-1] db=[unknown],user=[unknown] LOG:  could not accept SSL connection: TLSV1_ALERT_UNKNOWN_CA"

It looks the root certificate is not sent


texts: Adding a feature count to a label in the legend in ArcMap

texts: In ArcMap, I want the labels in the legend to include the count of that type of feature. For example, if I have a building feature class with a Type field where 900 buildings have a value meaning residential and 40 have the value for industrial, I want the labels to say "Residential - 900" and "Industrial - 40". What is the easiest way of doing this?


texts: I'm using nominatim as a forward/reverse geocoder. My search target is a city. Nominatim only provides the update mechanism for countries. My current solution is to use something like osmium to get diffs from the geofabrik country updates and then apply the diff to the old city extract. Am I doing it in a roundabout way? Are there any more straightforward ways of doing this?


texts: QGIS is closing when I try to load WIKIMaps/BING/GOOGLE maps loading with OpenLayers

texts:
I am using QGIS 2.14 on my Ubuntu 16.04.3 LTS.
My QGIS(2.14) is closing/quitting when I try to load Wikimaps/BING/GOOGLE maps loading with OpenLayers plugin.
How can I rectify this?


texts: How can I turn editable a shapefile that has been previously saved filtered?

texts:
After describing  this issue in Qgis:
Shape layer turns non-editable if saved filtered

(where the answer basically was "if you filter, don't edit"),

it has happened again (filtering  and editing  is soo useful that I forgot not to do it!)and I wonder if there is a way to turn editable a shape layer where the Toogle Editing  tool remains grey.

My only chance is unfiltering and saving the table with another name, deleting the original table and keeping on working with the new one. But this is not a very pro solution, isn't it?

Does enyone know another way around?


texts: Loading shapefiles into ArcSDE geodatabase using web application?

texts:
We are developing a web GIS using Flexviewer and we want to allow users to upload a shp file and insert it into a ArcSDE geodatabase.  I think there are two steps at least.

upload the shp file (zip) to a directory on the server
call some function to load it to the ArcSDE geodatabase

How can I do this using ArcObjects or Geoprocessing service (ArcPy)?


texts: api-ms-win-crt-runtime-1-1-0.dll is missing after QGIS 2.14 installation

texts:
After I install QGIS version 2.14 ESSEN on Windows 7,error message shows,"The program can't start because api-ms-win-crt-runtime-[1-1-0.dll is missing from your computer.Try reinstalling the program to fix this problem."
I have tried reinstalling the program inumerable number of times but it does not work.Whereas QGIS 2.18 Las Palmas version is working very well.
Kindly help solving this issue.


texts: Can ArcHydro tools be used on a network in an SDE enterprise geodatabase?

texts: Most of the documentation for ArcHydro tools discusses their use on data in ESRI's personal and file geodatabases. Is it possible to use ArcHydro tools to set network flow direction on a networked data set in an 'enterprise geodatabase' (SDE on Oracle)? Both the client and server are at version 9.3.


texts:
How to make a high-resolution map in ArcGIS?

I have a map scale on 1:1.000.000, I have a label in the annotations scale (1:2.500) and set data frame scale on 1:1.000.

I tried to export it on resolution 1.000 dpi (jpeg/pdf/png) but when I zoomed the map it's still blurry, I can't see my labels on it.

I want to make a map similar to this

Does anyone know how to make it?


texts: How do you generate an animation of Sentinel 5P data over a specific time period and area of interest in GEE?

texts: I looked at the examples given in the GEE website for other satellites but they do not work for  Sentinel 5P.


texts: How to calculate stream order from vector data in ArcGIS 10

texts: I know that NVision Solutions has a tool for calculating stream order from vector hydrography data for ArcGIS 9.3.1. Is there something similar that works with ArcGIS 10?


texts:
I am attempting to merge 153 geotifs into one large geotif of the Portland area. This is what my rasters looked like before the merge:
Here is a zoomed in image:

This is what the raster looks like after the merge:
Here is a zoomed image after the merge:

Here is the console call:

gdal_merge.bat -ot Byte -of GTiff -o D:/Desktop/483/startingover_merge.tif --optfile C:/Users/Travis/AppData/Local/Temp/processing_58647a944689450d886257b54252fa5d\mergeInputFiles.txt


texts: ArcMap definition query for a string contained in a longer string not working

texts:
I am trying to query records that contain a specific string within a longer string of an attribute table but it is not working.  There are a few extra details to the situation.

DQ is going onto a polygon layer "Parcels"
"Parcels" is joined to a non-spatial table coming out of an Access DB (MDB)
DQ is applied to values within the joined MDB fields

Syntax for DQ is:

Roll Detail Extract_1.FIX-LGL LIKE '%Section 80%'

The hope is that the DQ will leave only polygons where the term "Section 80" is in the value for the field FIX-LGL.

I have tested the query on an un-joined version of the "Parcels" layer, using an already existing field on a value I know to exist, and it worked. I also tested the query on a joined version of "Parcels" on the same existing field and it still worked, but when i switch the query syntax to the joined attribute and different value, it finds no records.

Is there something with the join table name or field name that is confusing arcmap, ie. spaces in file name or the "-" in the field name?


texts: Extract IDs from polygons intersecting with other polygons

texts:
I'm currently building a model where I have to extract all IDs from the red outline polygons containing the green polygon.

In the provided example I would like the green highlighted polygon to get a row containing the 2 IDs of the red outline polygons.


texts: Implementation of routing with multiple stops using Pgrouting

texts:
I have a network and pgrouting works well on it. I want to route from a start point and car stops at set of points as {point 1, point2, point3 etc}. there is not any order in destination points. well one way is as below:

find nearest point to start point and call it p1
route from start point to p1
find nearest point to p1 and call p2
route from p1 to p2
etc

well I could do this but results are not optimize.

Is there any better solution for this problem?


texts: Understanding difference between Extract Values to Points and Sample tool of ArcGIS Pro?

texts:
Both are in the Spatial Analyst tool kit and report values of rasters based on a given point feature class. I have read both of their web pages on the ArcGIS website and I am still confused about the difference!

Extract Values to Points:

  Extracts the cell values of a raster based on a set of point features
  and records the values in the attribute table of an output feature
  class.

Sample:

  Creates a table that shows the values of cells from a raster, or set
  of rasters, for defined locations. The locations are defined by raster
  cells or by a set of points.

  The input rasters can be two-dimensional or multidimensional. The
  structure of the output table changes when the input rasters are
  multidimensional.


texts: World file in georeferenced file does not work in QGIS?

texts:
In QGIS 2.18.3 on Windows 7, I have geotiffs that were created by georeferencing rasters in ArcGIS 10.2. These rasters do not project correctly in QGIS.

Opening images georeferenced (.aux.xml) in ArcMap in QGIS? refers to the world file associated with each raster as being necessary for proper projection. The world files exist for each raster. However, the created files end with x, as in ".tfwx" and the question I referred to said that x at the end is approximate and that could cause a problem for QGIS.

I closed QGIS, edited the .tfwx to .tfw, and opened the file in a new instance of QGIS, but the raster still doesn't project properly. I tried warping the file and that didn't fix it. (I have several rasters and would rather not reproject them anyway.)

In Properties>Metadata>Origin the values for one of the rasters is 0,601, which doesn't seem correct. They all project fine in ArcGIS. How do I get the files to project correctly in QGIS?

Edit per request for tfwx and aux.xml:

Here is the tfw(x):

1.3609845791669659
0.013698921791920324
0.013698921791776822
-1.3609845791682424
633950.54191826226
3686572.8424918195

aux.xml:

<PolynomialOrder>1</PolynomialOrder>
<SpatialReference xsi:type="typens:ProjectedCoordinateSystem">

PROJCS["NAD_1983_UTM_Zone_14",
GEOGCS["GCS_North_American_1983",
DATUM["D_North_American_1983",
SPHEROID["GRS_1980",6378137.0,298.257222101]],
PRIMEM["Greenwich",0.0],
UNIT["Degree",0.0174532925199433]],
PROJECTION["Transverse_Mercator"],
PARAMETER["False_Easting",500000.0],
PARAMETER["False_Northing",0.0],
PARAMETER["Central_Meridian",-99.0],
PARAMETER["Scale_Factor",0.9996],
PARAMETER["Latitude_Of_Origin",0.0],
UNIT["Meter",1.0],
AUTHORITY["EPSG",26914]]
    <XOrigin>-5120900</XOrigin>
    <YOrigin>-9998100</YOrigin>
    <XYScale>450445547.3910538</XYScale>
    <ZOrigin>-100000</ZOrigin>
    <ZScale>10000</ZScale>
    <MOrigin>-100000</MOrigin>
    <MScale>10000</MScale>
    <XYTolerance>0.001</XYTolerance>
    <ZTolerance>0.001</ZTolerance>
    <MTolerance>0.001</MTolerance>
    <HighPrecision>true</HighPrecision>
    <WKID>26914</WKID>
    <LatestWKID>26914</LatestWKID>
  </SpatialReference>
  <CoeffX xsi:type="typens:ArrayOfDouble">
    <Double>633958.08762850869</Double>
    <Double>1.3609845791669402</Double>
    <Double>-0.013698921791817191</Double>
  </CoeffX>
  <CoeffY xsi:type="typens:ArrayOfDouble">
    <Double>3685755.5644025686</Double>
    <Double>0.013698921791817191</Double>
    <Double>1.3609845791669402</Double>
  </CoeffY>
  <InverseCoeffX xsi:type="typens:ArrayOfDouble">
    <Double>-493017.26483040652</Double>
    <Double>0.73468774906927536</Double>
    <Double>0.0073949625660466234</Double>
  </InverseCoeffX>
  <InverseCoeffY xsi:type="typens:ArrayOfDouble">
    <Double>-2703191.3629040248</Double>
    <Double>-0.0073949625660466077</Double>
    <Double>0.73468774906927548</Double>
  </InverseCoeffY>

Sorry this edit took so long. It took me a while to figure out how to get the XML to show up without being garbled.


texts:
Is it possible to generate a boustrophedon pattern which fill a vector Region with GRASS?

I have a vector region in GRASS. And I would like to generate a vector (a polyline) being the bostrophedon pattern which fill this region.

In the picture bellow, I have the "gray" area, and I would like to generate the red polyline automatically.

Is it possible with GRASS or another GIS tool ?

In a perfect world, the function will take the following argumets :

A starting point (which must be anywhere in the area)
A distance between two parallels lines
A direction (a vector which will help to draw the pattern)


texts: Connecting vertices to form a line without using an attribute

texts:
I have a number of multipoint features that form lines (hiking trails). I would like to connect the dots of each multipoint feature to generate line features, but there are no attributes with which to sort the points in the correct order. The trails vary quite a bit in shape and direction so I can't assign an id from West to East for example.

Is there any way to accomplish this other than manually tracing the lines on top of the points? The main tools at my disposal are FME Professional Edition and QGIS.


texts: How to merge all layers in QGIS project using 'merge vector layers' within PyQGIS

texts:
I am trying to use the 'merge vector layers' tool in PyQGIS to merge all of the layers in my QGIS project into one shapefile, as part of a larger script I am developing. I tried this script, but I get an error related to the input layers. I am using QGIS 3.0.0.

processing.algorithmHelp('native:mergevectorlayers')
input = iface.mapCanvas().layers()
processing.run('native:mergevectorlayers',{'LAYERS':input,'CRS':4326,'OUTPUT':'C:\Users\hp\Desktop\COOPI\GIS_Data\CWW_Predictive_Targeting\GEOREF_HH_GPS_Data\Shapefiles\Merged\TA_Mponda_GVH_Bolera_malunga.shp'})

This is the error I get:

Traceback (most recent call last):
  File "C:\OSGEO4~1\apps\Python36\lib\code.py", line 64, in runsource
    code = self.compile(source, filename, symbol)
  File "C:\OSGEO4~1\apps\Python36\lib\codeop.py", line 168, in __call__
    return _maybe_compile(self.compiler, source, filename, symbol)
  File "C:\OSGEO4~1\apps\Python36\lib\codeop.py", line 99, in _maybe_compile
    raise err1
  File "C:\OSGEO4~1\apps\Python36\lib\codeop.py", line 87, in _maybe_compile
    code1 = compiler(source + "\n", filename, symbol)
  File "C:\OSGEO4~1\apps\Python36\lib\codeop.py", line 133, in __call__
    codeob = compile(source, filename, symbol, self.flags, 1)
  File "<input>", line 1
SyntaxError: (unicode error) 'unicodeescape' codec can't decode bytes in position 2-3: truncated \UXXXXXXXX escape


texts:
Is there a way to convert a shapefile into graph with nodes and edges(except networkx). I'm using qgis and I want to convert a road vector layer into a graph, so that I can use it to apply shortest path algorithms (astar).
Please help me to find a solution.


texts: I am embedding GeoServer's web UI in an IFrame and want to automatically log the user in.  Is it possible to include the login credentials in a query like this http//:host:port/geoserver/login?username=usr&password=pass?


texts: Refresh attributes of a feature type with GeoServer REST API

texts:
Using a PostGIS database, I have created a datastore and a featuretype. When I add or remove columns to/from the table, I need to update the featuretype accordingly. Fortunately, I can use the GeoServer webmin to “reload the feature type”. Unfortunately, this needs a user action.

Is this functionality available through the GeoServer REST API?

(according to this old post, it is not)

The featuretype REST API documentation does not seem to include such thing. This JIRA issue could be what I am looking for, but I don’t know if it has been implemented yet.

If the REST API does not have this functionality yet, I am considering creating a temporary featuretype in order to get the list of attributes. Let’s say I have a workspace ws, a datastore ds and a featuretype ft. For example I would perform the following REST API calls:

POST a new feature type __temp_ft__ to /workspaces/ws/datastores/ds/featuretypes
GET the featuretype from /workspaces/ws/datastores/ds/featuretypes/__temp_ft__.xml
extract the <attributes>...</attributes> from the XML file
PUT the attributes to the featuretype that needs to be updated /workspaces/ws/datastores/ds/featuretypes/ft
DELETE the temporary featuretype /workspaces/ws/datastores/ds/featuretypes/__temp_ft__

Do you think this approach would bring an equivalent functionality?

What kind of problems could arise with these steps?


texts: Setting scale of Data Frame in layout using ArcObjects?

texts:
I'm building an application to automate the production of MXDs for printing. I can add the layers from the database just fine, and configure the page layout. Now I want to set a particular scale to show in the layout data frame view, but I'm running into a problem.

Problem: the Data View scale is different than the Layout View scale.

I am successfully setting the scale of the map with:

map = map_doc.Map(0)
map.MapScale = scale_variable

However, that sets the Data View scale to whatever I want it to be, but in the layout the scale is changed, and I have no idea which scale to set in the Data View so I get the scale I want in the layout. I could do an ugly h4x and multiply my scale variable to the ratio of scale change between Data View and Layout, but this ratio changes with each MXD (as it focus on different regions).

Now, IMapFrame also has a scale property, but setting the scale through IMapFrame does nothing. This following code should work, but doesn't:

map = map_doc.Map(0)
page_layout = map_doc.PageLayout
graphics_cont = page_layout.QueryInterface(esriCarto.IGraphicsContainer)
frame_elem = graphics_cont.FindFrame(map)
map_frame = frame_elem.QueryInterface(esriCarto.IMapFrame)

map_frame.MapScale = scale_variable

It's not that it's not being assigned. There is no error message, and if I print the values it gives me the correct values:

>>> scale_variable = 100000
>>> print scale_variable
100000
>>> map_frame.MapScale = scale_variable
>>> print map_frame.MapScale
100000.0

However, this is not reflected in the map or the layout. Thus, I'm somewhat unsure of what to do next.


texts: Labeling polyline vertices with Z value using ArcGIS for Desktop?

texts: I would like to label all polyline vertices of a shape with the z value of each vertex without using export tools like Entity Vertices to Point and without creating new field in the attribute table


texts: GIS - Coordinates / Projections, Google Maps

texts:
I'm currently working with OS data and have just tried to import a table of google maps coordinates/point data into ARC. The data frame (& layer) coordinate system is as follows:

British_National_Grid
WKID: 27700 Authority: EPSG

Projection: Transverse_Mercator
False_Easting: 400000.0
False_Northing: -100000.0
Central_Meridian: -2.0
Scale_Factor: 0.9996012717
Latitude_Of_Origin: 49.0
Linear Unit: Meter (1.0)

Geographic Coordinate System: GCS_OSGB_1936
Angular Unit: Degree (0.0174532925199433)
Prime Meridian: Greenwich (0.0)
Datum: D_OSGB_1936
  Spheroid: Airy_1830
    Semimajor Axis: 6377563.396
    Semiminor Axis: 6356256.909237285
    Inverse Flattening: 299.3249646

I've imported the coordinates from google maps as a table, exported it as a shapefile defining the coordinate system to Spherical Mercator (ESPG 3857). I've then imported the data, transformed it to the above projection, in the hope that the two would align. The area of interest is London, and the points appear far further south than they should be - several lengths of  the UK further south (but approximately at the right longitude).

Can anyone suggest what I might be doing wrong? I'm at a loss after having followed various online tutorials and how-tos. Sorry if this seems obvious - I'm a beginner.

Thanks -


texts: Error in WMS GetMap requests to cascaded ArcGIS WFS service

texts:
I have a remote WFS service from ArcGIS.
http://<my-domain>:6080/arcgis/services/MyFolder/WFSDemoService/MapServer/WFSServer?service=WFS&request=GetCapabilities&version=1.0.0

I create a WFS store in GeoServer and have no problem with publishing layer from this service to GeoServer. I have a problem with showing a WMS of these published layers.
When I call WMS request manually, GeoServer logs is something like these:
My request:
/geoserver/wms?SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&FORMAT=image%2Fpng&TRANSPARENT=true&LAYERS=my_layer&STYLES=SampleWorldCities_cities8367&DATE=1596003527070&TILED=true&WIDTH=256&HEIGHT=256&SRS=EPSG%3A3857&BBOX=-10018754.171394622%2C0%2C-7514065.628545966%2C2504688.5428486555

error log:
ERROR [data.ows] - Failed to execute request
http://<mydomain>/arcgis/services/SampleWorldCities/MapServer/WFSServer?PROPERTYNAME=Shape&FILTER=%3Cogc%3AFilter+xmlns%3Axs%3D%22http%3A%2F%2Fwww.w3.org%2F2001%2FXMLSchema%22+xmlns%3Agml%3D%22http%3A%2F%2Fwww.opengis.net%2Fgml%22+xmlns%3Aogc%3D%22http%3A%2F%2Fwww.opengis.net%2Fogc%22%3E%3Cogc%3ABBOX%3E%3Cogc%3APropertyName%3EShape%3C%2Fogc%3APropertyName%3E%3Cgml%3ABox%3E%3Cgml%3Acoord%3E%3Cgml%3AX%3E-90.17578125%3C%2Fgml%3AX%3E%3Cgml%3AY%3E-0.17578097424708533%3C%2Fgml%3AY%3E%3C%2Fgml%3Acoord%3E%3Cgml%3Acoord%3E%3Cgml%3AX%3E-67.32421875%3C%2Fgml%3AX%3E%3Cgml%3AY%3E22.105998799750562%3C%2Fgml%3AY%3E%3C%2Fgml%3Acoord%3E%3C%2Fgml%3ABox%3E%3C%2Fogc%3ABBOX%3E%3C%2Fogc%3AFilter%3E&TYPENAME=SampleWorldCities%3Acities&REQUEST=GetFeature&RESULTTYPE=RESULTS&OUTPUTFORMAT=GML2&VERSION=1.0.0&SERVICE=WFS
    2020-07-29 11:02:13,778 ERROR [geotools.rendering] - org.geotools.data.wfs.internal.WFSException: Exception parsing server exception report
    java.lang.RuntimeException: org.geotools.data.wfs.internal.WFSException: Exception parsing server exception report

full error log
When I use preview layer module in GeoServer I got these errors:
WARN [wms.map] - Failed to determine CRS axis order, assuming is EN
org.opengis.referencing.NoSuchAuthorityCodeException: No code "CRS:GC,DEF,CRS,EPSG,4326" from authority "Web Map Service CRS" found for object of type "CoordinateReferenceSystem".
        at org.geotools.referencing.factory.AbstractAuthorityFactory.noSuchAuthorityCode(AbstractAuthorityFactory.java:952)
        at org.geotools.referencing.factory.wms.WebCRSFactory.createCoordinateReferenceSystem(WebCRSFactory.java:207)
        at org.geotools.referencing.factory.AuthorityFactoryAdapter.createCoordinateReferenceSystem(AuthorityFactoryAdapter.java:802)
        at org.geotools.referencing.factory.AuthorityFactoryAdapter.createCoordinateReferenceSystem(AuthorityFactoryAdapter.java:802)
        at org.geotools.referencing.factory.FallbackAuthorityFactory.createCoordinateReferenceSystem(FallbackAuthorityFactory.java:640)
        at org.geotools.referencing.factory.AuthorityFactoryAdapter.createCoordinateReferenceSystem(AuthorityFactoryAdapter.java:802)
        at org.geotools.referencing.factory.AuthorityFactoryAdapter.createCoordinateReferenceSystem(AuthorityFactoryAdapter.java:802)
        at org.geotools.referencing.factory.ThreadedAuthorityFactory.createCoordinateReferenceSystem(ThreadedAuthorityFactory.java:731)
        at org.geotools.referencing.DefaultAuthorityFactory.createCoordinateReferenceSystem(DefaultAuthorityFactory.java:179)
        at org.geotools.referencing.CRS.decode(CRS.java:525)
        at org.geotools.referencing.CRS.decode(CRS.java:453)

I have GeoServer 2.11.0 and 2.17 and I test these in both of them I got same error.


texts: Is it possible to append shapefiles based on field values using ogr2ogr?

texts:
We have 5 "Master" shapefiles(A, B, C, D & E). Each will have different table structure. We have 100's of shapefiles which has features which has to go to A,B, C, D & E.

We have a column in all those 100's of shapefiles, like Column Name: Master & it's values are A,B,C,D,E.

Is it possible to run a query using ogr2ogr and filter out all the features with values "A" and append them to "Master table A" ?


texts: Quality estimation of a spline interpolation with barriers

texts:
I interpolated the age of minerals for the Inner Tien Shan mountains (Kyrgyzstan) using the ArcGIS tool "spline with barriers". I use ~ 200 points for my interpolation. The barrier that I use is a fault map (polylines, green in example). For my further interpretations, I would like to estimate the quality of my interpolation. At the moment I am using the point density for this purpose, but this is only correct for regions where no barriers are (see example below).

Is there a solution, where I can calculate the point density with taking my barriers into account or are there maybe different approaches to estimate the quality of my interpolation?

Example:

Here D1,D2 and D3 have all the same point density of 36.5%, whereas my interpolation quality is different in each region.


texts: topology: points closer each others + report

texts:
I have a point feature class structured like this:

CompanyID | PersonalID | Lat | Long |
1         | 1          | ... | ...  |
1         | 2          | ... | ...  |
1         | 3          | ... | ...  |
2         | 1          | ... | ...  |
2         | 2          | ... | ...  |
2         | 3          | ... | ...  |
3         | 1          | ... | ...  |
3         | 2          | ... | ...  |
3         | 3          | ... | ...  |

I would like to check if points of the same company are too much closer each others (50m buffer zone).
So i run this script applyng the topology rule "Contain one point":

import sys, arcpy, math
from arcpy import env

# Set the environment workspace
env.workspace = r"path"

# Create a file geodatabase

# Create a feature dataset using 4326 as spacial reference
farms = r"farmsPointPath"
outFdsPath = r"geodatabasePath"
outFds = "fds_test"
arcpy.CreateFeatureDataset_management (outFdsPath, outFds, 4326)

# Move the feature class to the feature datasets
fdsPath = r"featureDatasetPath"
arcpy.FeatureClassToGeodatabase_conversion(farms, fdsPath)

# Draw a buffer around the farm point feature according to the proximity check we want to do.
bufDistance = "50 METERS"
farmsBufFds = farmsFds + "_buff"
arcpy.Buffer_analysis(farmsFds, farmsBufFds, bufDistance, "", "", 'NONE')

# Create topology
topoName =  fdsPath+ "_topo"
arcpy.CreateTopology_management (fdsPath, topoName, "")

# Add feature class to topology
topology = r"topology_path"
farmsFds = r"Farms_feature_class_path(within the feature dataset)"
arcpy.AddFeatureClassToTopology_management (topology, farmsFds, 1, 1)
arcpy.AddFeatureClassToTopology_management (topology, farmsBufFds, 1, 1)

# Add rule to topology
topoRules = "Contains One Point (Area-Point)"
arcpy.AddRuleToTopology_management(topoName, topoRules, farmsBufFds, "", farmsFds, "") #error with toponame, have to insert manually

# Validate topology
arcpy.ValidateTopology_management (topoName) #always same error with topoName --> manually

# Export topology errors
out_path = r"outPath"
ExportTopologyErrors_management (topology, out_path, "out_basename")

Actually everything is working, but as output(errors), obviously I have all the points that are included in the feature class (regardless of companyID).
I would like to have different outptus, as:

Check the topology only within the same companyID rows and to have an output related to it
Check the topology with all the others companyIDexluding the same companyID and to have another output.

I don't mind to fix the errors, i just need one (or two) big report to higlight the points affected by these errors.


texts:
I want to simulate urban growth through MOLUSCE plugin in QGIS.

I  prepared my rasters  for two years, and a variable map by using vector format of roads and distance to road map by Euclidean distance and converted it to raster. When I add the data in molusce window the tabs are not active.


texts: How can I plot a map having on my left Asia and on my right America using Geopandas?

texts: The lowres map comes with the conventional mapamundi (America at the left and Asia at the right).


texts:
My mapbox app consists of a basemap and geojson polygon layer. I want nice labels on top of each polygon, but nothing I have tried works

I've tried adding label point layer as geojson and MBTile, making z-index = 99, but the geojson polygons are always on top.

Is there a sure fire method for forcing polygon labels to top layer?

Here is the map


texts: Looking for a free program for 3D visualization

texts:
Looking for a free program or web software for 3D visualization. I'd like to have an ArcScene like abilities.

No editing needed, just need to see the data in a 3D format. It'd be really helpful for presentations. I used ArcScene a ton in college and liked it, but budget constraints do not allow for it here.

Is there anything out there?


texts: Alternative file format to ESRI Shapefile?

texts:
ESRI Shapefile requires multiple file, such as .shp, .prj and etc.
 and it's not preferable for me.
So I'm looking for substitute to it.

Are there good file formats which meet the follows ?

Results in a single file
Can be loaded by both of QGIS and ArcGIS without any conversions
Keep the same kind of information as the ESRI Shapefile
(eg. Don't lose info on projection)
Don't seem to be deprecated soon
Don't get larger size than ESRI Shapefile


texts: Live re-projection of web services in ArcGIS Desktop

texts:
I am trying to provide an easy way for not-technically-literate members within my organization to access watershed subbasin delineations in ArcMap. My idea is to:

1) Pull this USGS web service into ArcMap,

2) Re-symbolize and add a definition query to narrow down to the area of interest,

3) Save to a layer file on our shared drive.

This whole process works fine, but we have a standard to transform our data into our local State Plane projection (EPSG:2876), and this one comes in as Web Mercator (EPSG:3857). I would like to avoid having to host the data ourselves because I'd like for the layer to just be updated whenever USGS refines their data.

Is there a way for me to reproject the service somewhere in this workflow so that the layer file still refers to the web service, but it's projected to our standards?

I am relatively new to working with services, so if there's a different workflow using REST queries that can get me to the same product (a properly projected layer file), I'm all ears.


texts: Store multiple rasters in an enterprise database, independently

texts:
I have a large number of rasters with discontinuous data representing management plans for urban green areas which I want to store in an enterprise database independently. I don't want to mosaic the plans as it should be possible to turn on the plans individually.

When just loading the rasters to the database I get the preferred result but I want to put them in a dataset. Raster mosaic and Raster dataset seems to merge the rasters, and I have understood that Raster catalog is an old format.

The data should be stored in the enterprise database wherefore I don't want to use file geodatabases. Is this possible?


texts: Iterate through dates and through columns using Model Builder of ArcGIS 10.x

texts:
I have a point feature class in a geodatabase (approx 11000 features) with 366 columns.  From these, one represents the targeted date and the other 365 represents the temperature of one day of the year. I added a new column named “Temperature” and I would like to fill in this column with the corresponding temperature value of my date field.
 Hard to explain, here's an example:

DateFin   b1_T_M01     b2_T_M0   Temperature
01-01-18    5.6           3.8       5.6
06-01-18    5.6           4.2
02-01-18    6.3           4.6       4.6
20-01-18    6.3           4.6
20-01-18    6.3           4.6
27-01-18    5.9           4.2

From bx_TMx:  the b refers to the day and M to the month. So for example b1_T_M01 refers to the temperature value of day one of the month one (first of January) .

As I said before, instead of having so many dimensions, I would like to have only the temperature field and add the corresponding value to this field DateFin.
As example: for date 27-01-18 I would need the band 27 from the month 01 (b27_TM01).

That said, does anyone have a clue of how to do that with a Model Builder in ArcGIS, in an iterative way (without having to do it one by one per hand).


texts:
I'm pretty new to maps and geographic systems. Maybe you can help me out with this OSM problems. I'm working on web application in WebGL, and I have two major questions:

I would like to know how to correctly query for city boundaries polygon. I've tried to use JSONP and http.get requests to get data form 'http://polygons.openstreetmap.fr/get_geojson.py?id=1582777' but I can't get it because of CORS policy. This server seems to be JSONP disabled. Does any other kind of similar API exists?

To get proper relations of my address I'm querying 'nominatim.openstreetmap.org' multiple times, every time constructing a new URL. Query for address to get city, city to get country, etc. This works but feels really tiresome. How do I query nominatim, knowing street address and getting all info that I need (city, country, continent)?

Can you give me some hints?


texts: How to set zoom on the basis of the road in Google maps API?

texts:
I worked out a similarity Project a Google traffic, I would implement the same Google traffic visualization about the Data Traffic What a depending on the zoom level displays them or hides them.

For example, the first picture with a zoom

as second picture with higher zoom


texts: How to change location of Google Earth Engine control?

texts:
I want to have the layer selection on the left instead of on the right, because on SplitPanel, the left map's layer selector gets hidden unless you wipe all the way to the right. I know that it is possible to just hide all controls and come up with my own controls, but I kind of like the controls google came up with and they also have quite some extra functionality regarding visualization on top of layer selection. For this reason, I'd like just to change the style of the layer selector. I imagine that if I right just the correct thing as a style, for the correct element, I could achieve that. As of now, I don't know what is the element I need to change the style of, nor do I know what style I should apply.

Here's some code to get 2 linked maps:

var splitPanel = ui.SplitPanel({
  firstPanel: ui.Map(),
  secondPanel: ui.Map(),
  wipe: true,
  style: {stretch: 'both'}
});
ui.root.widgets().reset([splitPanel]);

Any suggestions?


texts: Is a United States City considered an Administration Level 3? Why or why not?

texts:
I am trying to get shape files for all of the worlds admin levels, but also their localities. If you are familiar with the Google Places API, you will notice that they have Country, Admin level 1, Admin level 2... all the way to 5, and them some other place types that I presume they define.

Cities are considered localities by Google Places, but according to this list, cities could be an Admin level 3.

I have dowloaded the gadm.org shape files split into 6 layers and imported into PostgreSQL/PostGIS, and now I am viewing with QGIS.

The problem is, there are NO United States cities. Why is this? Is there a reliable place I can go to get GIS data for cities? I have checked out census.gov and I see incorporated places data, which seems to include cities, but the borders for Austin (as an example) look different than what Google considers "Austin".


texts:
I am trying to align various Sentinel 1 image over an area with a similar path and row as the S1 images provided on GEE  are not coregistered. However, the values change of each image when inspected for before and after registration changes. Have provided the code. Can someone provide clarity on this?

I am doing this step of alignment so that I could draw pixel values by Feature in time.

https://code.earthengine.google.com/?scriptPath=users%2Faakashchhabra7489%2FGEE%3ASentinel%201%20Align


texts: Point cloud query is very slow and not using parrallel scan or more memory

texts:
The following PostGIS query is running very slow (> 1 day for about 1 million records):

UPDATE test.table1
    SET area= ST_Area(PC_EnvelopeGeometry(pa))
    WHERE area is null or area <> ST_Area(PC_EnvelopeGeometry(pa))::real

pa contains the point cloud patch (600 points)
I had created following indexes:

spatial index (using Gist on PC_EnvelopeGeometry(pa)), binary index on area.

The PC has Core I7 processor, 64GB Ram with 1TB SSD & 4TB HDD. I've set the following Postgres server configuration:

shared_buffers = 13280MB, temp_buffers = 4096MB, work_mem = 2047MB, max_stack_depth = 3MB,
min_parallel_table_scan_size = 5MB,
min_parallel_index_scan_size = 128kB
effective_cache_size = 18GB,
force_parallel_mode = on

Please advise what can be done to reduce the query processing time.


texts: how to create a new postgis database for usage with QGIS under windows

texts:
After installing postgresql with postgis under windows xp I have 3 databases:

postgis
postgres
template_postgis

now I try to create a new spatially enabled database in pgAdminIII: I select template_postgis in the pgAdminIII object browser and via rightclick in the context-menu I select  new object ->new database.
when I try to use this  database in QGIS with the spit-plugin I get the following error message:

  the chosen database does not have PostGIS installed, but this is
  required for storage of spatial data.

How do I create a working spatial postgis database?
please help a rookie

thanks

Kurt


texts: Intersecting overlapping polylines using ArcGIS for Desktop?

texts:
I have 4 multipart polyline shapefiles (download a sample dataset here).

The geometries are identical, except for the length of segments in each shapefile. They represent 4 variables describing an identical shoreline. From these 4 overlapping polylines, I would like to create a resulting polyline that has segments for each unique set of the 4 input lines.

My approach was to use an intersect. This, however, breaks down when I try to run it on 4 huge feature classes.

My question is whether there is a more efficient way to combine polyline attributes stored in separate layers. I thought a union would be good, but it only works on polygons. Spatial Join seems inefficient, as I would end up with many to many joins for a single segment. Turning the segment end points into a point layer, then splitting the line at those points?

I have access to arcpy, PostGIS, and shapely.

Has anyone had experience doing this before on large datasets? Is there a more efficient way to store these variables?


texts: How to group by a category and sum multiple variables within that category using QGIS?

texts:
I'm using QGIS for some analysis. I have data in this format:
Name is a text variable and categories are all integers

And I want to generate a sum for all categories for each name, ie a table like this:

Is this possible is QGIS. I've been using the Group Stats plugin and I'm able to get the sums one category at a time but my data actually has many more categories than I've displayed here so it would take me a very long time to do.

I've been putting Name into the row box in Group Stats and in the value box I put Sum and CategoryA.

I'm quite new to QGIS so any help you could give would be great!


texts: Red boxes appear where annotation should be

texts:
I'm working on a project where the workflow includes performing routine data updates from clients. They send me their updated FCs, and I truncate and load all of their data into (what should be) an identical FC on our side. (Delete Rows tool -> Append tool) Recently I did an update where there were some discrepant attribute fields between the two FCs:

You will see that in the Field Map, the fields without links are at the bottom, in all caps. Mysteriously, our FC contains these fields, and theirs does not (I'm relatively new to this project, so I have no clue why). Usually I'm told to ignore fields with no link and append anyway, so I do so.

Afterward, I find that some, but not all annotations in the newly appended FC have been replaced with red boxes:

The red boxes should actually appear as this:

I've never seen anything like this before. I figured our FC was possibly corrupt, so I created an entirely new one, and appended all features into that, but still got the same problem.

Does anyone know what causes annotations to instead appear as red boxes? MY intuition is that it has something to do with the discrepant fields, but as far as I know, we don't use those fields for anything.

I am working in ArcGIS 10.3


texts: Imported CSV point data layer invisible over the raster map layer, though each layer shows good separately "zoom to content"

texts:
I am trying to import CSV file (comma delimited text) and get no reported errors but there are no points data visible overlays over raster data map
 1. QGIS 3.6
 2. the point data layer over raster map layer.
 3. I can zoom to extent for each layer properly, each visible separately.
 4. Each layer data is set the same datum GDA94 EPSG 4283
 5. QGIS

I have tried to find an answer here but none seem to fit my problem.


texts: How do I center OpenLayers map over New Mexico

texts:
I'm trying to center an OpenLayers map over New Mexico (NM) and having a hard time to get it working.

Here is my code. I'm trying to use decimal values and not the pixel values OpenLayers uses by default.

var init = function () {
    // create map
     map = new OpenLayers.Map({
        div: "map",
        theme: null,
        projection:"EPSG:4326",
        controls: [
            new OpenLayers.Control.Attribution(),
            new OpenLayers.Control.TouchNavigation({
                dragPanOptions: {
                    enableKinetic: true
                }
            }),
            new OpenLayers.Control.Zoom()
        ],
        layers: [
            new OpenLayers.Layer.OSM()
        ],
        center: new OpenLayers.LonLat(105.9372, 35.6869),
        zoom: 5
    });
};


texts: How to export a png file from a raster file? r.out.png has an error

texts:
I used r.out.png to export a png file, but it produced this error:

(Sun Feb 27 00:26:32 2011)
r.out.png input=raster_bulacanElevation output=png_bulacanElevation.png
rows = 382, cols = 819
ERROR: Not able to open file for [png_bulacanElevation.png]
(Sun Feb 27 00:26:33 2011) Command finished (0 sec)


texts: Getting and using XY centroids in ArcGIS Field Calculator using Python Code Block

texts:
I am writing a script to obtain Grid reference but cannot fathom how to get XY centroids without creating new fields (all examples I have found are based on writing XY centroids to field).

Here's my current code and the east / north on top is what I want to achieve:

def getParcel:

    ## east = SHAPE.CENTROID.X
    ## north = SHAPE.CENTROID.Y

    digits = 8
    # get the 100km-grid
    e100k = math.floor(east/100000)
    n100k = math.floor(north/100000)
    if e100k<0 or e100k>6 or n100k<0 or n100k>12:
        return ''
    # get numerical equivalents
    l1 = (19-n100k) - (19-n100k)%5 + math.floor((e100k+10)/5)
    l2 = (19-n100k)*5%25 + e100k%5
    if l1 > 7:
        l1 += 1
    if l2 > 7:
        l2 += 1
    # change to integer values
    l1 = int(l1)
    l2 = int(l2)
    # calculate alpha code
    letPair = chr(l1 + (ord('A'))) + chr(l2 + (ord('A')))
    # Strip 100km-grid indices from easting & northing, and reduce precision
    # and change to string value
    e = str(int(math.floor(east%100000 / 10)))
    n = str(int(math.floor(north%100000 / 10)))
    # Pad values if less than 4 digit ref given
    if len(e) < 4:
       e = e.zfill(4)
    if len(n) < 4:
        n = n.zfill(4)
    gridFull = letPair + e + n
        return gridFull


texts: What is the format for OSM lat/lon in .osm/.osh files and how to import correctly into PostGIS?

texts:
I've been having issues converting OpenStreetMap geometries to valid points in PostGIS.

MY PostGIS table has the this definition for the geom column:

SELECT AddGeometryColumn('node','geom',900913,'POINT',2);

Using ST_GeomFromText, all points I enter into this table have points out in the middle of nowhere.

I've isolated one geom for investigation (the lat/lon below is from the OpenStreetMap .osm file):

SELECT ST_GeomFromText('POINT(53.3302090 -6.2642905)', 900913);

I've also tried:

SELECT ST_SetSrid(ST_MakePoint(53.3302090,-6.2642905),900913);

The point produced is no nowhere near land (when displayed in QGIS, also using 900913).

However, I have tested the lat/lon above in a website: http://www.latlong.net/Show-Latitude-Longitude.html, and it correctly identified a point in Dublin, Ireland.

I am assuming the lat/lon is in the WKT format?


texts:
I have a raster slope layer and a raster aspect layer. The slope layer is single band grey, that is black to white, with white as the steepest slope. The aspect layer is pseudo colour.

I'd like to make a combined slope + aspect view.

Using the transparency parameters within either layer's properties ends up mixing the two layers, so I get white+colour and black+colour, which gives a sludgy effect.

I think I need to use the slope layer as a mask over the aspect layer, with white having the maximum transparency.

How do I do this?


texts: How to return a map catalog from a WPS process with GeoServer?

texts:
I guess everything is in the title, I would like to be able to select a map and have it returned.
This is because I'd like later on to return a GeoTIFF which would have been generated by the server, but I'm making small steps.

I already have a process that works fine on the GeoServer WPS demo, but when trying with other clients, instead of offering me a list of maps I have hosted to choose from, I am asked for a URL. (I have no idea as to what I should enter as URL, and I don't care since it is not very user-friendy...)
I think that the problem here is that GeoServer "cheats" when using his own demo client, because it knows the hosted maps, which other clients don't. I am therefore looking for the way, either for the process to return a catalog of the available maps (if it works), or another solution you could offer me.


texts: NASA World Wind Pick From A Group Of Icons

texts:
I am using an IconLayer in NASA World Wind (Java SDK) to show positioned items on a map which works fine except for when several are in a group and start to overlay one another. The display of the icons in a group like that is not so much of a problem as I need to allow the end user to easily select the correct icon in that "pile" easily. A workaround would be to allow for a right-click on an icon in that group would generate a context menu that shows a list of the items in that "pile" so they could select the one they want.

Creating a context menu on right-click is easy enough but I do not see how I could create a buffer from that point to do a simple geometry relationship like "contains" in NASA World Wind. I have searched the API a bit and have not seen anything that looks like it has potential so I was hoping that someone here might now of either how to do this in NASA World Wind or a simple way of including or achieving this functionality.


texts: Extract raster values from RASTERSBRICK using XY coordinates falling within each polygon R

texts:
I have a rasterbrick s (10000 raster layers) that actually is a time series raster. Accompanying s is a SpatialPolygonsDataFrame with many sub-polygons which can be identified by names and/or IDs etc. All data sets have the same spatial extent.

Now I have a SpatialPointsDataFrame wherein lat and lon for each site is specified. The sites are also named.

r <- raster(ncol=36, nrow=18)
r[] <- 1:ncell(r)
s <- stack(r, r*2)
cds1 <- rbind(c(-180,-20), c(-160,5), c(-60, 0), c(-160,-60), c(-180,-20))
cds2 <- rbind(c(80,0), c(100,60), c(120,0), c(120,-55), c(80,0))
polys <- SpatialPolygons(list(Polygons(list(Polygon(cds1)), 1),
                              Polygons(list(Polygon(cds2)), 2)))
plot(r)
plot(polys,add=T)

#v <- extract(s, polys, small=TRUE)

Assume site names in SpatialPointsDataFrame can be generated using:

site_names=replicate(nrow(SpatialPointsDataFrame), paste(sample(LETTERS, 5, replace=TRUE), collapse=""))

Then values for a particular xy pair can be obtained using for e.g:

vals <- extract(s, matrix(c(-106.43, 52.1), ncol = 2))

How can I extract raster time series values for each xy point in SpatialPointsDataFrame per sub-polygon. The column names of the extracted values per point should be identified by the names of sites in SpatialPointsDataFrame. Extracted values for each polygon should be saved in a separate dataframe which I can use to form a list object containing data for all sub-polygons then do further analysis on it. So if I have 2 sub-polygons then I will have 2 dataframesbut with different # of sites (columns) depending on how many sites are in each sub-polygon.


texts: Convert SqlGeometry data from SQL directly to Tiles / PNG file

texts:
I would like to know a basic method for 'exporting' SqlGeometry data from SQL Server 2008 to a Tile Server / PNG File generator. I already have a Silverlight Map application which uses a Custom Tile layer to view Tiles. I would like to add another (this generated tiles) transparent layer of PNG files on top of my base layer.

The 'Spatial Results' in SQL Server Management Studio 'picture' needs to be converted to 256x256 tiles made up of PNG files, with a quadkey so it knows where to place it on the map.

Is there a basic 'quick and easy' way to do this? And, as data gets updated in SQL, the tiles needs to be re-generated.

Thanks a lot.

Cheerio


texts: OpenLayers 3 zooming to layer when source is ImageWMS

texts:
When I create a layer using a vector source I can force the map to zoom to the layer using

//stuff including setting up defaultStyle
var vectorSource = new ol.source.GeoJSON({ url: geoJsonUrl, projection: 'EPSG:3857', extractStyles: false });
var vectorLayer = new ol.layer.Vector({ source: vectorSource, style: defaultStyle });
//stuff
vectorLayer.addEventListener("change", function () {
    map.getView().fitExtent(vectorSource.getExtent(), map.getSize());
}

This works well for me so I'm trying to replicate it with my WMS layer. The layer and source are defined as:

//stuff including setting up params and geoserverRoot
var wmsImageSource = new ol.source.ImageWMS({ url: geoserverRoot + 'wms', params: params });
var wmsImageLayer = new ol.layer.Image({ opacity: 0.3, source: wmsImageSource });

An ImageWMS source doesn't have a getExtent() method (unlike vector sources such as a GeoJSON source), but an Image layer does. So I tried to replicate the way I zoom vector sources with:

wmsImageLayer.addEventListener("change", function () {
    alert('change');
    map.getView().fitExtent(wmsImage.getExtent(), map.getSize());
});

However, the alert doesn't trigger. I have no idea why this is the case.

What is the best way to automatically zoom into an image layer produced by GeoServer? Additionally, why does the change event not fire for the image layer?

EDIT: I tried using GetCapabilities as per advice from user30184
My layer is generated from a parametrized PostGIS view and the URL I used was:

  /geoserver/wms?service=WMS&request=GetCapabilities&layers=areas_for_type&viewparams=typecode:AdministrativeWard

This gives an XML response including:

<EX_GeographicBoundingBox>
    <westBoundLongitude>-1.0</westBoundLongitude>
    <eastBoundLongitude>0.0</eastBoundLongitude>
    <southBoundLatitude>-1.0</southBoundLatitude>
    <northBoundLatitude>0.0</northBoundLatitude>
</EX_GeographicBoundingBox>
<BoundingBox CRS="CRS:84" minx="-1.0" miny="-1.0" maxx="0.0" maxy="0.0"/>
<BoundingBox CRS="EPSG:4326" minx="-1.0" miny="-1.0" maxx="0.0" maxy="0.0"/>

so the BB for the passed in parameter is not being calculated.


texts: Understanding syntax for CalculateField_management in ArcPy?

texts:
I need to use CalculateField_management in a Python script, and the calc expression needs to concatenate a runtime variable, a string, and a field (string) value.  I have searched all over for examples from which to learn.  I have tried many things I have found in my research and some of my own ideas, too, but I just can't get the calc expression to evaluate.  I am now completely confused and need some help to dig out of this hole.  So, here's the scoop.

The field to be calculated, DESCRIP, is a string field in a fGDB feature class with one feature.  The runtime variable, enddesg, will be a string.  Essentially, what I want to do is concatenate the runtime value with a (longish) string and prepend that to the existing value in the DESCRIP field, storing the result back in the DESCRIP field.  I have tried many forms of an "all-in-one" string for the expression, but nothing worked.  I came to the conclusion that a pre-logic approach would be better, but I haven't been able to get that to work either.  Below is what I have right now.

expr = "getDescription(str(enddesg))"
codeblock = """def getDescription(endstr):
    exprstr = endstr + ' Touchdown Zone Elevation;duplicated from existing control point;' + !DESCRIP!
    return exprstr"""
arcpy.CalculateField_management(perpptsSX_src, "DESCRIP", expr, 'PYTHON_9.3', codeblock)  #Should be only one feature in the source

(BTW, the enddesg variable is either a string or an integer but is converted to a string as input.)  I get the following error:

<class 'arcgisscripting.ExecuteError'>: Failed to execute. Parameters are not valid.
ERROR 000989: Python syntax error: Parsing error <type 'exceptions.SyntaxError'>: invalid syntax (line 2)

Failed to execute (CalculateField).

The following code variation produces no errors, but I also get nothing entered in the DESCRIP field.

expr = "getDescription(str(enddesg), !DESCRIP!)"
codeblock = """def getDescription(endstr, desc):
    exprstr = endstr + ' Touchdown Zone Elevation;duplicated from existing control point;' + desc
    return exprstr"""

I am now thoroughly confused.


texts: Using ArcPy SearchCursor/InsertCursor to split attributes in fields?

texts:
I download data daily and create a geodatabase table.

Fields: Company, Customers_Active, Customers Total, County, State

Example Data

{CompanyA,CompanyB}, {141,424}, {1000,3421}, Cook County, Illinois
{CompanyA}, {145}, {1000}, Warren County, Illinois
{CompanyA,CompanyB,CompanyC}, {11,421, 30}, {100,900,60}, DuPage County, Illinois

I want to be able to "break up" the Company data into its own record if more than one company is listed. So I want the following. I don't want the brackets, but I can figure out how to remove them later.

CompanyA 141 1000 Cook County Illinois
CompanyB 424 3421 Cook County Illinois
CompanyA 145 1000 Warren County Illinois
CompanyA 11 100 DuPage County Illinois
CompanyB 421 900 DuPage County Illinois
CompanyC 30 60 DuPage County Illinois

with arcpy.da.SearchCursor("national_outage_map_merge",
                       ["Company", "Customers_Active", "Customers_Total", "County", "State"]) as cursor:
for row in cursor:
    if any(',' in code for code in row):
        print (row)

This code prints the records that have more than one company.

Do I also need to incorporate an InsertCursor?

Would it be easier to create a new record in the same table (and delete the original combined record) or just insert all records into a new table?


texts: How to show layer when select another layer?

texts:
I have an OpenaLayers/gxp application. I have base and overlay layers groups. I want to show some overlay layers when select base layer. For example when I choose BaseLayer1 I want to show OverlayLayer1 and OverlayLayer2 automatically.

I seek something in gxp but can't find anything interesting in class gxp.plugins.LayerTree.  So I think only OpenLayers can help me. Is it possible to show layers when I choose another layer?


texts: How to fix polygon with ring self intersection in QGIS?

texts:
When I polygonize a raster, I very often have 'pixels' holes/interior ring that touche exterior ring.
Example :

POLYGON((5 0, 10 0, 10 10, 0 10, 0 0, 5 0, 3 3, 5 6, 7 3, 5 0))

It's defined as an invalid geometry in the OGC geometry specifications.

In QGIS, how to fix this ?


texts: Esri ArcGIS server: Errors rendering annotation features

texts:
We have a ArcGIS Server 10 site based in Map Services that has some annotation feature classes. Most of the time the annotation text displays correctly. However, sometimes the rendering of annotation features 'goes crazy'.

![Bad Annotations][1]

Our annotation features are plain text, usually numerical or alphanumerical codes with 1 to 4 characters, using font Arial. When the rendering 'goes crazy' the annotations are rendering instead with:
•   A different font (ex. Wingdings), or
•   A different text angle, or
•   A different size, or
•   with underline or strikeout, or
•   with the characters overlapping each other's

These errors occur randomly and can affect all the features displayed from a feature class or not. The way each feature is displayed just keep changing as we zoom and pan. In attachment I'm sending two images that illustrates the phenomena.

Have you ever experienced the same problem? Do you have an idea what could be the source of the problem?

The client side of our site is developed with the Flex API 2.4. The Map Services are not cached as the data used is subject to updates every night. Our data is an ArcSDE 9.3.1 geodatabase, based on Oracle 10.

Your feedback will be welcome.

Regards,
Paula Lopes

Altran.com
Belgium


texts: Error when updating a raster file with a numpy array using rasterio.write

texts:
I'm iterating through thousands of polygons. Each polygon clips a tif raster producing a numpy array. Then the numpy array values are manipulated and the tif raster is updated with those new values. The script works for a few thousand polygons but at some point it throws an error.

import os, sys, datetime, time
import geopandas as gpd
import rasterio
from rasterio.mask import mask
from fiona.crs import from_epsg
import numpy as np
import pycrs
from rasterio.windows import Window
from shutil import copy
import petl as etl
import shapely.wkt

alr_path = os.path.join(working_folder, 'INPUT', 'NI_COMBINED_ALR.tif')
alr = rasterio.open(alr_path)

out_meta = alr.meta.copy()

alr_avg_path = os.path.join(working_folder, 'OUTPUT', 'alr_avg.tif')
copy(alr_path, alr_avg_path)

feature_gdf = gpd.read_file(gdb_polygs, driver='FileGDB', layer=fc_vml_polygs)

for index, row in feature_gdf.iterrows():
    window, out_transform = mask(alr, row, all_touched=True, crop=True)

    x = alr.meta['transform'][2]
    y = alr.meta['transform'][5]

    width = window.shape[2]
    height = window.shape[1]
    col_off = out_transform[2]
    row_off = out_transform[5]

    window_orig = alr.read(1, window=Window(int((col_off - x)/5), int((y - row_off)/5), width, height))

    if np.all([window < 0]):
        avg = -1
    else:
        avg = window[window != -1].mean().item()
    window_avg = (np.where(window!=-1, int(round(avg,0)), window_orig))

    with rasterio.open(alr_avg_path, "r+", **out_meta) as dest:
        dest.write(window_avg, window=Window(int((col_off - x)/5), int((y - row_off)/5), width, height))

This is the error that the last line of the code above throws:

ERROR 1: LZWDecode:Corrupted LZW table at scanline 10808
ERROR 1: TIFFReadEncodedStrip() failed.
ERROR 1: .\OUTPUT\alr_avg.tif, band 1: IReadBlock failed at X offset 0, Y offset 10808: TIFFReadEncodedStrip() failed.
Traceback (most recent call last):
  File "footprint.py", line 277, in <module>
    process_raster(con, pg_table_name, path_tif)
  File "footprint.py", line 206, in process_raster
    dest.write(numpy_array, window=Window(x, y, width, height))
  File "rasterio\_io.pyx", line 1409, in rasterio._io.DatasetWriterBase.write
OSError: Read or write failed

Any ideas about what this error means and how it could be overcome?


texts: Scatterplot, correlation of dependence and r-square between two rasters

texts:
I have two rasters of the same region and of same resolution and I want to plot in a scatter plot their values and compute the best - fit line with r-square. Any idea how can I do this on Arcmap or QGIS or extract the values of the rasters in a file so I can import it into Excel?
Unfortunately I don't use R neither python.


texts: Creating a .rul file for topology in ArcGIS?

texts: How can I create .rul file for topology in ArcGIS?


texts:
I was wondering if it's possible to make a layout, with rounded measured grids, in zoom maps.

I'm using ArcMap and rounding the corners of data frame on layout mode to give a magnifying glass effect over the map.

Here's a picture example of what I've already got.

Thanks.!

Rounded dataframe ArcMap


texts:
I have a QGIS program that is hosted on a SharePoint server that is used by multiple users. I have entered a rule for a point layer that when new features are created default values are added under certain fields. One such is generating the PARCEL_ID of the Field Parcel that the point feature is located in.
To generate this I have set the default value in the fields widget properties to:

if (geomwithin('Field Parcels' , 'PARCEL_ID') is not NULL,geomwithin( 'Field Parcels' , 'PARCEL_ID'),geomintersects('Field Parcels' , 'PARCEL_ID'))

This works for me normally, but for other users and if I use it on a remote assess computer it produces

  syntax error, unexpected '(', expecting COMMA or ')'

All users are using QGIS 2.18.20

Whats going on?


texts: Converting polyline to polygon without 'Feature to Polygon'

texts:
I am attempting to create a script tool to automate some processes at work. Part of this process was to manually cut (using Cut under Editor) a buffer. My approach to this has been to try and split the polygon instead by using the line layer that we were tracing. However you cannot split using a polyline AFAIK, and I would prefer to avoid using 'Feature to Polygon' as Advanced licences are limited at my workplace. For the same reason I am not sure if I can use the 'buffer to one side technique'.

I have been reading up on how to convert polylines to polygons and I have completed part of the process:

# Feature Class to Points
#
# Paul Smith (2012) paul@neoncs.com.au

# Imports
import arcpy
import numpy
from arcpy import env

env.workspace = "U:/ArcGIS"

#Inputs from user parameters
InFc  = "U:/ArcGIS/Default.gdb/clipdefenses"
OutFc = "U:/ArcGIS/Default.gdb/clipdefpoly"

# Spatial reference of input feature class
SR = arcpy.Describe(InFc).spatialReference

# Create NumPy array from input feature class
array = arcpy.da.FeatureClassToNumPyArray(InFc,["SHAPE@XY"],
spatial_reference=SR, explode_to_points=True)

# Check array and Exit if no features found
if array.size == 0:
arcpy.AddError(InFc + " has no features.")

# Create a new points feature class
else:
arcpy.da.NumPyArrayToFeatureClass(array, OutFc, ['SHAPE@XY'], SR)

So I have the points of my polyline, but I am unsure of the next steps to create a polygon out of them.


texts: Rasterising multiple shapefiles into a single raster

texts:
I'm not even sure you can do this, but I want to convert multiple shapefiles into a single single-band raster.

I'm pretty sure its going to be easier to first create a multi-band raster and then use each combination of band values and a lookup table of some sort to assign a pixel a single value, producing a single-band raster. But I'm stuck on creating the initial multi-band raster.

Lets say I have 4 shapefiles (shp1, shp2, shp3, shp4).
The attribute table in each shapefile has a shp_V column that has an integer ranging from 5-8.
I want the raster to have 4 bands, one for each shapefile that equals the shp_V value of one of the shapefiles.

All examples I can find with rasterio or gdal only rasterises one shapefile, which I know how to do.

Currently I have this, but it obviously only converts one shapefile to one single-band raster:

attribute = 'shp_V'
shp1 = "Shapefiles/shp1.shp"
shp2 = "Shapefiles/shp2.shp"
shp3 = "Shapefiles/shp3.shp"
shp4 = "Shapefiles/shp4.shp"
out_tif = "Rasters/raster1234.tif"

cmd = 'gdal_rasterize -a {} -co COMPRESS=DEFLATE
       -co BIGTIFF=YES -tr 0.3 0.3 {} {}'.format(attribute, shp1, out_tif)
print(cmd)
os.system(cmd)

The shapefiles and resulting rasters can get quite large, so any solutions which are a bit easier on processing power would also be a benefit.


texts: How can I decrease size of geodatabase "Nominatim"?

texts:
I downloaded Nominatim http://wiki.openstreetmap.org/wiki/Nominatim/Installation. The database size is more than 400 GB. The speed of geocoding is very slow. I want to increase it.

For example, I need to geocode only to the level of the village or town. Can I remove some not usefull data from planet.osm.pbf or from geodatabase "Nominatim"? I do not need streets, suburbs, buildings, etc. The structure of geodatabase "Nominatim" is quite difficult, so I think that the best way is to change osm.pbf file.

How can I remove data from geodatabase "Nominatim" or from osm.pbf file? Is this the right way?


texts: I have multiple shape files in a folder. I want to assign projection to all the files inside the folder at one go. In batch file i am able to assign projection to all the shapefiles, but i want to select the assign projection to that folder so that all the files get projected.


texts: Densify and smooth out categorical raster with diagonal neighbors

texts:
I'm classifying vegetation types in QGIS (3.10). I have a pretty coarse raster that I want to smooth out. The standard approach seems to be to polygonize and then apply various smoothing algorithms (e.g. Chaiken or sliding_averaging in Grass' v.generalize). This works fine for rectangular raster polygons, but doesn't smooth out thin diagonal features well. For instance, if the raster is

then the output has 3 islands of A's and 4 of C's. Technically, they may be one polygon pinched where the A's (resp C's) meet at the diagonal, but they appear as islands, not truly connected. Smoothing keeps sharp angles at the diagonal touchpoints rather than "building a bridge".

I'm thinking of a solution which involves building these diagonal bridges in a densified raster (for computational efficiency) prior to polygonizing. The output would be something like

It's a 4x4 nearest-neighbor densification of the raster, with then 2:2 and 2:1:1 diagonal touchpoints identified (marked with outline) and smoothed out by changing the cells in bold. Note the categories need to be prioritized, since in the example above at the centre-right, A "beats out" C in terms of building a bridge. This densified raster can then be polygonized with enough breadcrumbs provided for vector smoothing algorithms (if desired) to further adjust the curves while keeping desired connectedness.

Before I muck around in python or build some fine-tuned combo of QGIS/Grass processing algorithms to implement this, am I reinventing the wheel? Or are there equivalent simpler solutions?

I'd be OK doing the diagonal smoothing and bridge-building after polygonization rather than in raster form, but solutions with e.g. + and - buffering the staircase-shaped polygons seem unnecessarily complicated, especially in then filling gaps and resolving overlaps.


texts:
I'm studying how to check and repair geometry correctly and found the "check geometry" tool which seems to work fine.

However there are also tools called "check validity" and "fix geometry" which seem to do the same. How do they work? Are they doing exactly the same as the "check geometry" tool? Is there a situation when one should use these two tools rather than "check geometry"?

After all, "check geometry" gives me a list of repair strategies to choose from, "fix geometry" just chooses a strategy to fix the geometry.


texts:
Does anyone know how to build a mapnik with an ECW driver on Ubuntu.

There is a successfully installed GDAL 2.2.1 with ECW support.

Install ECW 3.3

wget http://www.kyngchaos.com/files/macosxport/libecwj2-3.3-2006-09-06.zip &&\
    unzip libecwj2-3.3-2006-09-06.zip &&\
    wget http://trac.osgeo.org/gdal/raw-attachment/ticket/3162/libecwj2-3.3-msvc90-fixes.patch &&\
    patch -p0< libecwj2-3.3-msvc90-fixes.patch &&\
    wget http://osgeo-org.1560.x6.nabble.com/attachment/5001530/0/libecwj2-3.3-wcharfix.patch &&\
    wget http://trac.osgeo.org/gdal/raw-attachment/ticket/3366/libecwj2-3.3-NCSPhysicalMemorySize-Linux.patch &&\
    cd libecwj2-3.3/ &&\
    patch -p0< ../libecwj2-3.3-NCSPhysicalMemorySize-Linux.patch &&\
    patch -p1< ../libecwj2-3.3-wcharfix.patch &&\
    ./configure &&\
    make &&\
    make install

Compile and install GDAL

wget -c http://download.osgeo.org/gdal/2.2.1/gdal-2.2.1.tar.gz
...
./configure --with-python --with-spatialite --with-pg --with-curl --with-ecw=/usr/local
make
make install
ldconfig

Install mapnik. Im tried:

https://github.com/mapnik/mapnik/wiki/UbuntuInstallation
And that:

curl -L https://github.com/mapnik/mapnik/releases/download/v3.0.12/mapnik-v3.0.12.tar.bz2 | tar -xj -C /tmp/ && \
cd /tmp/mapnik-v3.0.12 && python scons/scons.py configure && \
cd /tmp/mapnik-v3.0.12 && make JOBS=4 && make install JOBS=4

Result

from mapnik import Gdal
Gdal(file='ololo.ecw', nodata=0, shared=True)
ERROR 4: `ololo.ecw' not recognised as a supported file format.


texts:
I have a DEM of a region (metadata here) that I got from converting 16 las files using pointcloudviz. Now I have that loaded up in QGIS, and I have a shapefile I made from a hiking trip in that same area (metadata here). When I add both layers, the coordinates of the dem ranges in the 1000s (presumably meters) while the shapefile is in lon/lat.

I am somewhat aware of CRS, but changing both layer's CRS and the project CRS to the same thing doesn't really help (but I have no idea what I'm doing here). I've tried toggling OTF projection, and that changes things, but not in a useful way. Step by step instructions would be very helpful, since I know nothing about GIS.

I should note that my dem is in .vrt format and I have a large .vrt.ovr file and a small .vrt file. Idk the difference between these two but I have been loading the .vrt.ovr into qgis

EDIT: I think the issue is that I am not sure what CRS my dem is in. I think I have to figure that out before I can continue


texts: pgRouting: pgr_KSP() does not return shortest path when k > 1

texts:
I am trying to get an arbitrary number of paths between two points in a network.  I am using pgr_KSP — Returns the “K” shortest paths.  My aim is to get the K cheapest routes between two points.

However when I implement it with k = 1 it returns the same result as pgr_dijkstra which is the shortest path, this is the expected and desired result.

When I use it with k=2, I would have expected it to return a paths the included the two cheapest one of which should be the path returned when k=1.  The screen shot below shows what is returned for cases of k = 1, 2, 3
The version of pgRouting used is 2.6
The costs returned are

k = 1 - 88899.77

k = 2 - 91758.40:  91980.81

k = 3 - 93202.83: 93922.34: 94149.52

What can I try to get the result where if I used k =3 it would return the paths I am getting from the case where k = 1 and k= 2

edit: What looks like you get for a given k is
a value of k you get the Xth to Yth shortest routes where

x = ((k-1)*k/2) +1

y = ((k+1)*k/2)

So for the example where k = 6

X = 16

Y = 21.

The Code I have to generate these paths is:

DROP TABLE if EXISTS ksp_results;

SELECT pk.seq, pk.path_id, pk.path_seq, pk.node, pk.edge, cost, agg_cost,e.geom
INTO ksp_results
FROM
  pgr_KSP(
    'SELECT gid as id, source, target, length as cost, length as reverse_cost
    FROM edges',
    1,
    37,
    1,  --This is k value

    directed:=true

  ) AS pk
    LEFT JOIN
  edges as e
  ON pk.edge = e.gid
ORDER BY pk.seq

--Now use this to build paths.
DROP TABLE if EXISTS ksp_paths;
CREATE TABLE ksp_paths AS (
  SELECT --cost.seq,
  cost.path_id,
  g.no_links,
  ST_Length(g.geom) as length,
  cost.agg_cost AS cost,
  g.geom
  FROM
    (
      SELECT pk.seq, pk.path_id, pk.path_seq ,agg_cost
      FROM ksp_results AS pk
      WHERE pk.edge = -1
    ) AS cost

    INNER JOIN
    (
      SELECT
      pk.path_id,
      MAX(pk.path_seq) AS no_links,
      ST_Union(geom) AS geom
      FROM ksp_results AS pk
      WHERE pk.edge != -1
      GROUP BY pk.path_id
    ) AS g
    ON cost.path_id = g.path_id
)


texts: get path from where QGIS is loading the plugin dll

texts:
In my QGIS plugin I want to get the path of my plugin. As yet I done this so:

QDir appDir=QDir(QCoreApplication::applicationDirPath());
appDir.cdUp();
QString sPluginFile=appDir.path()+"/apps/qgis/plugins/plugin.dll";

Since QGIS LTR sometimes the plugin.dll is located at

  /apps/qgis-ltr/plugins/plugin.dll

so my code above is not working anymore. Has anobody an idea how to get the path from where QGIS my plugin has loaded?


texts: error of randomforest classification with Landsat data

texts:
I want perform randomforest classification on Landsat data. Here is some code I found from a paper. But when I run into the "classification" part, there is an error said "Error in eval(expr, envir,enclos): object not found t.1". How to solve this error?

###Invoking the R packages
library(sp)
library(maptools)
library(randomForest)
library(raster)
library(rgdal)
###Reading the .txt files of each band and the training data
train = scan('E:/toronto/2009/txt/roi.txt')
b1 = scan('E:/toronto/2009/txt/2009b1.txt')
b2 = scan('E:/toronto/2009/txt/2009b2.txt')
b3 = scan('E:/toronto/2009/txt/2009b3.txt')
b4 = scan('E:/toronto/2009/txt/2009b4.txt')
b5 = scan('E:/toronto/2009/txt/2009b5.txt')
b7 = scan('E:/toronto/2009/txt/2009b7.txt')
###Store all layers into one dataframe
data_all = data.frame(class = train, t.1 = b1, t.2 = b2, t.3 = b3, t.4 = b4, t.5 = b5, t.6 = b7)
###Eliminate all data with value of 0
data <- data_all[data_all$class != 0, ]
data$class <- factor(data$class)
###Training using random forest
rf <- randomForest(class~., data = data, ntree = 500, mtry = 3, importance = T, proximity = T)
importance(rf)
print(rf)
varImpPlot(rf)
###Read the image and store the classification data as .tif file
satImage <- stack('E:/toronto/2009/2009.tif')
outImage <- 'E:/toronto/2009/2009_RF.tif'
###Classification
predict(satImage, rf, filename=outImage, progress='text', format='GTiff', datatype='INT1U', type='response', overwrite=TRUE)


texts: How to find lost information about the projection?

texts:
I have problems to find lost projection and hope someone could give a hint. I think there is no way to find the projection automatically so I did try a semi-automatic approach shown below.

My initial data are ESRI Shapefiles (only .DBF, .SHP, .SHX). The only information I have is that these Shapefiles should be somewhere in Germany. So my idea was to collect first all Projections for Germany. I did that using the following filter on this website:

https://epsg.io/?q=Germany+kind%3APROJCRS

There are about 71 projections. I wrote a small application using Openlayers, OpenStreetMap and proj4js and it looks like this:

A click on a Listitem registers a new projection and I am able to tell at a glance if the projection is correct or not. That works for a small set of Data unfortunally. Most of them apprearing somewhere in the water. Here is the JS Method that load a new projection:

function switchProjection(stringProj4JSDefinition) {

        var arrTmp = stringProj4JSDefinition.split('#');

        var epsg = arrTmp[0];
        var projection = arrTmp[1];

        proj4.defs('' + epsg + '', '' + projection + '');

        ol.proj.proj4.register(proj4);

        map.setView(new ol.View({
            projection: '' + epsg + ''
        }));

        var vectorSource;
        map.getLayers().forEach(function (layer) {
            if (layer instanceof ol.layer.Group) {

                layer.getLayers().forEach(function (groupLayer) {

                    if (layer instanceof ol.layer.Vector) {
                        vectorSource = layer.getSource();
                    }

                });
            }
            else if (layer instanceof ol.layer.Vector) {
                vectorSource = layer.getSource();
            } else if (layer instanceof ol.layer.VectorImage) {
                vectorSource = layer.getSource();
            } else if (layer instanceof ol.layer.Image) {
                vectorSource = layer.getSource();
            } else {
                vectorSource = layer.getSource();
            }

        });

        map.getView().fit(vectorSource.getExtent(), map.getSize());
    }

Any ideas what else I could try?

I am new to GIS.

Most of my shapefiles are already projected to Cartesian coordinate system, so maybe that's problem here but I'm usure.


texts: QGIS - Exporting a raster layer without smearing individual rasters (so you can zoom)

texts:
I've been working with a raster map in QGIS (a shape file made up from squares) and now I'm trying to export it. I want to create an overview map for Germany that allows some magnification to look at individual rasters of 1km². My problem is that QGIS is sort of smearing the individual squares, smoothing the edges.

It also looks like this after exporting it even on 1200 DPI. If I set the scale of the map to a very small scale this doesn't occur, but obviously I don't want it to occur at all.
How can this be done?
Edit: Export Settings:

I've reset it down to see if that changes anything, but it also happens on 3k DPI.
Edit2: This is how it looks in a pdf:

As you can see now we have sharp edges (still no squares) and the squares are also missaligned.


texts: ogr2ogr: How to get center point of line?

texts:
How to get the center point of a line using ogr2ogr command line?

I am unable to get a proper script to get this work.


texts:
I need a coordinate grid per country, 0.08 degrees (?) apart from each other.
It looks like this:

I would want the coordinates of the intersections.
Instead of a perfect rectangle, I would want it to stop at the border of a specific country. Is there any way I can get an array of coordinates this way?
I have made a Python script myself which just makes the rectangle bigger than the country itself, but this isn't the best way.
def get_coords(top_left, bottom_right):

coord = [top_left[0], top_left[1]]
lat = top_left[0]
lon = top_left[1]
coord = [lat, lon]
coords = []
while True:
    index = 0
    if lat > bottom_right[0]:
        lat = coord[0] - 0.08
    else:
        if lon < bottom_right[1]:
            lon = coord[1] + 0.08
            lat = top_left[0]
        else:
            break

    coord = [lat, lon]
    coords.append(coord)

return coords

Any tips or solutions?


texts: How to iterate multiple raster dataset to calculate raster in ModelBuilder?

texts: I would like to iterate multiple raster images using ArcGIS ModelBuilder. The Raster Images include Modis13Q1 NDVI and pixel reliability. I want to perform ArcGIS raster calculator to compute the multiplication of each image (totally 453 images on automated basis rather than doing one by one) as seen in the images below.


texts: How to sample a 3D raster along a 3D path (or 4D raster along 4D path)?

texts:
I have a 3d raster dataset of atmospheric variables (lat x lon x altitude). I would like to sample this raster dataset along a 3d path.

Any idea how I can achieve this, preferable with qgis or postgis or rasdaman?

Edit: Actually, my data is 4d (time x lat x lon x altitude), but having a 3d solution would already be great.


texts: QGIS 3.16 Rasterize (Vector to Raster) Parameters

texts:
Barely started using QGIS last week so I'm still pretty confused with the nodata and burn parameter in raster conversion tool. I tried reading other people's question in this forum and through YouTube (all are previous versions) but I'm still pulling blanks, all I do get are:

Nodata and burn value shouldn't be the same;
Burn means giving the pixels value;
Nodata kind of reads the pixel value;
I tried doing random combinations of nodata and burn value like is shown in the picture, for example burn = 500 and nodata = 0, burn = 500 and nodata = -999, burn = 1 and nodata 0. Some really random numbers basically. But their differences aren't much different with each other so I got confused on which to use.
So ultimately I'm curious on how others decide what values they put in in their nodata and burn parameter when they're trying to rasterize a vector to get the best results they wanted


texts:
I need some help finding circle shapefiles (polygons) within a database, something like tanks and silos, but they weren't classified right in OSM so looking for attributes won't help (doh!).

Is there anyway to select circles in ArcGIS or QGIS?


texts: Avoiding GeoServer/Java out of heap space error?

texts:
I'm new to Geoserver, Ubuntu and Java, but have downloaded a virtual machine from gisvm.com and am getting up to speed. I got as far as configuring it with some fairly large shapefiles from a project I've worked on previously.

My question is related to a problem which I see if I show the shapefile using the OpenLayers layer preview option. I see an error:

  OpenLayers map preview code="internalError" Rendering process failed.
  java.lang.OutOfMemoryError: Java heap space

Googling has led me to plenty of Java command line options to increase heap space, but I have no idea if this should be applied to an environment variable, in a startup script or as a part of the Geoserver config.

Can you help me to understand what I need to edit to get this working?

I'm also wondering if I should be splitting my shapefile into smaller pieces.


texts: Does the $length function in QGIS field calculator support LinestringZ?

texts: I have a shapefile with a MultiLinestringZ geometry type that I created using the grass v.drape function to attach z-values from a DEM. In order to calculate the 3d length of a line segment, can I use the $length function in the QGIS field calculator, or will this function only report the 2d length?


texts: Changing projection of shapefile in QGIS

texts:
How do I change the projection of my shapefile using QGIS?
The default projection is set to EPSG:3003, I want to change it to WGS84 EPSG:32632. When I do a reprojection or I change projection in the properties mask I don't have the correct result.
On the left, I have the plan coordinates but they are not correct for the EPSG shown on the right:


texts: Accurately interpolating elevation from contour lines in QGIS?

texts:
I’m trying to interpolate elevations from contour lines (~70000 at 1m intervals, spanning ~20km^2) in QGIS.

When I interpolate using Inverse Distance Weighting (IDW) in the Interpolation plugin it uses all contour lines rather than just the adjacent contours?

Using the GDAL Raster plugin (inverse distance to a power) also appears to use all contours and treat the contour lines not as lines but as points?

I also tried Triangular interpolation (TIN), treating the contour lines as structure lines, but it crashes QGIS, probably because of the number and density of contour lines.

For powers below 5, using the GDAL raster plugin, the resulting DEM is very noisy, and the contour points tend to become local maxima or minima

i.e. for grid cells that have a contour line running through them, but do not contain the points that define the contour line, the interpolated elevations are shifted towards the average elevation of all contour data, rather than being assigned the elevation of the contour line running through them; this is because the interpolation is unnecessarily accounting for the elevations of every single contour point (though weighted by inverse distance).

While this problem is reduced when I use IDW with higher powers, relatively flat regions start to look like nearest neighbor interpolations - like steps, rather than like a large region with a small incline.

From other posts, it sounds like the GRASS r.surf.contour algorithm has similar problems.

Using a search radius is better, but it still smooths out features smaller than the search radius, and still allows contours not relevant to a grid cell to impact its interpolated elevation.

Is there a way in QGIS to interpolate using only the contour lines immediately adjacent in all directions?

What’s the best practice for interpolating from contour data in QGIS such that contour data is in fact treated as contours, rather than as one large set of points or lines?


texts: Joining HDF files with multi bands and no attribute tables to a shapefile

texts:
I have a HDF file showing raster data and I want to join this to a shapefile (in ArcGIS) so I can aggregate the values to areas rather than cells.  Unfortunately the HDF has no attribute table so I cannot join. When I try to create one it says it can't be done as isn't a single band file.

Is there a way around it so I can join the data?


texts: Open source approach to classifying and removing LiDAR points from overlapping scans?

texts:
The attached screenshot shows a LiDAR point cloud with significant overlapping scans. This point cloud tile is part of the MN statewide LiDAR dataset which has 25% scan overlap and 1.5 NPS. For many applications, these overlaps pose no problem. However, when calculating LiDAR metrics such as forestry grid metrics, these overlapping scan can pose a major problem by creating biases in the resulting products.

There are several commercial applications that can deal with these overlapping scans such as Esri's Classify LAS Overlap and Lastools lasoverage. Is there an open source approach to detect, classify, and remove overlapping points from multiple flight lines?


texts: Clipping raster file in ArcMap over the extent of polyline shape

texts:
I have a calculated raster file, overlayed on a shapefile with lines (like on figure)

I would like to clip the raster in ArcMap (since it is an interpolation and out of the vector shape it is mostly meaningless) to the extent, where the available data are on vector shape (on figure this is a clearly seen wedge or triangle, where the lines are present only). Of course, this can be done manually drawing appropriate polygon around the vector file. But the question is, whether it is possible to derive automatically a polygon mask from this vector layer for clipping the raster? The standard tools seem to suggest only clipping around a rectangular area, containing whole vector data.


texts: Scripting customization for ArcPad AXF file?

texts:
I have an ArcPad AXF file, created from a geodatabase.  I've opened the AXF in ArcPad Studio and created some a custom edit form.  It's all great.

But if the original geodatabase ever gets destroyed and I have to create a new AXF file, I will lose my custom forms, right?

Is there a way to script all that XML that makes up the custom form?

I really don't want to have go through all that ArcPad Studio GUI editing again.

I can see that it is possible to save the XML of a layer (including the custom edit form XML) as a .apl in ArcPad studio.

But how can I push that saved XML back into a new XML file?

Is there a way to script this?


texts: Updating raster or feature layers being used by ArcGIS Server using Model?

texts:
I have a mxd file that contains a few raster and feature layers that was published by ArcGIS Server as a map service. Meanwhile, a model was created to update the contents of the raster and feature layers in the mxd file. The intention is that all the users across my organisation will be able to see changes every time the model runs.

The problem that I have is the geoprocessing by the model failed due to the error that is the data "May be locked by another application." every time I ran it.

Could anyone please help with offering an idea how to get around this issue?

I also tried to use File Geodatabase to store the layers but ended up with the same error.


texts: I have a point layer of households and a polygon layer of tax lots. I want to assign a nearest tax lot polygon to each household point. Note that as the household points are coded on road, not within a tax lot, ST_Within does not work properly. Do you have any suggestions?


texts:
I know this might be a similar question to many others. I have been prompting multiple questions here and on Google but I can't figure it out. I have a raster and I have uploaded it into PostGIS with the following code:

raster2pgsql -s 28992 -I -M r69hn1.tif -F -t 1000x1250 schema.table | psql -d database -U postgres -p 5432

I am not using -C because it gives an error and I learned in other questions that removing it from the code makes it work; and it does load it in.

Then I am trying to open the data in QGIS using DBManager. There it gives me the message that I have the wrong raster type. The error message is in Dutch and not searchable. I have not found anything similar to my exact problem. The data I am trying to use is downloadable from the following link:
http://geodata.nationaalgeoregister.nl/ahn2/extract/ahn2_05m_ruw/r69hn1.tif.zip

The data is coming directly from the PDOK governmental services and is thus a trustworthy link. It is a tile from the AHN 2 height data.

My question is: Is there something wrong with my way of loading in the raster into PostGIS or with visualising it into QGIS. And further, what is wrong and how do I fix it?


texts: Vectors created in OTB Monteverdi are displaced

texts:
I'm using the Orfeo Toolbox Monteverdi in combination with QGIS in order to produce classified vectors of vegetation cover which I then use in maps in QGIS. (I can't use the QGIS OTB plugin because as far as I can tell it doesn't run on OS X - I'd be pleased to be corrected on that though).

Versions:
Monteverdi 1.10 (edit - just updated to 1.12 and tried that, same result)
QGIS 1.80
OS X 10.6.8

When I load the resulting shapefile into QGIS, it always appears displaced by about 75km to the north east of where it should be. As far as I can tell it's not a datum/projection issue as the CRS is the same for all layers.

My work flow is as follows:

produce georeferenced raster image from scanned jpg in QGIS
load raster in Monteverdi
perform supervised classification using SVM Classification tool
export resulting classified vectors as shapefiles
open vectors in QGIS

The vectors then display at a location different to the original georeferenced image.

I'll post the files of interest if anyone would like to check them out to see if I'm doing something wrong.
http://bit.ly/NCXasm - original GeoTiff raster (270kb zipped)
http://bit.ly/TDFWiS - final vectors produced in OTB and then merged to a single shapefile in QGIS (119kb zipped)

Any assistance or advice would be much appreciated.


texts: Change overview map postion in map from right to left

texts:
I am using OpenLayers 2.12. I have added an OpenLayers OverviewMap control. This overview map tool is located on the right side of the map panel by default. I would like it to be positioned on the left side instead (see image).

How can I achieve this?


texts:
Im trying to create a script that will copy the polygon files from one gdb into a newly created gdb but the files always are generated into the default gdb. ive played around with the "outputs" of CopyFeatures but they always still are generated in the default gdb.

My script is below

import arcpy, os
from arcpy import env
env.overwriteOutput = True
env.workspace = "C:\EsriPress\Python\Data\Exercise06\Results\NM.gdb"
fclist = arcpy.ListFeatureClasses()
arcpy.CreateFileGDB_management("C:\EsriPress\Python\Data\Exercise06\Results", "test.gdb")
outputgdb = "C:\EsriPress\Python\Data\Exercise06\Results\test.gdb"
for fc in fclist:
    fcdesc = arcpy.Describe(fc)
    if fcdesc.shapeType == "Polygon":
        arcpy.CopyFeatures_management(fc,os.path.join(outputgdb, "polygon_"+fc))
    else:
        print fcdesc.basename + " is not a polygon"