texts: Unable to add a polyline between points of equal latitude using the Mapbox Static Images API?

texts:
Whenever I try to add a polyline between points of equal latitude, Mapbox responds with an error (example shown below). If I change the latitude of one of the points by 0.00001 degrees, then it works fine.

Does anyone understand why?

Example:  Try to add a polyline between the points (40, -90) and (40, -91).  Using Google’s Interactive Polyline Encoder Utility, as shown in screenshot below, gives the following encoded polyline:  _ocsF~bidP?~hbE

But when using this encoded polyline in a request, Mapbox responds with an error as described below:
Request:
https://api.mapbox.com/styles/v1/mapbox/streets-v11/static/path-1+e00-2(_ocsF~bidP?~hbE)/-90,40,7,0,0/500x400?access_token=YOUR_MAPBOX_ACCESS_TOKEN
Response:  {"message":"Not Authorized - No Token"}

Now change the latitude of the second point to 40.00001.  Using Google’s Interactive Polyline Encoder Utility gives the following encoded polyline:  _ocsF~bidPA~hbE

And using this encoded polyline in a request returns the image I expect:
Request:
https://api.mapbox.com/styles/v1/mapbox/streets-v11/static/path-1+e00-2(_ocsF~bidPA~hbE)/-90,40,7,0,0/500x400?access_token=YOUR_MAPBOX_ACCESS_TOKEN
Response:


texts: How to load the full content of a vector layer using PyQGIS?

texts:
I am trying to load a geojson file using the following code:

capa = QgsVectorLayer("C:/...file.geojson","capa","ogr")
if not layer.isValid():
  print("Layer failed to load!")

QgsProject.instance().addMapLayer(capa)

The geojson file contents a FeatureCollection, some points feature and a multiline.
But just the poinst are rendered. Cannot load the multiline.

How can I do that stuff?


texts: What is road ahead for ArcGIS Runtime .NET/WPF SDKs?

texts:
Visiting Esri Developer website, I see two different sdks for .NET Developers.

ArcGIS Runtime .NET SDK
ArcGIS Runtime WPF SDK

I need to know the main differences between them.  If the first one (.NET) was introduced after the second one (WPF), is the wpf sdk going to be deprecated in the future?


texts: I am preparing a printable map in leaflet. For this purpose I would like to make a map zoom to a particular scale eg. 1:500 or 1:1000. But Leaflet zooms a map based on the extend and zoom level only. Is there any way the map can be zoomed to the scale as desired?


texts:
I'm using OSM-Data via PostGIS (osm2pgsql). To draw a highway with and without bridges or tunnels, I use the following expressions:

The highway:

"highway" = 'motorway' AND "bridge" IS NOT 'yes' AND "tunnel" IS NOT 'yes'

To draw bridges and tunnels for their own I use that expressions:

"highway" = 'motorway' AND ("bridge"  =  'yes' OR "bridge"  =  'suspension')

"highway" = 'motorway' AND "tunnel"  =  'yes'

So far so good. Was working fine for me for years, but now I kept an eye on my PostGIS log and found that:

Result: 7 (ERROR:  syntax error at or near "'yes'"
        LINE 1: ...h') OR ("highway" = 'footway')) AND ("bridge" IS NOT 'yes'))

It seems, there is something working wrong, but what? Bridges, tunnels and the highways are drawn as they should - I can not see any error, but PostGIS is telling me something else.
Is there a better way to draw OSM-Roads with bridges and tunnels?

ADD: the whole line:

2016-11-08T22:12:43 1   1 cursor states lost.
        SQL: DECLARE qgis_9 BINARY CURSOR FOR SELECT st_asbinary("way",'NDR'),"gid","route"::text,"service"::text,"name"::text,"aeroway"::text,"tunnel"::text,"highway"::text,"bridge"::text,"railway"::text,"ref"::text FROM "public"."planet_osm_line" WHERE ("way" && st_makeenvelope(426908.41326787287835032,5749630.84951520431786776,436049.76744384958874434,5752852.1516018183901906,32632)) AND (((((((("highway" = 'path') OR ("highway" = 'footway')) AND ("bridge" IS NOT 'yes')) AND ("tunnel" IS NOT 'yes')) AND ("bridge" IS NOT 'suspension')) OR ((("highway" = 'path') OR ("highway" = 'footway')) AND ("tunnel" = 'yes'))) OR ((("highway" = 'path') OR ("highway" = 'footway')) AND (("bridge" = 'yes') OR ("bridge" = 'suspension')))))
        Result: 7 (ERROR:  syntax error at or near "'yes'"
        LINE 1: ... OR ("highway" = 'footway')) AND ("bridge" IS NOT 'yes')) AN...
                                                                     ^
        )


texts: Where I can find a good regexp to extract sector, district and area from a UK postcode?

texts: I am looking for official/reliable regexps to extract sector/district/areas from a UK postcode. I have found a few online but I am not sure how good they are, and I have not a way to test them properly. Do you know any of those regex that can be considered valid?


texts: What is the PostGIS equivalent of ArcGIS "Intersect" tool?

texts: I am looking for the PostGIS ST_function that will give me the same results as running the ArcGIS Analyis-> Overlays-> Intersect tool on two polygons.


texts: What's the convention for "upside down" labels?

texts:
I have line features that are to be used for labels. These lines are correctly oriented(rotated) and positioned such that text placed on them should give a good cartographic output.

Compare the below two images and their labels. In particular the word "Waverton" and its orientation in relation to "Avenue".

The left is what all labelling engines seem to produce by default. The right is what it "should" look like according to the Ordnance Survey.

So my main question is - Which one is correct and why? What is the cartographic convention here?

I ask this because out of the 5 labelling engines split between: QGIS (has two), MapInfo and ArcGIS (has two), only Maplex has the option to allow the user to insist the label follow the orientation of the line. The others all assume that upside down text is undesirable and automatically "correct" it.

(Maps: © Crown Copyright and database right 2012. Ordnance Survey 100019520.)


texts: Using ST_Difference and preserving attributes in PostGIS

texts:
What is proper syntax for ST_Difference between tables?

I have two huge polygon datasets from shapefiles, call them 'baselevel' coverage and 'deluxe' coverage. They each contain over 2.4 million polygons, with many polygons describing complex coastal features.

The baselevel dataset has complete coverage over the area of interest but less detailed attributes, while the deluxe dataset has gaps and less extensive spatial coverage but more detailed attributes.

I want to create a third dataset which uses all the deluxe polygons, and uses the baselevel coverage to fill the gaps and extend the coverage.  I've been approaching this as two steps: produce a 'fill' dataset which is just the base-level coverage minus all the deluxe polygons, then union the 'baselevel-fill' dataset with the 'deluxe' dataset.

I can do this fairly efficiently in grass via v.overlay (create the fill layer with operator=NOT, then union with deluxe using operator=AND), after dicing the shapefiles into much smaller tiles.

I'm pasting my final syntax below) based on @PaulRamsey's answer, as well as some other Stack Exchange answers and PostGIS tutorials

GRASS is pretty excellent, and could handle my polygons where a lot of other python, ogr, and qgis approaches.  However the fastest I could complete a similar operation on my (very small) test dataset was 154 seconds.  Code below completed in under 28 seconds.

CREATE TABLE deluxe_plus_baselvl AS
WITH deluxe_cutter AS (
 SELECT ST_Union(d.geom) AS geom, b.gid
 FROM baselevel b JOIN deluxe d ON ST_Intersects(b.geom, d.geom)
 GROUP BY b.gid
 ),
baselevel_cut AS (
 SELECT ST_Difference(b.geom, d.geom) AS geom, b.gid
 FROM baselevel b JOIN deluxe_cutter d ON b.gid = d.gid
 )
SELECT 'baselevel' AS type, b.geom, b.gid
FROM baselevel_cut b
UNION ALL
SELECT 'deluxe' AS type, d.geom, d.gid
FROM deluxe d
UNION ALL
SELECT 'baselevel' AS type, b.geom, b.gid
FROM baselevel b LEFT JOIN deluxe d ON
ST_Intersects(b.geom,d.geom)
where d.gid is null;


texts: Is there a way to add an autoincrementing primary key field to an existing spatialite table/layer?

texts:
I know that I can add a new column to a spatialite table either through the Database Manager  (Edit Table – Add Column) or through the Attribute Table – Field Calculator (create new field) but I don´t know which settings to use there to get a new column populating itself with an autoincrementing primary key. I guess it must be possible...

I also have the Plugins Table Manager and QSpatialite installed.

I´m still a bit confused. Why should the newly created "Test_table2" have a geometry column after I create it with above code. The symbol is just the normal table symbol and the field list looks like this

If I then try the "Add Geometry" code above I do not see a "0" im the output but this:

After this step the symbol changes from table symbol to the points symbol in the browser tab but its still not a valid table with geometry. I cant see the two points from my original layer when I show this table as a layer.

If I try to "create the spatial index" by clicking next to the attention sign and confirm the following window with "Yes" (Do you want to create..." then something happens, the red sign disappears but still I can´t see my points in the map canvas.

If I try your other tipp and use the Recover Geometry though I get the following result and not even the symbol changes from the tables symbol into anything else:

I´d just like to understand whats going on here.

It´s probably got to do something with me being new at using code that I don´t recognize a mistake in dealing with this problem.


texts: attribute between points to tiles in Google earth engine

texts:
I have layer of points that are located in Brazil area in very big area.
I would like to get list of the unique tiles of sentienl2 that those points are relate to.
Is something like this possible in GEE?
now I do it manually, using remotepixel and insert the name of the tile manually:

var dataset = ee.ImageCollection('COPERNICUS/S2')
                  .filterDate('2019-11-01', '2019-11-30')
                  .filterMetadata('MGRS_TILE', 'EQUALS','21KXS')
                  // Pre-filter to get less cloudy granules.
                  //.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 100))
                  .select('B1','B2','B3','B4','B5','B6','B7','B8','QA60')
                  .filterBounds(geometry)
                  .map(maskS2clouds);

My end goal: to get list based on my points-feature collection of the tiles those points are belong to


texts:
I have a GRID raster dataset with no calculated statistics. I have tried building pyramids and calculating statistics with no luck. What am I missing?

EDIT:
For those who don't understand:
1. ArcToolbox-->Data Management
2. Raster-->Build Pyramids and Statistics

I have also tried to build pyramids and calculate statistics separately, from the same toolbox.

SOURCE:
http://www.nceas.ucsb.edu/globalmarine/impacts


texts: How to use GeoSeries.distance to get the right answer

texts:
When I want the distance between two points [(117.454361,38.8459879),(117.459880 ,38.846255)] (longitude,latitude) on the earth, I take the GeoSeries.distance method, but the method does not give me the right answer.

The real distance is about 479 meters, however the method give the result far from it, why?

Code:

import geopandas as gpd
from shapely.geometry import Point

geom=[Point(xy) for xy in zip([117.454361,117.459880],[38.8459879,38.846255])]
gdf=gpd.GeoDataFrame(geometry=geom,crs={'init':'epsg:4326'})
dis=gdf.distance(gdf.shift())
print(dis)


texts: Clustering spatial data based on categorical attribute

texts:
For a project I'm doing, I want to cluster parcels together based on their specific land use category. As an example, a small part of the data looks like this.
Spatial Data with different categories of land use highlighted in different colors

There are many more categories than what you see in this picture. I want to know if their is some kind of smart method to "cluster" these parcels based on the land use attribute. You can see from this picture that there is some spatial dependency and I want to create a new column which contains the number of the cluster a parcel belongs to.
I could probably one-hot-encode the category, but that would create a lot of columns and could therefore be suboptimal. Also I'm thinking of writing some simple code based on the 'DBSCAN algorithm' that checks for similar neighbors.
Do you have any suggestions for existing methods using QGIS to evaluate this problem?


texts: Prevent geometry features being drawn if it's outside of a given layer

texts:
I just started experimenting with OpenLayers 4, and I have the following problem:

I have a map with a default raster(OSM), and on top of it a vector layer for drawing features, that covers only a small area. Now when I start drawing (e.g.) a polygon within this area(orange) I don't want the user be able to pull the blue form(form to be drawn) out of this while editing.
Or just disable drawing geometry features on the raster layer.


texts: Google Geocoding API returns "OVER_QUERY_LIMIT" for some addresses, seemingly at random

texts:
I'm using a Python script to send a list of addresses to the Google Geocoding API and parse the JSON results to give me a text list of coordinates.  My problem is: For some addresses in the list, the API fails to return coordinates and instead returns the status "OVER_QUERY_LIMIT", but then it will continue to process more addresses successfully.  For example, in a list of 200 addresses, I will get "OVER_QUERY_LIMIT" for roughly 40 or 60 of them, in no particular order, and the rest will be fine.  But if I run the script again with the same list of addresses, some of those which failed the first time will return a result with coordinates like they should, with status "OK", but randomly different addresses in the list will fail.  Does anyone have a clue why this would happen?

The script uses a loop to send each address as an API request one at a time, including a pause between each request to satisfy Google's limit on requests per second.  And the total number of requests I'm sending should be well below Google's cap on free requests per day.  So I don't believe I'm really going "over query limit" - also note that I am able to continue getting successful results after the ones that fail.  So what could be making some addresses return this status at random?

Here's the script for anyone who's curious:

import urllib, json
import sys, time
import pprint

# specify input addresses as a list (for simplicity of demonstration):
addresses = [
   'address 1'
   'address 2'
   'address 3 etc.'
    ]
# I removed the real addresses for privacy, but you get the idea

for add in addresses:
    time.sleep(.1)
    prefix = 'https://maps.googleapis.com/maps/api/geocode/json?'
    data = urllib.urlencode({"address" : add})
    url = prefix+data
    gresp = urllib.urlopen(url)
    jresp = json.loads(gresp.read())
    if jresp['status'] == 'OK':
        lat = jresp['results'][0]['geometry']['location']['lat']
        lon = jresp['results'][0]['geometry']['location']['lng']
        print (str(lat)+"; "+str(lon))
    else:
        print jresp['status']

My results typically look something like this (for example):

33.782041; -84.4161232
33.8111927; -84.3640217
OVER_QUERY_LIMIT
33.7980347; -84.3691964
33.7396419; -84.3717205
33.7601601; -84.3955309
OVER_QUERY_LIMIT
33.7794269; -84.3675836
33.760549; -84.387061
33.80361; -84.3941337
33.814729; -84.3910812
OVER_QUERY_LIMIT
OVER_QUERY_LIMIT
33.7977534; -84.4076812
33.811963; -84.395412
OVER_QUERY_LIMIT
OVER_QUERY_LIMIT
OVER_QUERY_LIMIT
33.7682217; -84.3849992
33.8480993; -84.4284378
33.9229844; -84.3982854


texts: Simplify multidimensional list of points into track/linestring

texts:
I have a bunch of points representing the movement of a target:

Track_123 = [Point_1, Point_2, Point_3 ... Point_N]

Where each point is a multidimensional vector:

Point_N = [Lat, Lon, Time, Speed, Bearing, ... Var_X]

What I would like to do is to turn the list of points into a simplified track. This is sitting in a PostgreSQL/PostGIS database. I'd like to use something like ST_Simplify but that's limited to 2D or maybe 3D data. Should I do somekind of simplify in numpy and then store those as a linestring in PostgreSQL?

EDIT: What I'm hoping to achieve:
So the points are not evenly distributed along the linestring and some don't add any meaningful extra information. I'm hoping to simplify the line but maintain the other info in the vector. I am leaning towards using some kind of spline to evenly space the positions and time (Lat, Lon, Time_x) and then do individual interpolation along the other vectors: interp(Time_x, Speed), interp(Time_x, Bearing) etc


texts:
I want to select all tenth in a field! :( I know this way; "Km" = 10 or "Km" = 20 ... and so, but I am pretty sure theres other way.

I am not SQL literate, but with your help in this forums I have managed to do some things, however the most basic are missing in different posts and I would like to find a place with basic (and medium) SQL expressions with practical examples to look for.
I don´t usually post, but this could be of interest for some people outside.

Para los que hablan español, esto es una pregunta muy básica y hasta tonta como dice el titulo, pero quisiera un sitio para encontrar expresiones SQL que son básicas, y que quizás para los que no somos letrados en esas "ciencias" podriamos usar de referencia, recuerdo la ayuda de ESRI tenía unos ejemplos de expresiones muy útiles y ultimamente con el Xtools resolvía.
Mi problema es sencillo:
Seleccionar todas las decimas de un campo de Numero de Kilometros. Se hacerlo seleccionando uno por uno ("Km" = 10 or "Km" = 20, y asi) pero debe haber una forma más sencilla.... que pena preguntar algo tan básico pero tanta búsqueda en foros y demás, todos piensan que ya uno maneja estas cosas básicas.
saludos!


texts:
I am trying to get my map to appear in layout view but when i click on the icon (at the bottom left) the layout view is just blank.  Then when i go back to data view all my layers appear again.

How do I get a layer to translate to layout view?


texts: How to update the source code for apps built with ArcGIS Web AppBuilder widgets?

texts:
What is the procedure for updating a widget's code, in an app built using the ArcGIS Web AppBuilder?

I'm editing a widget which is working correctly. When I reload an app which uses this widget, I'm not seeing the changes. This includes using CTRL-F5 and > History > Clear All.

I can see from the file system that the changes I made in the the source directory C:\arcgis-web-appbuilder-1.0\client\stemapp\widgets.... have not flowed through to the app's directory (C:\arcgis-web-appbuilder-1.0\server\apps\31\widgets...) which strongly suggests that a Node-level process is required (and explains why F5 has no effect).

Deleting the app from the Web AppBuilder main page, and creating a new app, does use the new code. Presumably a new app always fetches the latest source code.

Is there a shortcut method? I'm hoping there might be something in Node.js to tell all apps to update their source code, but I'm a Node.js novice and don't know where to start.

I'm trying nodemon which says "nodemon will watch the files in the directory in which nodemon was started, and if any files change, nodemon will automatically restart your node application."

The startup.bat file launches node in the /server subdirectory, which might explain why the widget files in the /client subdirectory are not being tracked. I can't figure out how to start nodemon one level up, so that both subdirectories are monitored.

(Cross-posted to https://geonet.esri.com/message/450805)


texts: reshape coastline based on higher resolution ocean layer

texts:
I am trying to match one layer representing the land with another representing the ocean, they have the same coordinate system but they don't match because the resolution of the "ocean layer" is higher than the "land layer".

I would like to know if there is some easy way to modify the edges of the "land layer" reshaping them based on the "ocean layer" boundaries. I need to do it on a global scale so I need to find some automated way to do it.

I would like to keep the information related to my "land layer" since it is composed by several regions with related features.
Here is an image (yellow = low res land layer, blue = high res ocean layer):

ArcGIS


texts: Problems with plugin/installer manager in QGIS 2.6.1 Debian

texts:
An error has occured while executing Python code:

Traceback (most recent call last):
  File "/usr/share/qgis/python/pyplugin_installer/installer_data.py", line 460, in xmlDownloaded
    if isCompatible(QGis.QGIS_VERSION, qgisMinimumVersion, qgisMaximumVersion):
  File "/usr/share/qgis/python/pyplugin_installer/version_compare.py", line 188, in isCompatible
    if len(maxVer)TypeError: object of type 'NoneType' has no len()

Python version:
2.7.6 (default, Mar 22 2014, 23:03:41)
[GCC 4.8.2]

When trying to open the plugin manages on Qgis 2.6.1 I get the following error:

QGIS version:
2.6.1-Brighton Brighton, exported

Python path: ['/usr/share/qgis/python/plugins/processing', '/usr/share/qgis/python', u'/home/joseph/.qgis2/python', u'/home/joseph/.qgis2/python/plugins', '/usr/share/qgis/python/plugins', '/usr/lib/python2.7', '/usr/lib/python2.7/plat-x86_64-linux-gnu', '/usr/lib/python2.7/lib-tk', '/usr/lib/python2.7/lib-old', '/usr/lib/python2.7/lib-dynload', '/usr/local/lib/python2.7/dist-packages', '/usr/lib/python2.7/dist-packages', '/usr/lib/python2.7/dist-packages/PILcompat', '/usr/lib/python2.7/dist-packages/gtk-2.0', '/usr/lib/pymodules/python2.7', '/usr/lib/python2.7/dist-packages/ubuntu-sso-client', '/usr/share/qgis/python/plugins/fTools/tools']


texts: How to create a new shapefile from JTS geometries?

texts:
What I use:
Java 1.7, JTS 1.13, geotools.api 14

What I have finished:
Creating correct JTS geometries from OSM relations, ways and nodes.

What I want to achieve:
Create a CSV file which contains the following data:
the_geom, name, postal_code, country, osm_id

"the_geom" should be a JTS polygon geometry or a null geometry.(Complains ESRI Specification)

The CSV file makes it possible to easy fill up data.

After those CSV file is completely created I want to read that file and create a shapefile from it. Those attributes like name, postal_code and so on, should be written into a dbf file related to new shapefile.

Following this I want to import the new shape and dbf file into PostgreSQL.

Next step to accomplish requirements written above:
How can I create a new shape and dbf file from any collection containing JTS polygon geometries? There are many classes involved but none of them is explained clearly in geotools.org like creating file headers and so on.

Something like:

public class Shapefile
{
  public Shapefile(String filename)
  {
      //..creates necessary stuff
  }

  public void loadFromCSV(String filename, String delimiter)
  {
     //use addGeometry for each row in csv ...
  }

  public void addGeometry(Geometry g, int esriType, Map<String, String> attributes)
  {
     //esriType = 5 for polygons and so on....
     //attributes = attributes for dbf file. Like name, postal_code ...
  }

  public void write() throws IOException
  {
     //Writes output to shapefile and dbf file ...
  }
 }

I hope someone can explain it.


texts: How to set up Geoserver GWC with ArcGIS compact cache layer?

texts:
I have a huge problem setting up geoserver's geowebcache to serve arcgis compact cache layer and I am getting a bit confused. I have installed geoserver 2.8.0 which contains geowebcache 1.8.0.

According to official geowebcache manpage (http://geowebcache.org/docs/1.8.0/configuration/layers/arcgistilingschemes.html), serving of arcgis compact cache is not possible.

On the other hand, I have found this thread (http://sourceforge.net/p/geowebcache/mailman/message/31884822/) which says that it is possible to serve this data source. I have also found ArcGIS compact cache support in geowebcache RELEASE_NOTES.txt file on github under version 1.6.x. I can also see some code that should cope with arcgis compact cache in the git directories.

I tried to add following section to the geowebcache.xml:

<arcgisLayer>
      <name>naturalearth</name>
      <tilingScheme>/home/geoserver/geoserver2.8/data_dir/gwc/arccache/conf.xml</tilingScheme>
</arcgisLayer>

But geoserver failed to start. Does anyone successfully serve arcgis compact cache? Can you please share the information about how to configure such a service with geowebcache?

Thanx a lot.


texts:
Any known sources of Vector Basemaps for QGIS, excluding "OpenMapTiles.com" that preset upon the installation of the plugin "Vector Tiles Reader Plugin"?

Really looking for a hillshade vector map and a contour line vector map that I can add to QGIS.


texts: Exporting Maps to PDFs from Multiple mxd's using arcpy

texts:
I am trying to export a pdf for each mxd in a folder. So far for code I have:

import arcpy, os
Workspace = r"FilePath"
arcpy.env.workspace = Workspace
arcpy.env.overwriteOutput=True
for file in arcpy.ListFiles("*.mxd"):
    mxd_path = os.path.join(Workspace,file)
    mxd = arcpy.mapping.MapDocument(mxd_path)
    outpdf = file + ".pdf"
    arcpy.mapping.ExportToPDF(mxd,outpdf)
    mxd.save()

I receive the error:
Runtime error
Traceback (most recent call last):
  File "<string>", line 5, in <module>
  File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\utils.py", line 182, in fn_
    return fn(*args, **kw)
  File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\mapping.py", line 1156, in ExportToPDF
    layout.exportToPDF(*args)
AttributeError: PageLayoutObject: Error in executing ExportToPDF

I am confused as to why the ExportToPDF is failing.


texts: How to set different page orientation per page in QGIS Print Composer?

texts: I wanted to have portrait for page 1 and landscape for page 2. How could I achieve that in QGIS Print Composer?


texts:
I'm  to resampling about 30 raster images to different resolution (2, 3, 4, 5 meters) using the ArcPy. I wrote the code but its giving result for one image.
Here is the code:
import arcpy, os
arcpy.env.workspace = "C:/Users/test/input/Year2020/multispectral"
rasters =arcpy.ListRasters ("*", "TIF")
ab= r"C:\Users/test/input/year2020/multispectral"
ab_raster= os.path.join (output_folder, "a.tif")
arcpy.management.Resample(in_raster=raster, out_raster=ab,cell_size="10 10",resampling_type="CUBIC")


texts:
I am very new to GeoServer and currently trying to render WMS tiles on the top of the map that uses Leaflet. When tile loads, it does not load all the tiles. In other words, if it requires, 4 1024x1024 PNG8 tiles, it only loads two and the rest to vnd, or loads none, or very inconsistent.

I think SQL query and everything is fine and when I check log, it complains about error transforming BBOX. However, I have another layer with same setup (just different SQL query) and followed same as it is including EPSG:3395 and BBOX. I did not add cache for these layers.

What would cause this issue and how can I solve it?


texts: Why is Spatialite delivering "Invalid SRID" error message

texts:
Qspatialite in QGIS 1.8 (and also Spatialite-gui) delivers an "Invalid SRID" error message when I try to load shapefiles with a SRID value of 3814 (Mississippi Transverse Mercatur).  The SRID is within the database file spatial_ref_sys table.

For troubleshooting purposes I tried loading the shapefiles directly into a blank QGIS session.  QGIS builds its own custom CRS proj4 code and correctly reprojects the shapefiles on the fly against WGS84 Orthophoto (Bing Maps), but it doesn't attribute the projection as EPSG 3813, instead it's called "Custom 100003".  I figured it should find EPSG 3814.  Also, using the QGIS plugin, Qspatialite has a dialog that allows the user to scroll through a list of available projections, using both common and EPSG names, with SRID value listed as well.  SRID 3814 is not listed in dialog, but I double checked and it actually is in the database's spatial_ref_sys table.

Thanks for the user who had helped me discover that the SRID is already in the table, and doesn't need to be added. Can someone help identify the problem here?  It's not getting read.

EDIT:
Here is the data you requested regarding my prior question (I couldn't copy it into a comment field)

This shapefile correctly reprojects on the fly to WGS 84. EsriShapefile loaded into QGIS 1.8 (Ubuntu 10.04 -- UbuntuGIS Repositories version): USER:100003 - * Generated CRS (+proj=tmerc +lat_0=32.5 +lon_0=-89.75 +k=0.9998335 +x_0=500000 +y_0=1300000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs)

myEsriShapefiles.prj ------ PROJCS["mstm",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",1300000.0],PARAMETER["Central_Meridian",-89.75],PARAMETER["Scale_Factor",0.9998335],PARAMETER["Latitude_Of_Origin",32.5],UNIT["Meter",1.0]]

compare this to http://spatialreference.org/ref/epsg/3814/ 3814.prj -----

PROJCS["NAD83 / Mississippi TM",GEOGCS["GCS_North_American_1983",DATUM["D_North American Datum 1983",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse Mercator"],PARAMETER["central_meridian",-89.75],PARAMETER["latitude_of_origin",32.5],PARAMETER["scale_factor",0.9998335],PARAMETER["false_easting",500000.0],PARAMETER["false_northing",1300000.0],UNIT["m",1.0]]

comparison of prj files below: 1>Esri.prj 2>3814.prj (the order of the parameters is different)

PROJCS["mstm",GEOGCS PROJCS["NAD83 / Mississippi TM",GEOGCS

["GCS_North_American_1983",DATUM ["GCS_North_American_1983",DATUM

["D_North_American_1983",SPHEROID ["D_North American Datum 1983",SPHEROID

["GRS_1980",6378137.0,298.257222101]], ["GRS_1980",6378137.0,298.257222101]],

PROJECTION["Transverse_Mercator"], PROJECTION["Transverse Mercator"],

PARAMETER["False_Easting",500000.0], PARAMETER["false_easting",500000.0],

PARAMETER["False_Northing",1300000.0], PARAMETER["false_northing",1300000.0],

PARAMETER["Central_Meridian",-89.75], PARAMETER["central_meridian",-89.75],

PARAMETER["Scale_Factor",0.9998335], PARAMETER["scale_factor",0.9998335],

PARAMETER["Latitude_Of_Origin",32.5], PARAMETER["latitude_of_origin",32.5],

UNIT["Meter",1.0]] UNIT["m",1.0]]


texts: Load vector layers from PostGIS with a processing script

texts:
I'm trying to build a processing script to automatize the loading of specific vector layers from postgis db. I wrote the following code lines which worked fine in python console :
tablename = "SUIVI_LH"
user = "avteam"
mdp = "password"

from qgis.core import QgsVectorLayer, QgsDataSourceUri

uri = QgsDataSourceUri()
# set host name, port, database name, username and password
uri.setConnection("172.16.0.116", "5432", "foresale", user, mdp)
# set database schema, table name and geometry column
uri.setDataSource ("infra", tablename, "geom")
# define a name for your layer (tablename by default) and "postgres" as data provider
vlayer=QgsVectorLayer(uri.uri(), tablename, "postgres")
if not vlayer.isValid:
    print ('Couche non valide. Vérifiez la validité des paramètres de connexion du script')
else:
    QgsProject.instance().addMapLayer(vlayer)

However I can't find the right method to make it a processing script. I tried those code lines from here : https://stackoverflow.com/questions/62510556/pyqgis-adding-layer-from-postgis-database-and-work-whith-it-in-qgis-interface
from qgis.core import *
from qgis.core import QgsProject
from PyQt5.QtCore import QFileInfo
from qgis.core import QgsVectorLayer, QgsDataSourceUri
from qgis.utils import *

def run_script(iface):
    uri = QgsDataSourceUri()
    uri.setConnection("172.16.0.116", "5432", "foresale",
    "avteam", "password")
    uri.setDataSource("infra", "LH", "geom")

    layer = QgsVectorLayer(uri.uri(), "LH", "postgres")
    if not layer.isValid():
        print("Layer %s did not load" %layer.name())
    QgsProject.instance().addMapLayer(layer)

    run_script()

I also tried this solution for multiple vectors layers which is even better : Zoom on a PostGIS vector layer in a new project with PyQGIS?
from qgis.core import *
from qgis.core import QgsProject
from PyQt5.QtCore import QFileInfo
from qgis.core import QgsVectorLayer, QgsDataSourceUri
from qgis.utils import *

def loadPgLayer(instance, vl, path, name, data_filter=""):
    """load a postgis vector layer"""
        vl.setDataSource(path, name, "postgres")
        if vl.isValid() and vl.setSubsetString(data_filter):
            instance.addMapLayer(vl)

        canvas = iface.mapCanvas()
        instance = QgsProject.instance()

        uri = QgsDataSourceUri()
        uri_params = dict(aHost="172.16.0.116",
                  aPort="5432",
                  aDatabase="foresale",
                  aUsername="avteam",
                  aPassword="password")
        uri.setConnection(**uri_params)

        vl_LH = QgsVectorLayer()
        vl_SUIVI_LH = QgsVectorLayer()

        name = "LH"
        uri.setDataSource("infra", name, "geom")
        path = uri.uri(False)
        loadPgLayer(instance, vl_LH, path, name)

        name = "SUIVI_LH"
        uri.setDataSource("infra", name, "geom")
        path = uri.uri(False)
        loadPgLayer(instance, vl_SUIVI_LH, path, name)

        canvas.refresh()
        canvas.waitWhileRendering()

Both scripts run but do not produce outputs. Any ideas why ?
Log :
QGIS version: 3.16.6-Hannover
QGIS code revision: bfd36fddc9
Qt version: 5.11.2
GDAL version: 3.1.4
GEOS version: 3.8.1-CAPI-1.13.3
PROJ version: Rel. 6.3.2, May 1st, 2020
Processing algorithm…
Algorithm 'model' starting…
Input parameters:
{ }

Model processed OK. Executed 0 algorithms total in 0 s.
Execution completed in 0.02 seconds
Results:
{}

Loading resulting layers
Algorithm 'model' finished


texts:
I have a mapserver installed (Ubuntu Linux 14.04 LTS) with a world.map file which gives me a proper result when I access it as follows:

http://localhost:19090/wms?map=world&layer=world&mode=map&FORMAT=image

However, when I access it using:

http://localhost:19090/wms?map=world&&SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&BBOX=-90,-180,90,180&CRS=epsg:4326&WIDTH=834&HEIGHT=417&LAYERS=world&STYLES=&FORMAT=image/png&DPI=96&MAP_RESOLUTION=96&FORMAT_OPTIONS=dpi:96&TRANSPARENT=TRUE

then I get the following error:

<ServiceExceptionReport version="1.3.0" xsi:schemaLocation="http://www.opengis.net/ogc http://schemas.opengis.net/wms/1.3.0/exceptions_1_3_0.xsd"><ServiceException>
msProcessProjection(): Projection library error. proj error "major axis or radius = 0 or not given" for "init=epsg:4326:+epsgaxis=ne"
</ServiceException></ServiceExceptionReport>

Mapserver version: 7.0.0

I have already tried out the solutions that Googling this error generally gives me, namely:

Made sure that the epsg file is located in the PROJ_LIB folder
Exporting PROJ_LIB folder using an environment variable
Checking that the definition of epsg:4326 is in the epsg file

Here us the output of running gdalinfo world.tif:

Driver: GTiff/GeoTIFF
Files: world.tif
Size is 1350, 675
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (-180.000000000000000,90.000000000000000)
Pixel Size = (0.266666666666667,-0.266666666666667)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (-180.0000000,  90.0000000) (180d 0' 0.00"W, 90d 0' 0.00"N)
Lower Left  (-180.0000000, -90.0000000) (180d 0' 0.00"W, 90d 0' 0.00"S)
Upper Right ( 180.0000000,  90.0000000) (180d 0' 0.00"E, 90d 0' 0.00"N)
Lower Right ( 180.0000000, -90.0000000) (180d 0' 0.00"E, 90d 0' 0.00"S)
Center      (   0.0000000,   0.0000000) (  0d 0' 0.01"E,  0d 0' 0.01"N)
Band 1 Block=512x512 Type=Byte, ColorInterp=Red
Band 2 Block=512x512 Type=Byte, ColorInterp=Green
Band 3 Block=512x512 Type=Byte, ColorInterp=Blue

The full map file I am using is as follows:

MAP
  NAME   "world"
  CONFIG "MS_ERRORFILE" "/wmserver/data/logs/mapservice_errors.log"
  CONFIG "PROJ_LIB" "/wmserver/deps/proj.4-master/share/proj"
  CONFIG "ON_MISSING_DATA" "LOG"
  DEBUG 5
  EXTENT -180.000 -90.000 180.000 90.000
  SHAPEPATH "/wmserver/data/service_data/data/wms"
  IMAGECOLOR 255 255 255
  SIZE 1350 675

  WEB
    METADATA
      "wms_title" "Georbis WMS Server"
      "wms_author" "VizExperts"
      "wms_onlineresource" "http://localhost:19090/wms?map=world&"
      "wms_enable_request" "*"
      "wms_srs" "EPSG:4326"
      "wms_feature_info_mime_type" "text/html"
      "wms_format" "image/jpg"
    END
  END

  PROJECTION
    "init=epsg:4326"
    "proj=latlong"
    "ellps=WGS84"
    "datum=WGS84"
  END

  LAYER
    NAME world
    TYPE RASTER
    DATA "world.tif"
    METADATA
      "wms_title" "BlueWorld Map"
      "wms_srs" "EPSG:4326"
      "wms_server_version" "1.1.1"
    END
    PROJECTION
      "init=epsg:4326"
      "proj=latlong"
      "ellps=WGS84"
      "datum=WGS84"
    END
  END
END

Mapserver error log:

[Mon Sep 14 04:34:57 2015].834042 msWMSLoadGetMapParams(): enabling non-square pixels.
[Mon Sep 14 04:34:57 2015].834072 msProcessProjection(): Projection library error. proj error "major axis or radius = 0 or not given" for "init=epsg:4326:+epsgaxis=ne"
[Mon Sep 14 04:34:57 2015].834103 mapserv request processing time (msLoadMap not incl.): 0.000s
[Mon Sep 14 04:34:57 2015].834118 msFreeMap(): freeing map at 0x1fa2b60.

What am I doing wrong?


texts: Finding list of bus stops and what is next

texts:
First of all apilogies if this is already discussed. I am a low vision person and there is only so much time I can spend searching for information as it is stressful on the eyes. I spent considerable amount of time before posting here.

As I stopped driving and started taking public transit, I am feeling the need for an app that is customized for my needs. While on the bus, it is hard to either look out and find where I am and if the bus does announce the next intersection, it is most of the times, impossible to hear. So, what I need is an app that tells me what is the next stop as the bus travels along the route. It would also be helpful if it announced the next intersection/street. This would help me in figuring out when to pull the cord to request a stop.

What sdk/api would help in this? I am looking at google maps api or anything that I can use to integrate with the app on iphone. Given that this is for a low vision person, the visual map is less important. What I plan to do is, use the info that I get from the maps and feed it to the iOS accessibility framework to convert to voice. Any help and pointers will be greatly appreciated.


texts: GeoNetwork map viewer not displaying title name in Thai language?

texts:
I am using GeoNetwork 2.10.3 to Add WMS server from GeoServer and am unable to display its title using the Thai language.

I have tested an edit string with the Thai language and it worked.

It is only the map viewer which cannot display it.

Can you help me to fix this problem?


texts: How to display two line traffic road map cartographically?

texts:
I have two line traffic maps. I want to display these maps roads as image below. These images show that when I zoom out from map, roads are showing properly.. And two line road never disappear from map. And when I zoom in map, other roads appear properly and two line. How can I do that this trick with ESRI solutions ?

Example - http://www.arcgis.com/home/webmap/viewer.html?webmap=30c24162643b43c291c8151b24674c7f

I follow ESRI's representation tutorials from this link but I could not get this view.

EDIT 1 -
All I want to do is show 2 lines of road every scale. For viewing every scale width of road may be changeable.


texts: Storing distance between two points in KM

texts:
I am very new to PostgreSQL and PostGIS. I am filtering a collection of messages based on their distance from a point which works great, but I would like to be able to store the distance between the two points in km as well. I've tried a few different things with no luck.

This is the function I am using:

CREATE OR REPLACE FUNCTION public.search_posts_near_user(userid integer, distance_kms double precision)
 RETURNS SETOF geo_posts
 LANGUAGE sql
 STABLE
AS $function$
  SELECT  A.user_id, A.location,
  (SELECT json_agg(row_to_json(B)) FROM post B
   WHERE (
     ST_Distance_SPHERE(B.location::Geometry,A.location::Geometry) /1000) < distance_kms
   ) AS nearby_posts
  FROM user_location A where A.user_id = userid;

$function$


texts: Changes in MOD09A1 quality mask as processed with GDAL

texts: I'm analyzing MOD09A1 spectral data using GDAL/Python.  The quality mask is a uint32 dataset and want to pull out the value of the first bit to mask the spectral data (1 == poor data; 0 == good data).  When I read the HDF directly into a Numpy array without mosaicking the HDFs or applying gdalwarp to reproject/clip the matrices, and build my quality mask from the uint32 dataset (quality_mask = quality_data & 1), there are plenty of values in the matrix at which bit_1 == 0 and bit_1 == 1.  IF, on the other hand, I process the HDF data with GDAL (gdalwarp to reproject, transform to GTiff, clip), suddenly ALL of the bit1 values are equal to zero, i.e. all of the data values are flagged as good quality, which is not the case.  Any thoughts on how/why the uint32 data is getting transformed in this way as I run it through gdalwarp?


texts:
I am searching a plugin for Leaflet, to have an infobox at the border of the map. Please look at my screenshot.

A click on the info button should be open a window. In this window I would like to have some additional informations.

I have not found such a plugin within https://leafletjs.com/plugins.html
Maybe I have overlook a plugin with this function?

Background software: I am use the gis qgis with qgis2web plugin. With this plugin  I can export a fully leaflet website. After export I have the possibilty to change the code and maybe insert a leaflet plugin.


texts: In a recent post, a responder made reference to the Topology Framework .Net.  I'd never heard of it before, having gently used the NetTopolgySuite (and DotSpatial) before with satisfactory results.  Why would one use TF.Net over NetTopologySuite, or vice versa?  Why does TF.Net exist if NetTopologySuite has been around for quite a while?


texts:
Is it possible in CartoDB to have multiple data points for the same timestamp so that when I do things like Torque or other viz it shows up with all the data at that time point? In other words the data would look like:

time lat lon
xxx la1 lo2
xxx la3 lo4


texts: How to export a single feature (line) from multiple layers and export to individual layers?

texts:
I have multiple vector layers representing a time series. Each layer contains many lines each. I want to select a single specific line from each of these these layers and save it to another layer. This line represents the same feature (coastline) in every layer but it's position changes over time.

I can certainly do it 1 by 1 using "export - save selected features as...", but that would take forever. Also, I can't highlight all 40+ layers and use "Select Features by area" tool...

This is what 1 layer looks like:

This is all the layers with the info I want!


texts: How to verify a PostGIS installation on OS X?

texts:
I am trying to set up GeoDjango and installing PostGIS on my mac and when I am trying to

./configure

I am getting this error

checking PostgreSQL version... PostgreSQL 8.4.9
checking libpq-fe.h usability... yes
checking libpq-fe.h presence... yes
checking for libpq-fe.h... yes
checking for PQserverVersion in -lpq... no
configure: error: could not find libpq

Is this alright or should i go ahead and execute

make

??

I have installed GeoServer Suite which has Postgres, should first install Postgres separately as well?

UPDATE:

 PostGIS is now configured for x86_64-apple-darwin10.8.0

 -------------- Compiler Info -------------
  C compiler:           gcc -g -O2
  C++ compiler:         g++ -g -O2

 -------------- Dependencies --------------
  GEOS config:          /usr/local/bin/geos-config
  GEOS version:         3.3.3
  PostgreSQL config:    /usr/local/pgsql/bin/pg_config
  PostgreSQL version:   PostgreSQL 9.1.3
  PROJ4 version:        48
  Libxml2 config:       /usr/bin/xml2-config
  Libxml2 version:      2.7.3
  JSON-C support:       yes
  PostGIS debug level:  0
  Perl:                 /usr/bin/perl

 --------------- Extensions ---------------
  PostGIS Raster:       disabled
  PostGIS Topology:     enabled

 -------- Documentation Generation --------
  xsltproc:             /usr/bin/xsltproc
  xsl style sheets:
  dblatex:
  convert:
  mathml2.dtd:          http://www.w3.org/Math/DTD/mathml2/mathml2.dtd

Does this mean that PostGIS is now installed? Why is the PostGIS Raster disabled?

Also when i follow the tutorial I am getting this error

createdb -E UTF8 template_postgis
createdb: could not connect to database postgres: FATAL:  role "root" does not exist

Then i checked for the command

sudo su - postgres

and it returns 1 for echo $?


texts: How can I label both coordinates in mmQGIS? Everything I try leads to only one of the 2 coordinates (lat or lon) to be labeled. I cannot get both labeled at the same time.


texts: How can I change the CSS style of the GeoServer rendered map produced when I click on "Open Layer "Tab under "Common format"?

texts:
I have tried many times with left no stone unturned, but all efforts came in vain.  I might be of less knowledge in this area.  Then through surfing I found this site as a newbie. I will be very obliged if I get a solution to the below mentioned scenarios.

In GeoServer when I click on "layer preview" button then "OpenLayers" for preview, rendered map aligns to the left. I need to align center and give some CSS style and others to the rendered page as below.

Left side panel for navigational purpose through layers made from PostGIS data.
Middle side is intended map.
Right side is for showing some attribute table for the left panel chosen layer

How can I do this?


texts: Coordinate formatting for swath data - to polygon

texts:
I've got a CSV of Lidar coverage swaths in the below format. Each run or strip is tied together via run_id. I'm assuming L is for left, R is for right and M is for middle.

EDIT: These are coordinates extracted from our flight software's database to an Excel table of the swath boundaries; this is not actual lidar data, just the outlines of what was flown and recorded.

How would I format this for easy import into QGIS, Arc, or Global Mapper as polygons?

    Polygons from Left, Right and middle point coordinates-swath.jpg


texts: Kernel density estimation for frequency data on polygon level

texts:
I would like to make kernel density estimations for a variety of influencing factors on crime, but I do only have frequency data on polygon level, e.g. number of drug markets on county level.
For all kinds of kernel density estimations (KDE) or hotspot analysis, point data are required as input. Do you think I can nonetheless make such analysis and if yes, how?
My workaround would be to use the centroid of each polygon as point input features and set a weighting scheme to the number of observations. Do you think this is a valid approach in terms of usefulness of the models?
For the implementation I can work with any kind of GIS software or programming language (R, python).


texts: I am presently using LAStools to analyse LiDAR data. Can LAStools classify the points into various features such as vegetation, roads, buildings etc? If so, how is it done?


texts: What is the best approach to interacting with map layers after a JTX Job Change event when a new MXD is loaded?

texts:
To further clarify, my question is pertaining to a WMX/JTX job activation from within ArcMap, not Workflow Manager.

I have a JTXJobListener class where with a JobChanged method and that is successfully called.  My problem is with the event order.  I can do whatever I want with the map, but if that WMX/JTX job loads a MXD...I lose those changes (changes are on the original map not the new one that is loaded).  I'm in need of an event that fires after the job type MXD is loaded into ArcMap.

Am I missing an event that I can use?  If not, what are the typical workaround(s) to this issue? Hack it using the document events??


texts: Meaning of number after username in file geodatabase schema lock?

texts:
I'm using the ISchemaLock interface to return schemalock info on featureclasses in a file geodatabase:

// If schemalock count is >=2, then we have a problem, something else has a lock on the FC
                if (slc >= 2)
                {
                    c.msg = "Schemalock count on FC '" + fc + "' in FGDB '" + kvp.Key[0] + "': " + slc.ToString();
                    Console.WriteLine(c.msg);
                    Messaging.Log(c.msg, c.lw);
                    ISchemaLock sl2 = (ISchemaLock)ds;
                    IEnumSchemaLockInfo enumSLI2;
                    sl2.GetCurrentSchemaLocks(out enumSLI2);
                    ISchemaLockInfo sli2 = enumSLI2.Next();
                    while (sli2 != null)
                    {
                        c.msg = "  Schemalock info: " + sli2.TableName + " : " + sli2.UserName + " : " + sli2.SchemaLockType;
                        Console.WriteLine(c.msg);
                        Messaging.Log(c.msg, c.lw);
                        sli2 = enumSLI2.Next();
                    }
                    keysToDelete.Add(kvp.Key);
                    break;
                }

sli2.UserName returns something along the lines of:

HOU-AGSPRD01.10768

HOU-AGSPRD01 is my server name obviously, but can anyone tell me what the "10768" is? I'm thinking it's the PID of the process that has the lock on the particular featureclass. Any ideas?

EDIT: Kirk and Matt are both correct (can I not set 2 answers to correct?):


texts: Reading ushort and uint data using gdal C# Band.ReadRaster method

texts:
I am using C# bindings for gdal in order to get elevation data from GeoTIFF files. The API contains a Band.ReadRaster method. This method supports a lot of data types, but for reading unsigned integers, only contains a method compatible with byte data. However, GeoTIFF files also support data being in ushort or uint format.

How can I use the byte version of ReadRaster to read these formats?

My assumption is to create a byte array that is of length [img width * img height * 2] for ushort and [img width * img height * 4] for uint, call the ReadRaster method, and then convert every two bytes in the array to a ushort (or 4 bytes to a uint).

However, if this is the correct path, I am not sure if I need to worry about byte ordering (little endian vs big endian)?

What would be the method I need to use to convert the bytes to ushort or uint?

Will BitConverter.ToUInt16 and BitConverter.ToUInt32 work?


texts: Is PgRouting DARP Function still being supported?

texts:
I'm asking this, because I need some tool to solve a DARP problem.

After researching I found PgRouting had included this functionality in one of its branches. But after installing it, getting stuck and looking for other people problems similar to mine it seems as if no one is answering DARP Releated questions. (Not just here but in general, PgRouting mailing list included)

Example:
https://gis.stackexchange.com/questions/18999/how-can-i-make-darp-function-of-pgrouting-working

And this one which has a year already hasn't been answered neither:
Help making the DARP example run (I have the same doubts)

So finally my question is this functionality being supported? Is anyone working with it? or it has just being deprecated?


texts: Reducing file size by not storing no-data points

texts:
I have to compute sea wave forecast for large period of time. Land information is of no use here.

Is there any file format that allows to not store no-data/fill values for land area to reduce file size? Chosen format for now is NetCDF, but it can be changed.

Other idea is to store only sea point values in binary file and write or read it using land sea mask. But it would probably have to be some custom data format.

Is it good or common practice to not include no-data points information in file?


texts: How to bind GraphicsLayer to RadGridView?

texts:
I'd like to bind an Esri WPF GraphicsLayer to Telerik's RadGridView control.

Does anyone have any examples showing how to do this?

Examples showing how to do it with Silverlight would also be useful.


texts: why does the North American Datum not actually include north america in it's boundary?

texts:
http://spatialreference.org/ref/epsg/4269/

That link shows the area that the North American Datum 1983 covers and one thing I note is that it doesn't cover North America. It covers all of Asia and Europe and north Africa but not Europe. Why is this?


texts:
How can I set the cursor to busy-mode/in progress from my qgis plugin (cpp)?

I couldn't find anything in the API documentation.


texts: How to remove or delete column fields in a shapefile using Python?


texts:
I'm trying to track down a UK health geography shapefile and have hit a wall. The name of the specific geography that I'm after is 'NHS Area Team'. However I am looking for a historical version as the more recent version of NHS Area Teams have seen lots of mergers due to the shenanigans of NHS reorganisation.

Normally I would go to the ONS geoportal but they have stopped hosting historical data (annoying). I have emailed them about it and they don't seem to have what I'm after.

Here are the areas which I'm trying to obtain (25 in total) from an old ONS article.

I may be able to construct these by dissolving CCG boundaries to match but I'd quite like to have the proper ONS codes for reference too.


texts: I have a file that contains four rows. The first row indicates the longitude, the second row indicates the latitude, the third row indicates the colour of a specific longitude and latitude and the fourth row indicates the transparency of that specific location. I know how to input the first 3 layers and make it work, but I am unable to adjust the transparency of every single location. The data looks, for example, like this:


texts: How to change the symbology on many image files within a group?

texts:
I have a group of raster images from lidar data (so a DEM).  When loaded in, they all have different symbology - so their color scheme and color ramp are all different.  I would like to create a group with images that all have the same symbology, so that I can visually compare topography between the adjacent images.  There are approximately 600 images within this group, so I don't want to change the color scheme for each image individually.

I've tried changing the symbology of one image, and using the "apply symbology to layer" function, but this does not work because it seems my layer or group of images is not considered a raster layer... I also tried first making this layer a raster layer with "make raster layer", but the layer or group that I have is not considered a valid type of input.

I'm using arcmap 10.0.  Any advice or suggestions on how to do this would be most helpful.


texts: How do I post an issue ticket to the GitHub OSGEO/Ggdal repository?

texts:
After finding an issue with the gdal.py script, I would like to report the issue to Gdal. I looked to the GitHub repository and did not see where to post the issue.

The issue


texts: Choropleth raster to polygons with values

texts:
I'm trying to turn a PNG that looks like this to a shapefile with polygons that keep the categories marked by the colors on the map.

I have several maps to process, so manually tracing the polygons is not a valid option.

I tried to do this with arcpy.RasterToPolygon, from here :

import arcpy
from arcpy import env
env.workspace = "C:/data"
arcpy.RasterToPolygon_conversion("zone", "c:/output/zones.shp", "NO_SIMPLIFY",
                                  "VALUE")

In this case "zone" is the input PNG, this layer is already georeferenced to its actual location. zones.shp is the filename where I want to store the new features.

I get the following error:

Runtime error  Traceback (most recent call last):
File "<string>", line 5, in <module>
File "c:\program files (x86)\arcgis\desktop10.4\arcpy\arcpy\conversion.py", line 337, in RasterToPolygon
raise e ExecuteError: ERROR 010151: No features found in c:\gis\output\control_test.shp.
Possible empty feature class. Failed to execute (RasterToPolygon).

Is it possible to fix this error or to use another function to get to my goal?

I have the original PNG as an SVG with vector drawn polygons if that helps.


texts: How can I make a heatmap of polygons using a graduated color scheme in QGIS?

texts:
I want to make a heat map of species biodiversity like the people did in the picture below. I can sort of get it to work by setting my layer with multiple polygon subclasses in the Layer styling panel to "Single" instead of "Graduated", and then modifying the opacity. However, this only allows me to use one color, and the areas with few species are really dim. I want to make a graduated color scheme of multiple colors. I have hundreds of species, so it isn't feasible to go and change the color of every single polygon layer. I also don't want to rasterize the polygons into a grid to make a series of grid squares. I want it to be in vectors

This is what my POLYGON layer looks like. It is not a series of points, it is a series of polygons. The example map I attached is made up of overlapping polygons too, not points. Though some spots may look like points in their map, they are actually circular polygons


texts: Speeding up searches on layer with many attributes?

texts:
I'm using QGIS 3 with my data in a PostGIS server on the same local network. I want to quickly search for a feature based on a certain attribute or attributes of my layer. This layer has ±2 million features and 74 columns. I don't want to search in most of the columns and would be happy specifying which columns to search in.
I've added indexes to the columns that I want to search on. Searching via a SQL query returns practically immediately (±100ms).
I've tried a few methods of searching in QGIS, but keep running into performance much slower than it could be:

Locator bar
I've tried the locator bar per this answer. I select the layer I want in the Layer panel, then type f search string into the locator bar. I see some results quickly when I start typing the search string, but when I type more it stops showing results, and I have to wait 15-30 seconds before it finally finds the match.
Logging SQL statements on my server, I can see why: it creates a query of the type
SELECT [each of the columns]
FROM thetable
WHERE ((((textcolumn1) ILIKE '%search string%')
OR ((textcolumn2) ILIKE '%search string%'))
OR ((textcolumn3) ILIKE '%search string%'))
LIMIT 30

I can speed this up on a table with a small number of attributes, by adding trigram matching indexes with the trgm module, but with a large number of attributes, this becomes a sequential scan. It seems Postgres doesn't realize that, by limiting to 30 results, it could search the indexed columns first, then only look at the others if it didn't get to 30. Though even if Postgres could be that clever, that wouldn't solve my problem if I enter a search term that matches only one result!
Much better would be if I could specify which attributes to search.

Search Layers plugin
When I use the Search Layers plugin, I get the match(es) practically immediately if I specify both the layer and the attribute. Great! But if it's set to "<All Layers>", or "<All Fields>", it takes a long time to return anything, and, if I close the plugin window and open it later, it defaults back to "<All Layers>". I don't want to have to change the settings every time I open the window to avoid this "gotcha".

MMQGIS Search
I tried the MMQGIS Plugin from this answer to the same question. It uses a panel that remembers the previous settings and asks for what attribute(s) to search. But instead of running a SQL query to search those columns, it runs a SQL query to load the entire table into local memory and then searches the local memory copy. This is unusable for the use case of large tables - takes much longer than Locator bar, and freezes QGIS while it's running.

Question: Is there a current performant solution to easily search a layer, with many features and many attributes, for certain attributes?


texts: How to color different regions that are on one single layer

texts:
I received the following files:

.shp .shx .dbf

I opened the shapefile using QGIS 2.18 and I see there only one layer. I did Add vector layer.

I'd like to be able to color different regions with different colors.
Is this possible, having one single layer? How? (Even if by using another software)

I'm new to this. I did Add vector layer for importing. Maybe I could do this by importing the file differently, and not have one single layer.

Or maybe I could modify the .dbf file and add a new column for colors. I don't know.

My .dbf file looks like this:


texts: How to use PHP MapScript to generate a MapServer Mapfile

texts:
I just read from this question that the OP was

  using Mapscript and C# to generate and display a Mapfile

I was wondering if it is really possible to use MapScript to generate a Mapfile? And if so, is there an example of doing so with PHP MapScript?

P.S. I am new to MapServer, and I feel that manually writing every mapfile's is a pain even if I learn all its grammar. It would be nice if a bunch of mapfiles can be generated from a program by varying certain parameters.

If MapScript isn't the solution, is there any other way to generate mapfiles programmatically?

P.S.2. This post asks a similar question about doing it in python, but there seems to be no solution in the answers.

Thanks in advance,

--- EDIT ---

Just to clarify, I am interested in writing/serialization of Mapscript objects into a mapfile (in PHP), not writing a mapfile verbatim via a programming languages such as PHP or C++. From the SO question quoted, this seems possible in C#. I was asking whether there are specific examples of actually doing that (preferably in PHP).


texts: Plotting multiple raster stacks with rasterVis::gplot and ggplot2::facet_wrap

texts:
I'm trying to ggplot2::facet_wrap a stack of rasters, and then call geom_raster to add more rasters on top of the previously mapped raster stack.
Here is a very simple example dataset to re-create what I've tried:
    #Load packages
    library(raster)
    library(ggplot2)

    # create a raster
    set.seed(11)
    r <- raster(nrows = 10, ncols = 10, res = 30, xmn = 267195, xmx = 267375, ymn = 4016985, ymx = 4017135)
    r <- setValues(r, runif(ncell(r), min = -10000, max = 10000))

    # create a raster stack
    rr <- lapply(1:4, function(i) setValues(r, runif(ncell(r), min = -10000, max = 10000)))
    r_stack <- stack(rr)
    crs(r_stack) <- "+proj=utm +zone=43 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"
    prj <- crs(r_stack)
    plot(r_stack)

    # use gplot w/ facet_wrap
    library(rasterVis)
    gplot(r_stack) +
      geom_raster(aes(fill = value)) +
      facet_wrap(~ variable) +
      scale_fill_gradientn(colours = rev(terrain.colors(225))) +
      coord_equal()

   # create another raster stack
    ss <- lapply(1:4, function(i) setValues(r, runif(ncell(r), min = -1, max = 1)))
    s_stack <- stack(ss)
    crs(s_stack) <- "+proj=utm +zone=43 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"
    prj <- crs(s_stack)

# plot r_stack and call geom_raster to plot s_stack on top
    gplot(r_stack) +
      geom_raster(aes(fill = value)) +
      facet_wrap(~ variable) +
      scale_fill_gradientn(colours = rev(terrain.colors(225))) +
      coord_equal() +
      geom_raster(data=s_stack, aes(x=x,y=y))

#Error: `data` must be a data frame, or other object coercible by `fortify()`, not an S4 object with class RasterStack
#Run `rlang::last_error()` to see where the error occurred.

I realize that the new raster would essential cover up the previous raster plot, but the raster data that I'm calling in the 2nd geom_raster call will have different data to display, most of which will be 'gaps' or NA.


texts: Loading shapefile into PostGIS: ${TRAINING_ROOT} Directory

texts:
I try to follow this tutorial (link below: Loading shp file into PostGIS), however I dont understand how do I get to the $TRAINING_ROOT directory? It always tells me: No such file or directory.
https://docs.geoserver.geo-solutions.it/edu/en/adding_data/shp_postgis.html


texts: Geocoding on OSM shapefiles via Elasticsearch

texts:
I want to get a geocoding system based on OSM data and Elasticsearch engine. I've tried some projects like komoot\photon, but, as I'm a really newbie in this, I stuck on preparing of the data for my area and many troubles with a Nominatim server.

I could be using the whole world dump of indices but it is too big data for me(about 56Gb space). However, I found out I can use just Elasticsearch from scratch and some tool such as ogr2ogr which provide a support of a files converting from shape to elasticsearch. So, I downloaded shapefiles from geofabrik and converted some of them to Elasticsearch.

What kind of shapefiles (points, roads, places...) should I import to Elasticsearch for organize simple geocoding?

I'm thinking of using it on my Leaflet HTML page.

And what kind of requests should be to work with Elasticsearch engine then?


texts: Problems with rasterization through QGIS Modeler

texts:
My goal is to rasterize a vector polygon-based shapefile.

Whenever I run in the modeler or even from the Processing Toolbox the Rasterize tool of GDAL, I receive a raster layer that's not good for me, as in "NaN" values.

To be more precises, what does my Model do? It first takes in a vector, then through Advanced Python Calculator I've placed a function so that it reclassifies the string value to integer, which is added to a new field (I've tested the output right after the Calculator step and it's fine, however it does have a lot NULL values because that data is not of interest but perhaps that could be the problem).

When it converts the vector into raster based on this new field, it basically gives me a raster with "NaN" values, as in pure black. Also, I plan on obtaining raster oa cellsize of 10.

Any ideas on how to fix this?


texts:
I'm newbie to programming in pyqgis.

I need to perform a simple "GeoTiff * A" type of calculation(where A is a constant), but I am struggling to do so.

I tried processing.runalg('saga:rastercalculator') but failed 'coz my constant is NOT a raster.

I suppose I should be using 'qgismodelertools:calculator' or possibly something like to what they use here http://geoexamples.blogspot.it/2012/12/raster-calculations-with-gdal-and-numpy.html

Any clues?

Gene's answer iss succinct and logical but rasterio does not come out of the QGIS 2.8 box, which I'm using from OSGEO4W.

In the meantime, I used EikeMike's useful code by turning it into a working, callable function:

def raster_calc(RGBI_DS,outFile,constant):

  import processing, os, time, subprocess, array
  from osgeo import gdal
  from osgeo.gdalnumeric import *
  from osgeo.gdalconst import *

  # Read Raster dataset
  GDALRgbiDS = gdal.Open(RGBI_DS, GA_ReadOnly )
  band1rot = GDALRgbiDS.GetRasterBand(1)

  # Convert band to Array
  Rdata = BandReadAsArray(band1rot) * 1.0

  # calculate raster * constant
  dataOut = (Rdata * constant)

  # Export result as Geotiff
  driver = gdal.GetDriverByName("GTiff")
  dsOut = driver.Create(outFile, GDALRgbiDS.RasterXSize, GDALRgbiDS.RasterYSize, 1,
   gdal.GDT_Float32)
  CopyDatasetInfo(GDALRgbiDS,dsOut)
  bandOut=dsOut.GetRasterBand(1)
  BandWriteArray(bandOut, dataOut)

  #Deletion of variables, this actually triggers the writing of the output file
  dsOut =None
  bandOut = None
  driver =None

Later, I found myself a much easier and shorter solution.

Here it is:

Stick to 'saga:rastercalculator' (or similar algorithm). Simply build the mathematical formula therein as a new variable where your constant is itself a variable. In other words:

 import processing

 formula = ('a * %s' % (str(constant))) # where 'a' is going to point to yourraster   and 'constant' is the numeral you want to multiply your raster by
 processing.runalg('saga:rastercalculator',yourraster,None,formula,outFile)


texts:
I am trying to split a street polygon based on the segment ID of the street, giving them the same segment_ID in the attribute table of the street lines. I did not use the line feature class to derive the polygon but instead used a polygon feature class that included the width of the roads (but was not divided into segments).

I tried using the split tool, but it doesn't accept the street line segment as a feature.


texts: I'm new to QGIS, and I'm wanting to export a world map that I can then project onto a sphere in After Effects. I'm just wondering what kind of CRS projection I would need to make the 2D map in, so that it all lines up properly when projected onto a sphere.


texts: ogr2ogr from personal geodatabase to spatialite produces "non-spatial" tables

texts:
I have around 4 GB of data in 10 separate Personal Geodatabases, each of which contains multiple feature classes. I'm trying to convert the smallest one to SpatiaLite via ogr2ogr, in preparation for converting them all. The geodatabase contains two feature classes, both of which are linestrings. Here is my command line:

ogr2ogr -f SQLite -dsco SPATIALITE=YES Utility.sqlite Utility.mdb

But I get inconsistent results when I attempt to display the imported data in QGIS.

If I create a connection to the SpatiaLite database and then open a layer via the SpatiaLite node in the browser, I get non-spatial tables. Specifically, the DB manager tells me "This is not a spatial table".

But if I load the same layer directly from the database's location on the file system (within the browser) it displays perfectly.

I've been experimenting with different ogr2ogr arguments and I've also tried running RecoverGeometryColumn() on the imported data, but I can't make the thing display properly. (RecoverGeometry returned an error about "no such column", but it's actually referring to the table name. I'm slightly confused.)

I have created several other SpatiaLite databases within QGIS or Spatialite_GUI, which do work properly. But I'm not familiar enough with the database schema to identify any obvious differences between the ones that work and the ones that don't.

I'm migrating from ArcGIS (having not used it for a few years) and I'm very new to the rawness of relying on SQL commands (and command line statements in general) to manage data, so I'm happy to accept that I've missed something obvious. I'd just love to know what it is.


texts: Using ways table features of OpenStreetMap data

texts:
I have made this web-app http://160.75.133.173/openWebApp/routing.jsp
Users just have to click two points to figure out the shortest possible route available connecting them. I am using OSM data and pgrouting extension in postgresql.
But the OSM ways table is missing any feature which could be used to rank the roads (like Highways, one-lane, two-lane etc, or any sort of ranking depending upon road's width or any other feature).
How do I deal with questions like "I don't want the shortest route but the route using Highways only or as much Highways as possible", using OSM data?


texts: How to edit (particularly flip direction of) polyline obtained through UpdateCursor?

texts:
I want to flip the direction of a polyline feature obtained as a result of an Update query; e.g.

lines= arcpy.UpdateCursor(myFeatureClass,myQuery)
dsc= arcpy.Describe(myFeatureClass)
try:
    for ln in lines:
        # ? Check number of parts
        # ? Flip direction of each, or some of the parts
        lines.updateRow(ln)
except:
    print "Error:", sys.exc_info()[0]
finally:
    # releasing locks
    if lines: del lines
    if ln: del ln

BTW, any better way to know how many rows are returned than counting them?


texts: gdal_merge: ERROR 1: Attempt to create 590x-369 dataset is illegal,sizes must be larger than zero

texts:
I am attempting to merge 600+ .tif files using gdal_merge.py (GDAL 2.2.2, released 2017/09/15) but I'm getting this error:

  ERROR 1: Attempt to create 590x-369 dataset is illegal,sizes must be
  larger than zero.

I have created all the files with Python and individually they display correctly in QGIS. If I gdalinfo one of the .tif files I get this:

Driver: GTiff/GeoTIFF
Files: 561700_186000.tif
Size is 10, 10
Coordinate System is:
PROJCS["OSGB 1936 / British National Grid",
    GEOGCS["OSGB 1936",
        DATUM["OSGB_1936",
            SPHEROID["Airy 1830",6377563.396,299.3249646,
                AUTHORITY["EPSG","7001"]],
            TOWGS84[446.448,-125.157,542.06,0.15,0.247,0.842,-20.489],
            AUTHORITY["EPSG","6277"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4277"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",49],
    PARAMETER["central_meridian",-2],
    PARAMETER["scale_factor",0.9996012717],
    PARAMETER["false_easting",400000],
    PARAMETER["false_northing",-100000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH],
    AUTHORITY["EPSG","27700"]]
Origin = (561700.000000000000000,186900.000000000000000)
Pixel Size = (100.000000000000000,100.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  561700.000,  186900.000) (  0d19'58.22"E, 51d33'26.27"N)
Lower Left  (  561700.000,  187900.000) (  0d19'59.88"E, 51d33'58.61"N)
Upper Right (  562700.000,  186900.000) (  0d20'50.11"E, 51d33'25.23"N)
Lower Right (  562700.000,  187900.000) (  0d20'51.78"E, 51d33'57.58"N)
Center      (  562200.000,  187400.000) (  0d20'25.00"E, 51d33'41.92"N)
Band 1 Block=10x10 Type=Float32, ColorInterp=Gray

Which seems OK. Any ideas why I can't merge these files?


texts: Qgis: work in pyqgis from an external idle

texts:
I'm writing a python script but, at the moment, I do everything from the Qgis python console.
Is there a way to write and run a pyqgis script externally from Qgis?

I'm using Qgis 2.18 and python 2.7.12 on Ubuntu 16.04 OS.


texts: Adjusting automatically angle before printing Strip Map from Data Driven Pages?

texts: Is it possible using Data Driven Pages automatically adjusting angle to the print frame area before exporting to pdf? I have blocks (red marked, a3 format, 1:1000 scale) where each is in different angle; please see below:


texts: Part of layer disappears when printing ArcGIS 10.0 as pdf

texts:
There are a couple of maps i made that I want to save as pdf in ArcGIS 10.0. For three of the maps it worked fine, but there are problems with the other two.

The first one has a layer that consists of four topographic maps. Two of the maps are not printed completely. The second map also has a layer of four topographic maps, and a horizontal line is drawn where the upper two overlap the lower two.

I have no clue how to prevent this. I am printing with the 'windows printer' printer engine. I've tried both Adobe PDF and PDF creator.

An update, what I've tried now:

Someone told me to try Repair Geometry, but I think that's
impossible because the topographic maps are .tif files
MXDDOCTOR did not find any corrupt files, I am trying to make a new mxd in the
'treatment' part to see if it worked anyway, but it is taking ages
Making a new .mxd file, not entirely new but based on a different mxd does not work
Making the map on a smaller scale (zoomed out) does not work
Exporting or printing to pdf at a lower resolution does not work
I tried doing the same while turning of the working tifs to test memory, still the same problem.


texts: Projected MODIS data not aligning with ESRI basemap in ArcGIS for Desktop?

texts:
I can not get 1km MODIS data aligned with ESRI basemaps in ArcGIS 10.2.2 for Desktop.  See the photo below.

I reprojected MODIS17A3 data using the MRT (MODIS reprojection tool). I selected WGS84 under Projection Parameters. I tiled two files (shouldn't matter). My shapefile is in GCS_WGS_84 and it aligns perfectly over the islands.  It seems that the MODIS data should move south and east. Generally, it is aligned, but there are some grid cells completely in the water.

I have several study areas in the Great Lakes region. Another is pictured at the bottom.  Any ideas?  These are quite north of the equator.


texts:
I have a large database in Filemaker. To get this joined with a feature class in ArcGIS I have to make an ODBC connection. I do this in the 32-bit ODBC in windows because ArcGIS is 32-bit. In Arccatalog I then have to make an Oledb connection to the ODBC to get the table translated in ArcGIS.

This seems to work up until the point I want to add symbology based on the joined table. The text fields can't be used to add symbols. The number fields do work and the table does show the values correctly.

some thoughts:

the text fields sometimes contain 1000000 characters
could there be an office involvement?

What could be going wrong?

ArcGIS 10.3.1
Filemaker 14
Windows 7 (64-bit)


texts: How to replace the Nan or nodata values into 0?

texts:
I have a raster in .img format and i would trasform some values (Nan) into 0 with QGIS 2.14.3 (Windows).
I've tried different ways but I didn't have results:
1- with raster --> conversion --> conversion format
2- operations with raster calculator (to replace non-zero values).
With this operations, Nan values have been replaced by Nodata. After I have repeated the same operations but no data values haven't been replaced.


texts:
I write because I can’t use the plugin “Geotag and Import photos”. My version is 2.18.13 Las Palmas and it runs on WIN 7 Home.
First of all, when I click on “Import photos” and specify the path to photos, I obtain the following error message:

The details are:

File "C:/Users/Emazza/.qgis2/python/plugins\geotagphotos\importphotosdialog.py", line 88, in selectPhotosDir
            self.loadTags(dirName)
          File "C:/Users/Emazza/.qgis2/python/plugins\geotagphotos\importphotosdialog.py", line 125, in loadTags
            with et:
          File "C:/Users/Emazza/.qgis2/python/plugins\geotagphotos\exiftool.py", line 179, in __enter__
            self.start()
          File "C:/Users/Emazza/.qgis2/python/plugins\geotagphotos\exiftool.py", line 162, in start
            stderr=devnull)
          File "C:\OSGeo4W\apps\Python27\lib\subprocess.py", line 711, in __init__
            errread, errwrite)
          File "C:\OSGeo4W\apps\Python27\lib\subprocess.py", line 948, in _execute_child
            startupinfo)
        WindowsError: [Error 2] Impossibile trovare il file specificato

The pointer begins to rotate as if it is working; nevertheless, I can still use the software.
Even if I succeed to insert the output shapefile, I can’t see any EXIF tag.

I already installed Exiftool and I can’t understand which is the problem.
Any suggestions for these two problems?


texts: I need to change the name of a field in an attribute table. I just installed QGIS 2.18.3. I understand that I should now use "Refactor Fields" instead of Table Manager. I don't see it anywhere, including in my Toolbox under Processing. I understand that I may need to be in the Advanced Interface mode. However, I don't see how to get into it. It is not at the bottom of my Toolbox, as I have seen mentioned.


texts:
I'm serving a NDVI coverage via Mapserver. So, one band, float32 values, [-1.0, 1.0]. I want to display the NDVI in greyscale, and generate a legend. I've tried 2 things:

Create two dummy classes, just for testing:

So, my legend looks like this

And the image, obviously, is red and blue: http://dl.dropbox.com/u/6599273/capturas/ndvi_en_gvsig.png

But I'd like to see the image in greyscale... The I've tried the next thing

Pre-classify values with PROCESSING "SCALE=-1.0, 1.0". I don't know what Mapserver does internally with this processing order, but I can see the image in greyscale: http://dl.dropbox.com/u/6599273/capturas/ndvi_gris_en_gvsig.png

The problem is now I can't see any legend, because I've not created any class. I'd like to see all the colors of the grayscale image in my legend.

How could I do it?

Many thanks in advance


texts: Extruding buildings without them floating in qgis2threejs

texts:
I am trying to extrude building footprints relative to the DEM. But buildings seem to be floating in the air.

height_2 represents the heights of the buildings.


texts: What's the easiest way (style wise) to publish my MapInfo map as WMS?

texts:
I've downloaded data for my country from OSM and let a GIS expert to style it for me in MapInfo pro. (I have TAB files and .wor - workspace).

When I tried converting them to SHP (for QGIS) I've noticed that the styles don't transfer with the conversion (they are saved in the workspace in MapInfo Pro). So I manually re edited the SHP file to match the style.

Next , When I imported the SHP files to GeoServer , Again the styles from the SHP files didn't transfer - and here I'm stuck. I'm using windows based machine with QGIS and MapInfo. I have GeoServer installed (and can also install other WMS like MapServer if needed).

I'm looking for the easiest way to set styled layers in any MapServer so I could install "offline" map based web app on client site (using open source WMS)

Thx.


texts:
I'm working on QGIS with a database created in PostGIS and a shapefile. Inside the database, I have geometry type data about location of some points, instead on the shapefile there some selected areas.

With these 2 kind of data, I would like to know if there are some points inside or close to the shapefile area. For example: is there a point inside a selected area?

Can I select a part of the total area to see points also near to it?

Is there a query tool to use?

How can I connect these data created in different layers?


texts: Disk space and numbers of tiles for world layer

texts:
I'm thinking about testing tiles generation of a world layer for webmapping.

In this context, I want to know how I can calculate the number of tiles (256x256 for tilecache if I remember) I will need for n levels in WGS84 coordinates.

Maybe I will use a composite solution with pregenerated tiles and on the fly generation but I have to know at the end the available disk space if every tiles are generated and the number of files it will represent.

There are two goals, time it will take and disk space needed.

Any informations are welcome

Edit:

I've found this script to calculate the number of tiles to generate depending on zoom levels and your data extent. See this gist https://gist.github.com/1675606


texts: Error in Select Numeric type column in pyodbc with ArcCatalog?

texts:
I created a python script that uses "pyodbc" plugin to connect a SQL Server database and write out data. It was imported in an ArcGIS 10.1 tool and run from there.
try:

arcpy.AddMessage('FLAG 0')
cnxn = pyodbc.connect('DSN=' + DSN + ';UID=' + UID + ';PWD=' + PASS)
cursor = cnxn.cursor()

arcpy.AddMessage('FLAG 1')
# Retrieves all the records from the All Obs for KML view
if cnxn is None:
    arcpy.AddMessage('cnxn is None')
else:
    arcpy.AddMessage('cnxn is not None')

if cursor is None:
    arcpy.AddMessage('cursor is None')
else:
    arcpy.AddMessage('cursor is not None')

cursor.execute('SELECT LocationName, Lat from VW_Obs_for_KML_Exp')
for row in cursor:
     arcpy.AddMessage(row.LocationName + ', ' + str(row.Lat))

arcpy.AddMessage('FLAG 2')

except (pyodbc.IntegrityError, pyodbc.ProgrammingError, pyodbc.Error):
    arcpy.AddMessage("PYODBC Error. Details: " + str(sys.exc_info()[:2]))
except TypeError, e:
    arcpy.AddMessage("TypeError. Details: " + str(sys.exc_info()[:2]))
else:
    arcpy.AddMessage("No error occurred.")

The script accesses a Table and outputs 'Lat' as Numeric(38, 8) data type. It worked perfect the first time from ArcCatalog. However, it threw out an error the second time and onwards. Detailed error stack trace is as follow:

FLAG 0
FLAG 1
cnxn is not None
cursor is not None
Runtime error
Traceback (most recent call last):
  File "<string>", line 1, in <module>
  File "C:\Python27\ArcGIS10.1\Lib\decimal.py", line 543, in __new__
    m = _parser(value.strip())
TypeError: 'NoneType' object is not callable
TypeError. Details: (<type 'exceptions.TypeError'>, TypeError("'NoneType' object is not callable",))
Process is complete.

So I had to restart ArcCatalog to successfully run it again. By the way, it also worked well when run from Python in the command line prompt. So this error only occurred when run from ArcCatalog.

Is there anyone who has worked with "pyodbc" plugin and experienced similar problem before?


texts: LoCoH or TLocoH home range estimators in ARC?

texts:
Is there a way to create a Getz and Willmers Local Convex Hull in ArcMap? The "Minimum Bounding Geometry" tool is insufficient for my purposes, as it greatly overestimates the sizes of home ranges.

From my research, I've found that there used to be an extension for Arc 9, but I am using ArcMap 10.1.


texts: How do I convert a DEM GeoTIFF to Contour Shapefile

texts:
I have a digital elevation model as a GeoTIFF and I'm looking to generate a Shapefile of contours with 5m bands on the command line. Is this possible?

Ideally using GDAL and Python.


texts: GEE error when extracting image timestamp

texts:
I have filtered Sentinel 2 images to get suitable cloud free images from the time period I require, however I cannot figure out how to extract the date from the images, but with my current code I get the error:

  Date (Error) Date: Parameter 'value' is required.

//  Define the image collection
var image = ee.ImageCollection(sent2
// filter to get only images in the date range we are interested in
    .filterDate("2018-05-01", "2018-07-30")
// filter by polygon
    .filterBounds(tiwi))
// sort the collection by cloud cover
    .sort("CLOUD_COVERAGE_ASSESSMENT")
// select the first image i.e. the most cloud free image in the date range
    .min();

// print the image to the console.
    print("A Sentinel-2 scene:", image);

  // Define visualization parameters in a JavaScript dictionary for true colour rendering. Bands 4,3 and 2 needed for RGB.
    var trueColour = {
        bands: ["B4", "B3", "B2"],
        min: 0,
        max: 3000
        };

  // Add the image to the map, using the visualization parameters.
  Map.addLayer(image, trueColour, "true-colour image");

// Get the timestamp and convert it to a date.
var date = ee.Date(image.get('system:time_start'));
print('Timestamp: ', date);


texts: Using multiple map layers with different extents. Any suggestions.

texts:
I have three layers I need to try match up at the borders. This is not a projection problem, all layers are projected to Albers NAD83. The base layer (Green) is 30m DEM, with a US boundaries layer on top of the DEM and finally a polygon "layer of interest" on top (Red). As you can see in the figure below, none of these layers match up at the borders.

What is the best way to get these layers to match up in ArcMap 10.1?
This is mainly for aesthetics, as I have it the red layer is transparent and it just looks bad.

EDIT:

So figure 2 (bottom) is what I have done thus far. (Note: I have changed the colors around in figure 1 just to make the different layers easier to distinguish for you all.)

As you can see, the problem with the layers lining up isn't all that bad at this scale.. but I'm picky like that. Figure 1 is zoomed in to an area on the coast of South Carolina (apologies for no scale bar).

@Kevin, the state boundaries were requested because this project spans many states and the hillshade (while not required) was also requested and adds some pizaz.

Thanks everyone for comments so far.


texts: How to define a css class for a delete button in Openlayers?

texts:
Define a css class for delete button in openlayers?

        .olControlDelete {
     background-image: url("help.png");
     background-repeat:no-repeat;
     background-attachment:fixed;
     background-position:center;
     z-index:10;
      }

Or

  div.olControlDelete {
  top: 100px;
  left: 50px;
  position: absolute;
       }

  .olControlDelete div {
  background-image: url("edit_sprite.png");
  background-repeat: no-repeat;
  width:  22px;
  height: 22px;
  border: thin solid black;
  margin-top: 10px;
}

Don,t seem to work, Please help. Last question for openlayers. please

Below Code from Vadin:

                 deleteShape = new OpenLayers.Control.SelectFeature(
polygonLayer,
{
    displayClass: "olControlDelete",
    title: deleteFeature,
    eventListeners: {
        featurehighlighted: function overlay_delete(event) {
            var feature = event.feature;
            if( confirm(strConfirmMessage) ) { // prompt
                polygonLayer.removeFeatures( [ feature ] );
                // TODO other clean up - remove popups or send some ajax to server.  Up to you....

            }
        }
    }
}
       );

            map.addControl(deleteShape);


texts: GDAL error in R path.expand(path): invalid 'path' argument

texts:
Using windows 10.

I would like to stick to using R since that's where a lot of my pre-processing and analysis will take place for a bigger project.

I have downloaded a number of MODIS files for the MOD13A3 product, they are in hdf format. The end goal is to create a loop to convert the hdf files into a tiff, scale, clip and project the NDVI and EVI layers.

I am getting stuck before I even attempt those steps with this strange error when I try to use some of the gdalUtils functions:

Error in path.expand(path) : invalid 'path' argument

I get this when trying:

gdalinfo("C:\\Users\\mydirectory\\R\\MOD13A3.A2015335.h13v03.005.2016007192527.hdf")

I have also tried using different types of paths like:

gdalinfo("MOD13A3.A2015335.h13v03.005.2016007192527.hdf")

I even get this error when trying: gdal_setInstallation() and get_subdatasets

Additionally, gdal-config produces an error. I suspect this is due to a weird installation issue and environment path if someone could provide a bit of guidance!

I have my working directory set to the same folder where these files are found.

I haven't been able to find much on this error so I'm hoping somebody here can help me.


texts: Why are maps made on Google Fusion Tables, using KML file format, not rendering properly?

texts:
On this particular map of mine, the feature in Fusion Tables seems to only approximate the geo-locators that I have specifically inputted.

In fact, it takes every second point and plots it out, making some very weird shapes.

At first I thought it was an issue with the coding. But after repeatedly checking back, it turned out the points had no problem, and they aligned perfectly.

The display was messing up, especially the "preview KML" option, although things usually get blurred in there and I don't use it that often, so this should be fine.

It doesn't even go over Google's own maximum "limit" of a million characters per element's object.

To clarify, however, this happens even if/when I refresh, go full-screen, publish, or share the document-type.

I've been using Fusion Tables for a long time, yet it's never given me this problem before. What's going on, how can its "problem" be fixed, and what should be done, in order to return this issue to its original state?


texts:
Since QGIS 3, there are two vertex tools: current layer and new: all layers. Normally I was able to access all actions with the code below. But now I can't find and activate the vertex tool (current layer). Even not by searching its caption! I have no idea how an action could be added which "does not exist". Any ideas?

def findAction(objName):
    try:
        toolbars = iface.mainWindow().findChildren(QToolBar)
        for toolbar in toolbars:
            for action in toolbar.actions():
                if isinstance(action, QWidgetAction):
                    if action.defaultWidget() and action.defaultWidget().actions():
                        dwa = action.defaultWidget().actions()
                        for a in dwa:
                            if IsSearchedAction(a,objName):
                                print ("a",a.objectName())
                                return a
                else:
                    if IsSearchedAction(action, objName):
                        print ("action",action.objectName())
                        return action
        menubar = iface.mainWindow().menuBar()
        for action in menubar.actions():
            if action.menu():
                for action in action.menu().actions():
                    print("menu")
                    if action.menu():
                        for a in action.menu().actions():
                            if IsSearchedAction(a, objName): return a
                    else:
                        if IsSearchedAction(action, objName): return action
            else:
                if IsSearchedAction(action, objName): return action
        return None
    except Exception as e:
        print (e)

def IsSearchedAction(action, objname):
    objname = objname.lower()
    caption = action.text()
    caption = caption.replace("&", "")
    caption = caption.replace(".", "")
    caption = caption.strip()
    caption= caption.lower()
    return action.objectName().lower() == objname or caption == objname

print (iface.mapCanvas().mapTool().action().objectName())
action ="mActionVertexToolActiveLayer"
action ="mActionMeasure"
action ="mActionSelectFeatures"
action ="mActionVertexToolActiveLayer"#not found
print(findAction(action))


texts: Calculating Distance from Outlet to Dam(s) using ArcGIS Desktop?

texts:
I am trying to carry out what should be a relatively simple analysis of river network fragmentation by dams. Similar questions have been asked as:

Obtaining length of river upstream of points up to next point using ArcGIS Desktop? and
Finding upstream locations (polygons) using ArcGIS Geometric Network?

What I would like to achieve is:

Calculate the distance from river outlet upstream to first dam on each river/tributary
Calculate the cumulative number of dams from the outlet to the uppermost point

The data I am working with is:

hydrosheds 15 arc second river network. This data set includes IDs ("arcid") for numerous small segments of the polyline river network.
A point dam data set

Many answers, like a comment on Finding distance from river outlet to points, following river, using ArcGIS Desktop?, refer to a computer program, RivEX, and end there. I acknowledge the potential utility of this program, but I would like to work through the steps and run these calculations in ArcGIS Desktop.

Is creating a geometric network the right approach in this case or is there another tool that could use the built-in hydrosheds line ID and calculate distance from outlet to point (dam)?

In essence, I am unsure of how to work with the arcid field (specifying the relative position in the watershed) and represent the dams as nodes.


texts: I have a single shapefile full of overlapping polygons. The ArcGIS Union tool will generate individual, non-overlapping polygons with that geometry. Is there a way to do this with GDAL/OGR or another Python library that is not made by ESRI?  The ogr.union function in OGR seems to be more like a 'dissolve' in ArcGIS.


texts: what is the suitable projection coordinate system for Iran

texts:
I am doing analysis for province wise Iran. I am using WGS1984 as my geographical coordinate system, and I do not know which is suitable Projection coordinate system for entire Iran.

I feel the UTM will not be suitable because the country shares several zones in UTM system.


texts: How to calculate the angle at which two lines intersect in PostGIS?

texts:
I want to calculate the angle between two lines where they intersect in PostGIS.

The starting point for angle calculations in PostGIS seems to be ST_Azimuth - but that takes points as input. My first thought was to take the endpoints of the intersecting lines and performing an Azimuth calculation on those. That is not good enough, because most of the line features are not straight, and I am  interested in the angle at intersection. So what I came up with is a nested operation that goes through the following steps:

Identify all the intersections between the two line feature tables.
Create a very small buffer around the intersection point
Identify the points where the line features intersect the buffer exterior (taking the first point if there are more than one - I'm really only interested in whether the angle is close to 0, 90 or 180 degrees)
Calculate ST_Azimuth for those two points.

The full SQL is kind of long to post here, but I gisted it here if you're interested. (By the way, is there a better way than to carry over all the fields going down the WITH statements?)

The results don't look right, so I'm clearly doing something wrong:

EDIT I redid the calculations in EPSG:3785 and the results are a little different but still not right:

My question is where the flaws are in this process. Am I misunderstanding what ST_Azimuth does? Is there a CRS issue? Something else altogether? Or maybe there's a much, much simpler way to do this?


texts: Identifying number of overlapping polygons in QGIS

texts: I'm looking for a possibility to identify the intersecting parts in one single polygon layer with ten features, where at least four polygons overlap. How can I do it in QGIS?


texts: Esri Notebook - Select and append records to table

texts:
Using Notebook editor in ArcGIS Online i want to extract from a FC table, those records with a specific data and then insert them into a separate table. I am trying the code below, but getting an "TypeError: 'FeatureSet' object is not subscriptable" error.
from arcgis import features
import pandas as pd
from copy import deepcopy
FC_Item = gis.content.get("a1e115946edf47cb8128b8a6cae4b4ec")
ROAD_Records = FC_Item.layers[0]
Table_Records = FC_Item.tables[0]
Current_Sweep_Records = ROAD_Records.query(where = 'Last_Swept_Date > DATE \'2021-04-19\'')
#print(Current_Sweep_Records)

features_to_be_added = []

# get a template feature object
template_feature = deepcopy(Current_Sweep_Records[0])

# loop through each row and add to the list of features to be added
for row in new_rows.iterrows():
    new_feature = deepcopy(template_feature)

    # assign the updated values

    new_feature.attributes['Road_Name'] = row[1]['Road_Name']

    #add this to the list of features to be updated
    features_to_be_added.append(new_feature)

Table_Records.edit_features(adds = features_to_be_added)