texts: PostGIS-syntax: export query-result to QGIS / Shapefile

texts:
Using PostGIS: syntax-basics for function, I was able to solve a long-standing problem by creating a geometry, using a SQL-query with PostGIS. However, I want to use the created geometry in QGIS and have problems exporting it.

QGIS 3.8.2. Zanzibar, Windows 10, PostgreSQL 11.5.1., PostGIS 2.5.2

The query I used was as follows - it creates a Medial Axis of a Polygon (Shapefile bufferrail imported to PostGIS with shp2pgsql-gui.exe):

Select ST_ApproximateMedialAxis(geom) from public.bufferrail;

The result is shown in the Geometry-viewer-tab in pgAdmin4 and it looks like what I wanted to accomplish. However, I failed to export this result into QGIS. What I already tried:

1.) Using database-manager in QGIS, when I perform the query there and try to load the result as a new layer (load-Button), I get an error-message - see the protocol-entry (in german, unfortunately, it basically says Invalid PostgreSQL-Layer and WARNING: incorrect query):

Ungültiger PostgreSQL-Layer
WARNING    Fehlerhafte Abfrage: SELECT "" FROM (SELECT row_number() over ()
AS _uid_,* FROM (SELECT ST_ApproximateMedialAxis(geom) from public.bufferrail
             ) AS _subq_1_
             ) AS "subQuery_0" LIMIT 0 ergab 7 [FEHLER: Bezeichner in Anführungszeichen hat Länge
null bei »""«
             LINE 1: SELECT "" FROM (SELECT row_number() over () AS _uid_,* FROM ...
              ^
             ]
2019-08-25T23:47:57     WARNING    Ungültiger PostgreSQL-Layer

2.) Again in QGIS database-manager, I created a view from the query result. When I try to add this view to map (always from QGIS DB manager - right-cick on the view in pubic), QGIS always crashes.

3.) Still in QGIS database-manager, I try to export to file (Geopackage) - however, when I add the saved file as a layer to QGIS, an empty layer with no geometries is added.

4.) In pgAdmin 4: I run the SQL-query mentioned above and than click Download as CSV. I get a text-file with cryptic content (when opened in an editor). So no idea how to use it.

5.) I tried to google export shp from PostGIS and tried to use pgsql2shp - however, all of the many expressions to include pgsql2shp in my query failed - again, no idea how correct syntax works here.

6.) By the way: is there no way to use a simple GUI-method for export - similar to import shp to PostGIS with shp2pgsql-gui.exe? I tried the Export-tab in the GUI there, but export works only with tables. So probably a solution would be to convert the result of my query to a table - but unfortunately, again all my efforts failed - wrong syntax, again!

To come to my question: I don't care much how I can export the result form a relativele simple query in PostGIS to opening it in QGIS. Probably solution 5 is the most promising: so if someone could help with correct syntax, connection my original query and the export-function with pgsql2shp. But any other solution is good, too.


texts: Are Dimensions/Attributes Other Than Time and Elevation Possible

texts:
I have Time Dimension Working with an Image Mosaic layer. Is it possible to add dimensions other than time/elevation? For example the files in my mosaic are named as follows where rtma or pris represents a "climate" attribute:

pristmax_20150115.tif
pristmax_20150115.tif
rtmatmax_20150116.tif
rtmatmax_20150116.tif

In my mosaic directory, I now have two regex property files (one for time dimension and one for climatesource dimension):

timeregex.properties containing regex=[0-9]{8}
climateregex.properties containing regex=^.{4}

I modified my indexer.properties property file look like so:

TimeAttribute=ingestion
ClimateAttribute=climatesource
Schema=*the_geom:Polygon,location:String,ingestion:java.util.Date,climatesource:String
PropertyCollectors=TimestampFileNameExtractorSPI[timeregex](ingestion),StringFileNameExtractorSPI[climateregex](climatesource)

I also tried modifying postgis_rasters.properties to have a ClimateAttribute:

#-Automagically created from GeoTools-
#Wed Nov 25 13:16:00 MST 2015
Levels=0.026578191679850746,0.026578191679850746
Heterogeneous=false
TimeAttribute=ingestion
ClimateAttribute=climatesource
AbsolutePath=false
Name=postgis_rasters
TypeName=postgis_rasters
Caching=false
ExpandToRGB=false
LocationAttribute=location
SuggestedSPI=it.geosolutions.imageioimpl.plugins.tiff.TIFFImageReaderSpi
CheckAuxiliaryMetadata=false
LevelsNum=1

Geoserver created a table in Postgis as expected:

1 | 0103... |   rtmatmax_20150115.tif | 2015-01-14 17:00:00.000000 |    rtma

2 | 0103... |   pristmax_20150115.tif | 2015-01-14 17:00:00.000000 |    pris

But using &climate=rtma or &climate=pris in my url doesn't effect which raster is retrieved, whereas &time=2015-01-14 or &time=2015-01-15 selects the raster from that day. How can I get the climate attribute to function?

--- UPDATE ---

After following @rovo's answer, the new dimension appears in geoserver. Here is a screen cap.

And here is an example the postres side (note my regex was wrong so I currently don't have what I meant to have under climatesource but this is just a test):

The part I'm not yet able to get working is geoserver seems to be ignoring the CLIMATESOURCE param in my url.

In the getcapabilities I see

<Dimension name="time" units="ISO8601"/>
<Dimension name="CLIMATESOURCE" units=""/>
<Extent name="time" default="2017-07-06T00:00:00Z">2016-01-01T00:00:00.000Z,2016-01-02T00:00:00.000Z...
<Extent name="CLIMATESOURCE" default="agdd_20170706_base_thirtytwo">agdd_20160101_base_fifty,agdd_20160101_base_thirtytwo.....

When using &time=2017-07-06&CLIMATESOURCE=agdd_20170206_base_thirtytwo I get a valid map back because 2017-07-07 happens to match the default for CLIMATESOURCE. When I change time to be anything else it doesn't return a map because it's still using the default "agdd_20170706_base_thirtytwo" even if I try overriding CLIMATESOURCE with a matching date.

--- UPDATE2 SOLVED ---

RoVo's answer now clarifies to use DIM_CLIMATESOURCE in the WMS request.


texts: Error using arcpy.da.SearchCursor and MakeFeatureLayer

texts:
im trying to write a script that allows me to iterate through my map layers, and for each layer, iterate over each feature and add a field that gives the location of the feature using my index layer as reference. In other words, I want to know where every single feature is in relation to the index layer. The index layer is comprised of many polygons called "zones". So every feature should have a field that tells me which zone it is in. This is what I have so far :

import arcpy
arcpy.env.overwriteOutput=True
mxd = arcpy.mapping.MapDocument(r"C:\Users\Daimon Nurse\Desktop\Zonetrial.mxd")
df = arcpy.mapping.ListDataFrames(mxd,"Main Map")[0]
lyrList = arcpy.mapping.ListLayers(mxd,"",df)
fieldname = "Zone"
type = "TEXT"
inlayer = "layeritem2"
inlayer1 = "Zoneitem"

#for layeritem in lyrList:
    #if layeritem.name == "Zone of Interest":
        #Zones = set([r[0] for r in arcpy.da.SearchCursor(layeritem, "ZoneName")])
       #print Zones
for layeritem in lyrList:
    if layeritem.name == "Zone of Interest":
        with arcpy.da.SearchCursor(layeritem, "ZoneName") as cursor:
            for Zone in cursor:
                print Zone
                arcpy.MakeFeatureLayer_management(Zone,inlayer1)
                print Zone
                for layeritem in lyrList:
                    if layeritem.name == "AllZones_Topology":
                        print "Nothing will happen to " + layeritem.name
                    elif layeritem.name == "UWI Campus Boundary":
                        print "Nothing will happen to " + layeritem.name
                    elif layeritem.name == "Zone of Interest":
                       print "Nothing will happen to " + layeritem.name
                    elif layeritem.name == "polytrial selection 5":
                        print "Nothing will happen to " + layeritem.name
                    else:
                        print layeritem
                        print Zone
                        arcpy.AddField_management (layeritem,fieldname,type)
                        arcpy.MakeFeatureLayer_management(layeritem,inlayer)
                        #arcpy.MakeFeatureLayer_management(Zone,inlayer1)
                        arcpy.SelectLayerByLocation_management(inlayer,"intersect",inlayer1)
                        arcpy.CalculateField_management(inlayer,fieldname,Zone)

print "Script Finished!"

However, I am getting this error:

(u'1',)
Traceback (most recent call last):
  File "C:\Python27\ArcGIS10.1\Lib\site-packages\pythonwin\pywin\framework\scriptutils.py", line 326, in RunScript
    exec codeObject in __main__.__dict__
  File "C:\Users\Daimon Nurse\Desktop\DFMPROJECT\Scripts\Zonesort10.py", line 20, in <module>
    arcpy.MakeFeatureLayer_management(Zone,inlayer1)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.1\arcpy\arcpy\management.py", line 5774, in MakeFeatureLayer
    raise e
RuntimeError: Object: Error in executing tool

I think I am using the wrong SearchCursor or I am using it wrong. Help please.


texts:
I  wonder if it's possible to add an plugin button (QAction) to the print composer toolbar.
I only know options to add the plugin-buttons to the plugin,vector or database menu or to own toolbars in the main window of QGIS.


texts: PYQGIS, when we search based on the distance between two points, is the measurement in meters or something else?

texts:
I'm new in QGIS and in GIS.
I want to calculate distance between two points in meters using WGS 84 coordinates. I'm working with QGIS 2.0.1 Dufour and Python 2.7

The code is the following (based on http://docs.qgis.org/testing/en/docs/pyqgis_developer_cookbook/geometry.html?highlight=measureline#geometry-predicates-and-operations):

point1 = QgsPoint(-46.443077,-67.51561)
point2 = QgsPoint(-46.4446,-67.512778)

#Create a measure object
distance = QgsDistanceArea()

#Measure the distance
m = distance.measureLine(point1, point2)

Using these coordinates I'm getting 0.00321554863126
And in Google Maps it is approximately 290m.

Why is it not the same distance, or is the first distance not meters?


texts: What's the easiest way to move tile-based services?

texts:
We're still running ArcGIS 9.3.1 for a short while longer, and have been tiling more of our data, which is great until we have to deploy it. Deploying currently consists of transferring the entirety of the tileset to the production machines, which takes about half of forever, especially on a deadline. We're currently using tar.gz, but would welcome a better approach. Tiling in place is not an option as there is not a "spare" server at the production facility.

Edit:
I should mention that the destination server is a long way away, and behind a firewall, meaning that our access to it is limited, and mapping a drive or a share is Right Out. I'd be tempted to use external drives and FedEx, but that would require a physical person at the destination, which tends not to be the case for months at a time.


texts:
The data I'm interested in is found here and can be opened directly in your browser:

http://maps.greenmin.gl/geusmap/ows/32624.jsp?nocache=nocache&whoami=[email]&LAYERS=grl_geus_500k_geology_map&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&FORMAT=image%2Fpng&TRANSPARENT=true&STYLES=&SRS=EPSG%3A32624&BBOX=499980.0,7190220.0,609780.0,7400040.0&WIDTH=1500&HEIGHT=3000

Then I'm using this code in my terminal to output it to a GTiff with its reference system intact, which is based on information found here.

gdal_translate -of GTiff -co "TFW=YES" "http://maps.greenmin.gl/geusmap/ows/32624.jsp?nocache=nocache&whoami=130.225.0.251&LAYERS=grl_geus_500k_geology_map&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&FORMAT=image%2Fpng&TRANSPARENT=true&STYLES=&SRS=EPSG%3A32624&BBOX=499980.0,7190220.0,609780.0,7400040.0&WIDTH=1500&HEIGHT=3000" output_geo.tif

I then get an error

Input file size is 561895206, 1073741824
0ERROR 3: Free disk space available is 745445371904 bytes, whereas 2413321533549182976 are at least necessary. You can disable this check by defining the CHECK_DISK_FREE_SPACE configuration option to FALSE.

Any suggestions?


texts: Problem with SRTM 30 m data while creating watershed

texts: I have used SRTM 90 m in many of my hydrology projects for creating watershed and so far it has proven its worthiness. But, I haven't had any luck with its newer and high resolution SRTM 30 m prodigy. I  have attached an image comparing flow accumulation maps created from both these sources for same location. As can be seen in the image,  the flow accumulation model from 90 m is consistent and gives results as opposed to 30m in which there is discontinuity in the flow network. I was wondering if there is some workaround for this situation.


texts: How to clip a single layer of overlapping polygons in QGIS

texts:
I have a single layer of polygons. Some of the polygons overlap.

I wish to make a new layer of just the parts of the polygons that overlap.

What is the best way to do this in QGIS seeing it's all in one layer?

I tried the clip and intersect in the geo-processing tools, but they expect two layers.


texts: Symbology add field going to all other values

texts:
I am working with a complex symbology in ArcGIS Pro using two fields. The fields are for a cities parcels that have zoning and general plan designations. There are fields that have a null value in them and I have them included in the symbology:

I want this behavior but what is perplexing is that an attribute field is not showing up in the symbology but is instead going strait to the 'all other values' field.

There is nothing under the 'select values to add' area.

The fields in question are both 'text' and have a value (the two fields are shown in yellow highlighter)

What I would expect is the fields to appear as | Null, A | under the symbology but they just aren't and this has been frustrating me for a good long while.

I have tried to make a new field(s) and copy the values over etc to no effect.

*Edit: Another anomaly is that the 'count' field says that there are zero of the 'all other values' symbology but it accurately recolors what would be the | Null, A | parcels. I've still yet to solve this problem.


texts:
I have a table with points that has lon/lat, i am adding postgis to my app. I am doing some tests and I found type POINT and type GEOMETRY.  I want to convert my columns lon and lat to a point.  For what i was reading it made more sense to add column of type POINT, but most examples I have seeing use geometry to store points.

@Column(columnDefinition="Point")
    private Point locacion;

Should I use geometry or geography instead of points.  I am planning on using this points to find other points within a certain radius of meters


texts:
I try to see a point layer that been created in QGIS 2.12 in Google Earth and the problem is that I can't see the point in there - it only located me in the right frame with out the point layer.

I tried to turn on all the layers but had no success:


texts:
I have configured a date/time field with the default value "now ()" on a shapefile.
When I create an object in QGIS, the datetime is correct. In QField, there is a 2 hour lag. However, the tablet is in the correct time zone.
I don't understand where I can configure that. Is it a bug?


texts:
Is there any particular reason why I cannot clip in QGIS (2.8)?

I have a Vectorgrid which I created based on another shape - namely a country.
Now I want to clip the vectorgrid to the size of the country.
But, whenever I use the "clip" tool in QGis it gives me an empty output.
It worked just fine before and I do not know what I made differently there...

Here the two files I try to clip.
http://www.file-upload.net/download-10442231/shapes.zip.html

Other Vectortools do not work neither... Basically I have the same problem than this:
Clip in QGIS 2.4.0 resulting in empty output files

Here are some additional information:


texts:
I'm trying to change dtype of a .jp2 Sentinel 2 image, with this code:

    with rasterio.open(tmp_file) as src:
    meta = src.meta
    meta.update(dtype=rasterio.uint16)
    array = src.read()
with rasterio.open(filename, 'w', **meta) as dst:
    dst.write_band(array.astype(rasterio.uint16), 1)

But once is stored, if I check the values of the array they have changed.

Original array:

[[[10 10 10 ...  6  6  6]
  [10 10 10 ...  6  6  6]
  [10 10 10 ...  6  6  6]
  ...
  [10 10 10 ...  0  0  0]
  [10 10 10 ...  0  0  0]
  [ 4 10 10 ...  0  0  0]]]

Array when open new image: [[[21  8  8 ...  3 16 11]
  [ 8  0  0 ...  0  4  0]
  [ 8  0  0 ...  0  3  0]
  ...
  [ 7  0  0 ...  0  0  0]
  [20  7  7 ...  0  9  4]
  [ 9  3  1 ...  0  4  0]]]

And as you can see the transformation is not even consistent.

If I do the same but storing the image as a GTiff, then the values are preserved.

Do you why? Is there any way to avoid this?


texts: Running the Graphic Modeler from the Python Console in QGIS

texts: How can I run the graphic modeler from the python console in QGIS?


texts: Converting a point file, into a smoother set?

texts:
I am currently attempting to show how a coastal area was flooded from the 2010 Mentawai Ocean tsunami on the Pagai Islands. I currently have a DEM and a DEM which I have sourced from a point data set, and I am attempting to show how the coast was flooded using these data points. However, I cannot seem to find out how to show the area which was flooded using the specified water heights. To provide some information I have attached the image here.

.

I want to be able to 'smooth out' the point data sets, so that I can explicitly show which areas were flooded. I also have a horizontal inundation data set, so I can depict how far in land the tsunami water inundated the land. However, I am unaware as to how I can combine these data sets to show exactly how far flooding went and how tall flooding was.


texts: Why does map click listener not work with touch?

texts:
I deleted this question which I posted a few days ago, thinking that I would be able to solve it. But after struggling for some time, I think I really need some help. The problem is this: I am testing my map application on a touch device. There are a number of click listeners registered on the map and markers. I would expect these listeners to be able to respond to touch events that consist of a tap which is equivalent to a mouse click.

Some listeners work; some don't. Example, the base layer switcher button on the top right corner of the map respond to the click. However, the map and the markers don't. While examining the Openlayers javascript debug file. I notice that click handlers are assigned to the map, map navigation control and buttons. The buttons respond to the touch click event perfectly; the map is able to zoom in on double tap; which leads me to conclude that something must be wrong with the navigation control.

There is an awkward solution that suggested doing this:

// create click control class
OpenLayers.Control.Click = OpenLayers.Class(OpenLayers.Control, {
 defaultHandlerOptions: {
  'single': true,
  'double': false,
  'pixelTolerance': 0,
  'stopSingle': false,
  'stopDouble': false
 },
 initialize: function(options) {
  this.handlerOptions = OpenLayers.Util.extend(
   {}, this.defaultHandlerOptions
  );
  OpenLayers.Control.prototype.initialize.apply(
   this, arguments
  );
  this.handler = new OpenLayers.Handler.Click(
   this, {
    'click': this.trigger
   }, this.handlerOptions
  );
 }
});

// add click control
var clickControl = new OpenLayers.Control.Click({
 trigger: function(e) {
// do whatever
 }
});
map.addControl(clickControl);
clickControl.activate();

There are two related questions here and here. But till now, there seems to be no one able to point out what went wrong with the navigation control and how to solve this.

I spent a few days pouring through the code trying to decipher what is going on, but no luck so far. To me, the above solution causes a noticeable delay in the clicks which affected other function calls. Therefore, I am posting this to see if anyone can help. Thanks.


texts: How to avoid getting "float division by zero" error message when using hierarchical clustering in the scipy plugin?

texts:
I am trying to use the hierarchical clustering function in the plugin scipy point clustering in order to group individuals that are close to each into populations.  Here is a screenshot of my selections:

But every time I try, I get the following error message in the log:

  Algorithm  Hierarchical Clustering starting... Extracting geometries

  from the input layer

  float division by zero See log for details

I checked for duplicates in locations (perhaps there was a problem dividing by a zero distance).
I also made sure that the points are not multipart by applying the function : "split multipart features to single part" from digitizing toolbar.

Here is a sample of my shapefile:

>   wkt_geom    name    lat long
>   Point (286556 5028628)  2501    5028628 286556
>   Point (286479 5028418)  2502    5028418 286479
>   Point (286507 5028404)  2503    5028404 286507
>   Point (286516 5028412)  2504    5028413 286516
>   Point (286519 5028417)  2505    5028417 286519
>   Point (286517 5028428)  2506    5028428 286517
>   Point (286518 5028433)  2509    5028433 286518
>   Point (286510 5028426)  2510    5028426 286510
>   Point (286512 5028427)  2511    5028427 286512
>   Point (286503 5028423)  2512    5028423 286503
>   Point (286501 5028429)  2513    5028429 286501
>   Point (286495 5028419)  2514    5028419 286495
>   Point (286509 5028397)  2517    5028397 286509
>   Point (286513 5028392)  2518    5028392 286513


texts: How can I match the cells of two rasters?

texts:
I have many raster files each one containing information that is different from the others (e.g. one contains information on rain, the other on altitude, and another on land use). The cells in each raster file may or may not line up with the cells of the other raster files.

My goal is the following: I need to create a text file where each row represents a cell and, along the columns for that row, contains the information in that cell on rain, altitude, land use, etc.

What workflow would you suggests to "merge" these raster files together so that their cells line up and contain all the information available across all of the rasters?

I am using ArcGIS 10.1


texts: Dijkstra's own implementation using PostGres & PostGIS

texts:
I'm currently working on a travel routing application and would like to implement my own version of Dijsktras or A*.

I've been looking at PostGres & PostGIS for storing OSM data, added using tools such as OSM2PO.

I know there are tools such as pgrouting that can be used, but is there any resources/guides on how one would avoid pgrouting and create a custom route algorithm?

The end goal is to show a shortest path route from A -> B (possible waypoints) on a map such as OpenLayers.


texts: Near/Spatial Join picking wrong polyline segment as closest?

texts:
I have a set of points near a street centerline. I need to find the perpendicular point on the street centerline, and use their x y instead of the GPS information of the original points.

I tried Near tool and spatial join in Arcgis 10.1. They returned the same result of a segment which is far away from the right one. For example, for point (60.0001, 40.0001), I could see from the street centerline attribute table that there's a segment start from (60.0002, 40.0002). Yet the Near tool and Spatial Join tool return the same wrong result, which  a segment start from(60.002, 40.002) This could be a neighborhood away.

Anyone has any idea? Any other intuitive tools I could use?


texts: Add GPS tracks directly to a PostGIS table

texts:
Disclaimer:
I am a total GIS newb and I have a use case I have not seen elsewhere.

Background:
I drive country roads looking for old building, cemeteries, school houses, etc and I record the locations and document them for a website.  I want to use QGIS to keep track of the places in their categories (school, church, cemetery, etc) and also track the roads I have driven as a way to keep duplication of work to a minimum.

I have successfully (I think) setup PostGIS, imported all my items (places + tracks from GPS) from a KML and setup a bluetooth GPS with a small netbook.

I would like to use the built in GPS track maker to record where I have driven and append these lines to a PostGIS table which contains a list of other lines.

I have seen examples of exporting the data to a GPX or SHP file but according to this post ( How to add features to a PostGIS layer using Quantum GIS? ) I should be able to just the line to PostGIS.

According to the documentation I cannot add a GPS Track directly to the PostGIS DB, I must create a new layer and record it there then export / import the data to get it into the DB.  I have had no luck with this up to this point.  I get errors that there are missing columns, once those are added I get other errors about constraints.

How do I accomplish this?   I feel like I'm very close.


texts: Getting attribute table from WFS data source

texts:
I am attempting to get some WFS 2.0 INSPIRE data from the Finnish meteorological institute into a table format of some kind. Is there any software which can easily grab the attribute data for the WFS features I request? Currently, I am just running the queries from the browser but I am not very familiar with XML so I am having difficulties parsing the resulting document (I'm using R) into a useable format. I need a significant amount of historical information which generally precludes doing this in excel, although I may have to do that if I can't figure out another way to do it. I've also tried to grab the geo features directly in QGIS with the WFS 2.0 plugin, but always seem to get a 'bad request' server response when attempting to execute a stored query. Is there a simple way to get the data from the XML document into a tabular format, or a way to parametrize my own query to get the data I need and nothing else?

If you are interested in trying the service, you need to register an API key which is free.

http://en.ilmatieteenlaitos.fi/open-data-manual-fmi-wfs-services


texts: Control level of detail of OpenStreetMap/ESRI map using QGIS print composer

texts:
I am using QGIS 3.4 and tried to print a map covering a fixed region which has a fixed lat lon boundary. The print composer gives me too many details whereas I just want the main highways.

I am wondering if there a way to control the level of detail shown in the print composer so that the certain level of detail will be shown no matter how much I zoom in? For example, I want to print the region covered in Figure 2 with the level of detail in Figure 1.


texts:
Although there are some posts on this matter, there is no answer in any one of them. This is why I am asking it again.

One post I found was https://gis.stackexchange.com/questions/68032/raster-layer-invalid

I read information from the following link: https://hub.qgis.org/wiki/17/Arcgis_rest

I used the command:

gdal_translate "http://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer?f=json&pretty=true" s.xml -of WMS

And it generated the file successfully. However, when I try to open the file and assuming the provider is WMS, the code report layer is invalid.

The code I used is:

    file = QFileDialog.getOpenFileName(self,
               "Open WMS", ".", "WMS (*.xml)")
    fileInfo = QFileInfo(file)
    # Add the layer
    layer = QgsRasterLayer(file, fileInfo.fileName(),"wms")

    if not layer.isValid():
        print "Failed to load."
        return

I just choose the file from the dialog box.
I also tried the other command:

qgis.utils.iface.addRasterLayer("http://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer?f=json&pretty=true","raster")

by using the following code:

layer = QgsRasterLayer("http://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer?f=json&pretty=true", "layer")

if not layer.isValid():
    print "Failed to load."
    return

It also report "Failed to load". The original command can be run successfully in QGIS python command line. Also, if I try to enter the code in python console, the layer.isValid() would return true. It is just not working in standalone script.


texts: Transforming GPS tracks to line segments, and retaining start_time and end_time?

texts:
I have a long GPS track representing a route of boat survey, which looks like:

  Id    X   Y   time
  1 274604  2778803     2016/5/5 08:31:34
  2 274589  2778804     2016/5/5 08:31:49
  3 274564  2778811     2016/5/5 08:32:04
  4 274538  2778818     2016/5/5 08:32:19
  5 274511  2778822     2016/5/5 08:32:34
   ...

For analysis, I transformed this point layer into polyline, and then split the line into 100 equal segments.

How can I get the information of start time and end time of each segment?

Something like:

  SegmentId time_start time_end
  1 2016/5/5 08:31:34   2016/5/5 08:50:00
  2 2016/5/5 08:50:00   2016/5/5 09:12:25
  3 2016/5/5 09:12:25   2016/5/5 09:35:12
  ...

I approached it by Spatial Join, but it's erroneous, because the route has some overlapped parts (same coordinates but different time).


texts: Reading raster files by block with rasterio?

texts:
I find rasterio much easier to use that plain ogr/gdal, so I'd prefer using it. However, I need to process large rasters, so reading rasters by block as explained in the Python GDAL course is important for me.

Does rasterio support this or should I handle the large, tilable rasters with GDAL instead?


texts: Overpass API JSON output cannot be read by QGIS

texts:
My steps:
1.) Query Overpass API and download the JSON output by this URL
https://www.overpass-api.de/api/interpreter?data=[out:json];node[route=bus];way[route=bus];relation[route=bus](52.5,13.3,52.5,13.4);out;

2.) Loading this JSON into QGIS
Issue: QGIS 3.x says, that this JSON cannot be read.
My ideal approach is to not use QuickOSM or other tools in between.
Do I miss any ETL from JSON output to QGIS?


texts: Setting a custom projection in ArcGIS, from data provided by the Government of Mauritius?

texts:
I am trying to save a custom projection in ArcGIS. It is a local projection used by the Land Unit from the Government of Mauritius. Below are the information provided by the government:

My question is, what is the Latitude of Origin in the "New Projected Coordinate System" console in ArcGIS (shown below). Is it the latitude of the "Origin: Verdun PTo0391" ?


texts: PyQGIS QgsVectorDataProvider unexpected type in argument in method changeGeometryValues

texts:
When I try to modify a features geometry in a vector layer (linestrings) resulting from a processing.runalg() operation like

clip = processing.runalg('qgis:clip', ...)
layer = QgsVectorLayer(clip['OUTPUT'], 'TheLayerName', 'ogr')

# add layer to map registry, load named style, create spatial index

layer_dp = layer.dataProvider()

# assume feat is a feature of mem_layer
# assume pt some QgsPoint on feat.geomatry(), but not a vertex

dist, p2, to_id = feat.geometry().closestSegmentWithContext(pt)

# obviously, insertVertex() does not accept QgsPoint as first argument
newgeom = feat.geometry().insertVertex(pt.x(), pt.y(), to_id)

layer_dp.changeGeometryValues({feat.id() : newgeom})

as decribed in the PyQGIS Programming Cookbook 2.18, Section 5.4.3, p.21 (https://docs.qgis.org/2.18/pdf/en/QGIS-2.18-PyQGISDeveloperCookbook-en.pdf)

I get the following error:

  QgsVectorDataProvider.changeGeometryValues(unknown-type): argument 1
  has unexpected type 'dict'

What can be possibly gone wrong here?


texts: how can I get the vertex list from a shapefile line using python or R?

texts: I need to extract the vertex list from a line inside a shapefile. I want to know which would be a better option between python and R.


texts:
I have a dataset that shows the borehole yield for fractured aquifers in South Africa. Using the "yield_ls" column, my CARTO VL legend looks like this (notice how the 4th item should be the first).

const aquiferViz = new carto.Viz(`
  color: ramp($yield_ls, teal)
  strokeColor: rgba(0,0,0,0.3)
  strokeWidth: 1
`);

I noticed that each legend item is assigned an ID column in the database ("sayld_id") starting from 6 to 10. Therefore if this database is ordered by the sayld_id ASC, then the yield_ls features go into a logical order. I tested this using the following SQL command in my CARTO dataset:

SELECT *
FROM fractured
ORDER BY sayld_id ASC

I then tried to do this when adding my the layer to a CARTO VL map using:

const fracturedSource = new carto.source.SQL(`
  SELECT *
    FROM fractured
    ORDER BY sayld_id ASC
  `);

But the legend remains the same as before. Is there a better way to correctly arrange the legend, or possibly move an item up or down in the legend order?


texts: Exporting Multiple Rasters to a Table: working with large dataset

texts:
I have 20 rasters that cover the same geographic extent, are in the same projection, have the same cell size etc. Instead of a stack of rasters, I want a table where the columns are the different rasters and the rows are the pixels. I have tried multiple tools in ArcGIS, such as Extract Multi Values to Points and Sample, as well as using R. The issue is my rasters are very large, about 60 million pixels each. Those tools either don't work correctly and output incorrect values (Extract Multi Values to Points) or fail after a short time (Sample). I have tested them out with much smaller sub-samples (i.e. 100 points) and they work fine in that case. My R attempts have been met with errors due to lack of memory. Any tips or other ways I can approach this problem? I am willing to try different methods! I have read similar questions here and elsewhere in the past but it appears the large size of my dataset is leading to additional issues.

The goal is to have a table like this:

R1 - R2 - R3 etc.

P1
P2
P3
etc.

That I can bring into R and work with further.

EDIT: My R code... very simple at the moment (also down to 12 rasters).

Rasters <- stack(NDVI_C, Landcover, Elevation, Slope, Aspect, TPI, TWI, Lake_50, River_250, MajRiver, FTE, Precip, Temp)
d <- values(Rasters)
test <- values(NDVI_C)

A couple issues: memory and dealing with NA values. Want final table to only include pixel stacks that have no NA values (i.e. every Raster has a value for that pixel)


texts: Google Maps JS API v3 - Polygon formation

texts:
So my big question here is how does the Google Maps API determine the "inner" versus the "outer" for a filled polygon?

Playing around with the polygon creator @ http://www.the-di-lab.com/polygon/ there seems to be a preference to:

wrap the north pole;
take the "smallest" polygon possible?

The big reason I want to know is to assist with answering the question "is point (x,y) in multipolygon z?"

Thoughts?


texts: QGIS3 - how to define/change field type for non-spatial data table?


texts:
I have series of data in my GDB

I am running the program under VBA and what I want to do is using the data by the tool Kriging.

This is what I thought:

First, Kriging needs data that must be a point shapefile and I am getting the data out of my GDB (one point in one table but there is a lot of data)

I get the data I need and make it in a shapefile
I use the Kriging tool to calculate the shapefile
I use the raster to point tool to change the raster Kriging make to a point layer
I get the point data I want!

Here is my question:

When I get the data from the GDB (using sSQL), how to change it to a shapefile in VBA?

I have only found the page to display XY data using VBA: http://edndoc.esri.com/arcobjects/9.0/samples/tables/create_a_layer_based_on_xy_data.htm but it's can't save as shapefile!

How to save it automatically?

Is there any wrong with my steps?

If need more information please tell me!

@artwork21

An error occurred in the line

gp.FeatureClassToShapefile_conversion "'C:\Temp\New File Geodatabase.gdb\point'", Temp__2_

says

automation error

unspecified error

how to fix it?

Call tools from ArcToolBox in VBA form

in this post about my problem I fix it like follow

gp.Kriging "C:\temp\123.shp", "F30", "C:\AA\Kriging_123", "Spherical 249.101224", "250"

but this time no matter what I try it's doesn't work...


texts: Labels not displaying correctly in the centroid of the polygon in GeoServer SLD

texts:
When creating labels using centroids of a polygon layer, I'm facing the problem that some of them are missing. Eg:

Ok, I can make it reappear by using "maxDisplacement" vendor option, but, using this option labels appear displaced. Some to the right, some to the left.

How can I render all the labels propperly in the centroid of the polygons?

Here is the SLD code I'm using:

<?xml version="1.0" encoding="UTF-8"?>
<sld:UserStyle xmlns="http://www.opengis.net/sld" xmlns:sld="http://www.opengis.net/sld" xmlns:ogc="http://www.opengis.net/ogc" xmlns:gml="http://www.opengis.net/gml">
  <sld:Title/>
  <FeatureTypeStyle>
     <Rule>
       <PolygonSymbolizer>
         <Fill>
           <CssParameter name="fill">#40FF40</CssParameter>
         </Fill>
         <Stroke>
           <CssParameter name="stroke">#FFFFFF</CssParameter>
           <CssParameter name="stroke-width">2</CssParameter>
         </Stroke>
       </PolygonSymbolizer>
       <TextSymbolizer>
        <sld:Geometry>
          <ogc:Function name="centroid">
            <ogc:PropertyName>geom</ogc:PropertyName>
          </ogc:Function>
        </sld:Geometry>
         <Label>
           <ogc:PropertyName>cod_25000</ogc:PropertyName>
         </Label>
         <Font>
           <CssParameter name="font-family">Arial</CssParameter>
           <CssParameter name="font-size">11</CssParameter>
           <CssParameter name="font-style">normal</CssParameter>
           <CssParameter name="font-weight">bold</CssParameter>
         </Font>
         <LabelPlacement>
           <PointPlacement>
             <AnchorPoint>
               <AnchorPointX>0.5</AnchorPointX>
               <AnchorPointY>0.5</AnchorPointY>
             </AnchorPoint>
           </PointPlacement>
         </LabelPlacement>
         <Fill>
           <CssParameter name="fill">#000000</CssParameter>
         </Fill>
         <VendorOption name="autoWrap">100</VendorOption>
         <VendorOption name="maxDisplacement">50</VendorOption>
         <VendorOption name="repeat">1</VendorOption>
       </TextSymbolizer>
     </Rule>
   </FeatureTypeStyle>
</sld:UserStyle>


texts: Problem with variable for arcpy.select_analysis whereclause

texts:
I have a problem with the whereclause on python for arcpy.Select_analysis.

I try to use a variable for condition:

for row_stat in rows_stat:               % for calculated the variable"maxi"
    maxi= row_stat.MAX_GRID_CODE         % for calculated the variable"maxi"
del rows_stat                            % for calculated the variable"maxi"
where = '"GRID_CODE" = maxi'             % whereclause
arcpy.Select_analysis("point_acc.shp", "point_max.shp", where)

Ofcourse, I get an error message telling me that 'maxi' can't be found in "point_acc.shp" but I don't want 'maxi' but the number in the variable maxi!

Who can help me write this?


texts: Forcing ArcMap to Ignore Double Quotation Mark In CSV Text Field When Exporting to Shapefile?

texts:
In ArcMap 10.2 I'm following this workflow:

Add .csv to map.
Right click .csv, --> "Display XY Data"
Right click XY Event Layer, --> "Data" --> "Export Data" --> save it as a shapefile.

There are 1232 rows in my .csv but the resulting shapefile has only the first 235 rows.  The last row in my new shapefile (row 235) has one double quotation mark (") in the fifth (text) column.

This screen grab shows how the double quotation mark should appear in its cell.

This screen grab shows the result of bringing that row into a shapefile in ArcMap.  Notice the quotation mark is absent but the lat, lon and subsequent field values appear in the same text cell, while the subsequent lat, lon, and so on fields are empty.

I found out that you can escape a double quotation mark by adding a second right after it, like this... (""), so I made this change manually in Excel.  This was the solution.  Next time I followed the above workflow I had all 1232 rows in my new shapefile.

How can I achieve this escaping of the double quotation mark in ArcMap rather than doing it manually in Excel?

A related side problem: Although setting the " manually to a "" in Excel allows all 1232 rows to be converted in ArcMap, the cell in question in the resulting shapefile has "" instead of where it should only have ".  It look like this:


texts: GeoJSON - Obtain linestrings inside a given bounding box

texts:
I have two layers in an Argentinian railroad-related project: one is for train stations and the other is for railroad lines. The first one has about 1100 points and the second, about 250 finely-detailed linestrings. I use leaflet.js as my framework, and a Bing layer and a background satellite imagery.

Right now, I have a web application in which I charge straight into memory both GeoJSON files, and it works just OK, except when I need to zoom or pan to a specific feature: it works really slow. I tried removing almost all the features, leaving just two or three linestrings and a couple dozen stations, and it works fast enough to be usable. So, I guess I need to refactor my code so every time I pan or zoom the map, there's an Ajax query to a PHP script that filters my layers using a bounding box query.

It is trivial to filter train stations because they are just points. My problem arises when I need to filter railroad lines: as I said, they are finely drawn, so each one of them has hundreds if not thousands of points. I tried a "foreach" loop with a coordinate query inside, but it takes forever to iterate and filter every coordinate that falls outside the bounding box. So, my question is: is there any way to make this more effective, or at least less CPU-intensive? I need to return a geojson object containing trimmed lines that fit nicely into the provided bounding box.

Edit: please bear with me, as I am an absolute amateur on this subject.

Edit 2: to summarize what I'm trying to achieve: I need to filter a GeoJSON object containing linestrings to return a sub-set of linestrings that fit inside a given bounding box, trimming those linestrings that are bigger than it.

Edit 3: forgot to mention that simplified geometries aren't a valid solution for me, as I need the current level of detail for railroad routes.


texts: What solution is there to project google/bing layers to another EPSG (e.g. 3844)?

texts:
I am using Geoserver and Openlayers. The Openlayers map has all the layers in EPSG:3844.
I need to overlay Google/Bing layers.
I am currently trying to do it with proj4js,  mapproxy.org as openlayers doesn't have reproject support for WMS.

Any solutions?


texts:
Problem: I have a bunch of raster data (DTED) and I want to use this to make elevation queries using PostGIS (if it makes sense). E.g. I want to ask the question "what is the elevation at (lat1,lon1), (lat2, lon2)?" and the database would return the interpolated elevations based on the nearest neighbors where I have data.

Idea 1: Should I load the raster data into the data base in its raw form and make queries on that? Immediate downside sounds like the rasters are pretty big and maybe querying would be slower?
Idea 2: should I preprocess the raster data and load a bunch of points into the database (i.e. the database would have columns latitude, longitude, and altitude for all the measurements in the DTED data).

An example is show below.The intersection of the grid represents points that I have, the red stars represent points that I am querying for. So ideally, I would get a bilinear interpolation (or some other method) between my nearest neighbors. In the comments below, it was mentioned that the function http://postgis.net/docs/manual-dev/RT_ST_InvDistWeight4ma.html would help with this calculation. However, it is not clear to me how to use this function to get the answer I am looking for.


texts: QGIS 2.18.x - GRASS 7.2.x installation in Linux deepin 15.4: which repository to use?

texts: In the default linux deepin repository there is the QGIS LTS 2.14.x version. Which repository could I use (if possible) to install the stable version 2.18.x?


texts:
I am trying to add new  text ID's to a field using a basic counter function within the field calculator. I have two questions here 1. does the field calculator always have to use a function or can it use straight code that's not in a function? 2. I keep running into error 00539 which I believe is a syntax error. When I run this in a regular IDE it works out fine but when i throw it into the field calculator I always get that same error.

How do i fix this?

Below is the code I've written.

def countfunction(input):
count=0
for N in range (1600,2360,5):
    if N >=0:
        print "TBL01_01172017_"+format(N)+"00.JPG"
        count =count +1
input = countfunction(1)
print input


texts:
How to fix the raster.asc to work with Maxent?

Running Maxent with raster from WorldClim and others as Altitude, LandCover, etc. I have problems because shows up one error (have different geographic dimensions). I have been trying to fix that in several ways, probing that the spatial reference is the same, the last thing I did was to clip the raster with the same extent introducing manually the coordinates. I am not sure if is neccesary exactly the same coordinates in the four points for the whole set of rasters. Maxent works with the 19 worldclim variables but I need to include altitude, landcover, etc, but I have this problem.


texts: [Error0125] error when create new service with AXL file in ArcIMS 9.3.1[Help]

texts:
i have problem to create a new service with AXL file in ArcIMS 9.3.1,when i do create a new service error message appear (error0125). so can any one help me how to solve this problem please.

this is pic shown that error message:

The picture is also posted here if imgur.com is blocked from you site.

AXL file: [URL]http://www.mediafire.com/?3dp2sbg5yx1k4ac[/URL]

<?xml version="1.0" encoding="UTF-8"?>

<ARCXML version="1.1">
  <CONFIG>
    <ENVIRONMENT>
      <LOCALE country="US" language="en" variant="" />
      <UIFONT color="0,0,0" name="SansSerif" size="12" style="regular" />
      <SCREEN dpi="96" />
    </ENVIRONMENT>
    <MAP>
      <PROPERTIES>
    <ENVELOPE minx="-94.04178542666283" miny="28.939655671907317" maxx="-89.02180295398398" maxy="33.02342240774887" name="Initial_Extent" />
    <MAPUNITS units="decimal_degrees" />
      </PROPERTIES>
      <WORKSPACES>
    <SHAPEWORKSPACE name="shp_ws-0" directory="C:\Documents and Settings\Mustafa\Local Settings\Temp\LearnArcIMS9\VirtualCampus\LearnArcIMS9\Start" />
      </WORKSPACES>
      <LAYER type="featureclass" name="laparishes" visible="true" id="0">
    <DATASET name="laparishes" type="polygon" workspace="shp_ws-0" />
    <SIMPLERENDERER>
      <SIMPLEPOLYGONSYMBOL boundarytransparency="1.0" filltransparency="1.0" fillcolor="127,127,27" boundarycaptype="round" />
    </SIMPLERENDERER>
      </LAYER>
    </MAP>
  </CONFIG>
</ARCXML>


texts:
I'm trying to create a 3D scatterplot using an las file that contains xyz and return values. I want to keep all returns but would like to subset the cloud. I've learned how to subset using laspy so no issue there. The trouble comes when I'm trying to load the las to matplotlib.

How do I do I load the las to matplotlib? Is there a better way to view it? I'm doing it this way so I have scales around the outside of the las and 3d view capability.

The code below, loads two points that are at the extents of the file, upper left and lower right. I want to load those and all the points in between.

if I try to use the line

x,y,z = points

I get the error

Traceback (most recent call last):
  File "<stdin>", line 1 in <module>
ValueError: too many values to unpack

Code:

from liblas import file
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
fig = plt.figure()
las_file = r"path/to/las"
f = file.File(las_file, mode='r')
num_points = int(f.__len__())
PointsXYZ = np.empty(shape=(num_points, 3))
counter = 0
for p in f:
    newrow = [p.x, p.y, p.z]
    PointsXYZ[counter] = newrow
    counter += 1

ax = fig.add_subplot(111, projection= '3d')

ax.set_axis_bgcolor('black')
ax.xaxis.label.set_color('white')
ax.yaxis.label.set_color('white')
ax.zaxis.label.set_color('white')
ax.tick_params(axis='x',colors='white')
ax.tick_params(axis='y',colors='white')
ax.tick_params(axis='z',colors='white')

listmin = np.amin(PointsXYZ, axis=0)
listmax = np.amax(PointsXYZ, axis=0)
zmin = listmin[2]
xmin = listmin[1]
ymin = listmin[0]
zmax = listmax[2]
xmax = listmax[1]
ymax = listmax[0]

xs = (counter, xmin, xmax)
ys = (counter, ymin, ymax)
zs = (counter, zmin, zmax)

ax.scatter(xs, ys, zs, c='r', marker='o')
plt.show


texts: Eurostat shapefiles plotting with errors

texts:
I am trying to plot interactive maps at the NUTS3 level. Initially I was using shapefiles that I constructed myself for the data in QGIS, however, when I tried to use them in R they rendered poorly with lines or even entire polygons being warped out of place/shape. I then tried using the shapefiles provided by Eurostat directly (which my shapefiles were large based off of) using the Eurostat package, however, when I plot the map there are still lines/polygons that are being warped out of position (different ones from those of my original shapefiles, however).
I have what I believe should be a reproduceable example showcasing the issue below:
# Get relevant NUTS3 shapefile from Eurostat
NUTS3EU <- get_eurostat_geospatial(output_class = "df", resolution = "60")

# Generate random number vector
x <- rnorm(n = nrow(NUTS3EU), mean = 1)

# Merge random number vector into data as new variable 'RandVar'
NUTS3EU$RandVar <- x

# Subset data down to just Germany NUTS3 level
deuNUTS <- NUTS3EU[(NUTS3EU$CNTR_CODE == "DE"),]
deuNUTS <- deuNUTS[(deuNUTS$LEVL_CODE == 3),]

# Generate ggplot object
userMap2 <- ggplot(deuNUTS, aes(long, lat, group = NUTS_ID, fill = RandVar, data_id = NUTS_ID, tooltip = id)) + geom_polygon_interactive(colour = "white", size = 0.3)

# Generate Girafe interactive plot object
ggiraph(ggobj = userMap2, tooltip_offy = -20, zoom_max = 3)

# Output should be map of Germany with small number of lines/polygons being warped, may need to zoom in to truly appreciate the issue

I understand that the 'order' variable is important and have made sure it is correctly sorted, but the issue persists. How can I make the plot geographically correct? If I need to use a different source for the shapefile that is fine as long as it is accurate.


texts: Change label position for a specific feature

texts: I try to change label position to specific feature. I red  How do I choose a specific position for ArcMap labels , but i don't wont to convert the labels to annotation. Is it possible? I work with arcGIS 10.3


texts: get segments of constant length from Linestring MS SQL SERVER

texts: For Example I have simple linestring 'LINESTRING(25 50, 100 125, 150 190)' and I need to divide it to segments of length 10. In postgis i know there is a function ST_Segmentize but what about MS SQL Server?


texts: Set metadata for a shapefile (GDAL/OGR Python bindings)

texts:
I'd like to set a metadata value within a shapefile using Python bindings. I'm using the SetMetaData() / GetMetaData() methods within GDAL (gdal.org/python/osgeo.ogr.MajorObject-class.html). I've had no problem with raster files, but when I go to do the same with shapefiles, it will update the datasource while it's open, but won't save them if it's closed.

ds = ogr.GetDriverByName('ESRI Shapefile').Open('C:/Polygons.shp', 1)
print( ds.GetMetadata() )
ds.SetMetadata( {'ThisKey':'ThisValue'} )
print( ds.GetMetadata() )
ds.Destroy()
del ds

ds = ogr.GetDriverByName('ESRI Shapefile').Open('C:/Polygons.shp', 1)
print( ds.GetMetadata() )
ds.Destroy()
del ds

produces

{}
{'ThisKey': 'ThisValue'}
{}

I can't see what I'm doing wrong.


texts: How can I create an OmniSearch using OSM data?

texts:
I want to create a "search-for-anything" using the OSM data. I know the name column is the most relevant attribute, but it come empty sometimes.

What tags I must include in my query to return the maximum data possible of a search?

Example: When the user search for Peter I need to return roads, places, hospitals, etc... with this parameter in the name.

I'm using Java to query the planet_osm_* tables in a PostGIS database. I'm planning to create an Ajax autocomplete method in the HTML users interface to query the database for any word the user types.


texts: Extracting slope of aspect information using QGIS?

texts:
In What is SOA (slope of aspect)? @Spacedman explained to me what SOA means. Now I need help to understand how I could work with it and extract the data that I want.

I want to use Qgis (with Grass and Saga) and, if necessary, Python to this task. My goal is to delineate the rim of impact craters. I have been trying a bunch of methods and right now I want to try doing this by working with SOA information.

I know how to make an slope and aspect map from the DEM. But I don't know how to extract the regions of high values of SOA (crater's rim).


texts: OpenLayers 2: Can I delete vector layer or change its z-index when mouse hover on it?

texts:
I am adding two vector layers.

this.mapport.vector = new OpenLayers.Layer.Vector();
this.mapport.inside_vector = new OpenLayers.Layer.Vector();
map.addLayer(this.mapport.vector);
map.addLayer(this.mapport.inside_vector);

and I am selecting features of inside_buffer like

inside_vector.setZIndex(999);
var inside_feature = new OpenLayers.Feature.Vector(inside_geom, null, this.inside_buffer_style);
this.mapport.inside_vector.addFeatures(inside_feature);

and vector layer 1 features as

var feature = new OpenLayers.Feature.Vector(geom, null, this.buffer_style);
this.mapport.vector.addFeatures(feature);

this.request_select(feature, options, event, layers, behavior);

Initially I want to show user inside_vector layer but when user will hover over it Its z-index must be 0. Can I do it in OpenLayers 2?  till now I tried 4,5 examples on stack overflow but they are about features not about layer itself.


texts: Calculating from a location to another location (calculation must be on a line/road path)

texts:
I'm using QGIS. I'm new.

I have a CSV file, with the source point (the_geom) and destination point (the_geom).

I have a road network already imported into PostGIS with pgRouting, and I need to have the ability to calculate the distance from a source to a destination, by using the road path as the road. After calculating, it should save the distance back into the DB in a new column.

Is there any quick and easy way to solve this?

Thanks for all the help! :)


texts: Download data from webmap made with leaflet

texts: Can I download the data from a webmap made available from leaflet? I am referring to these dataset (https://www.vectorbase.org/popbio/map/), for which I wish to download the data for Europe? Is there a way to import these data into QGIS providing the link of the website?


texts: Remove Field Name from Pop Up Label in ArcGIS Online

texts:
I am trying to remove the field name in the pop up label in one of my map layers in ArcGIS Online in order to simplify the label.

In the below example, I don't want to show the words: Name, Address, ST_DIR, etc.

Screenshot:
https://i.imgur.com/16aKAbV.png


texts: How to Connect OpenLayers Map to non-GeoServer / PostGIS Setup?

texts:
I apologize in advance if my question is unclear - I am afraid I do not understand all of this quite well enough to form my question precisely, but I will give it a shot.

I recently managed to get an OpenLayers map set up and displaying a WFS layer (mostly by piecing together examples from others' questions), the data for which is stored in a PostGIS database and served with GeoServer. My next task is to figure out how to connect the map to our existing database, which is an Accumulo database. We have a Jetty server set up for communicating with the Accumulo database. To the best of my knowledge, what I want to do is to replace the WFS protocol on my vector layer with a different protocol that will communicate with this Jetty server. If we can get a request sent to the Jetty server in some easily parsed format, we already have the capability to turn this into a database query and return JSON with the results of the query.

I have looked at both the HTTP and Script protocols, and it seems to me that one of these may be what I want to use. However, Google searches provide little in the way of examples of how these protocols work.

Am I on the right track with this?


texts:
In my studies, I work with the internal dynamics of different snowpatches. In this study I wish to generate 3D models of snow patches.  I have some problems when it comes to the actual work flow of creating this 3D model:

Data I have available are GPS points (x, y, and z values) over large parts of the research area (these points are collected by using a Real Time Kinematic GPS, on the surface of the snow patches) and GPR (ground penetrating radar) data, which tells me what's under the ground.  For example, how deep is the ice, how thick are the snow layers, etc.

My thought is that I want to create a TIN from gps points using height values. The problem here is that I do not have GPS points for the whole research area, about 80% is covered.

Question 1: Are there any interpolation tools that can estimate the points where I do not have the point data, based on the points I already have?

When it comes to the actual 3D modeling, that I have a TIN of the glacier, but will now set the GPR data.

Question 2: Would it be best to create eg. three shapefiles one for ice, one for snow and one for firn? Then put different base heights on these files according to the Z value in the GPR data? Or can I do this in a simpler way?

Highly appreciate all the help I can get.

Regards,
Jorn


texts: Adding osm data to Postgres via Osm2pgsql : duplicate key value violates unique constraint

texts:
I'm trying to import osm global world data into a database created in postgreSQL but after a while of processing, I've got this error :

Reading in file: E:\planet-latest.osm\planet-latest.osm
Processing: Node(1873580k) Way(0k) Relation(0k)COPY_END for COPY planet_osm_node
s FROM STDIN;
 failed: ERREUR:  la valeur d'une cl├® dupliqu├®e rompt la contrainte unique ┬½
planet_osm_nodes_pkey ┬╗
CONTEXT:  COPY planet_osm_nodes, ligne 1758128898 : ┬½ 2147483647       68830134
7       -16699132       \N ┬╗

Error occurred, cleaning up

I read that I have to remove the slim mode but when I do that, another error appeared and tell me I don't have RAM enough.
Have you got any idea?


texts: Moving multiple polygon features simultaneously in arcgis js

texts:
Is there a way of moving two or more polygons simultaneously?

The following code enables the movement of a singular graphic:

           moveToolbar.activate(esri.toolbars.Edit.MOVE, tempMoveLayer.graphics[graphicNum]);

However I cannot seem to find a method that allows the movement of multiple selected graphics simultaneously.

EDIT

With the following code I manage to activate multiple toolbars for multiple graphics at the same time, however they still move separately rather than being moves as a whole when one of them is clicked and dragged.

//Loop to get all the layers on the map
for (var gCount = 0; gCount < app.map.graphicsLayerIds.length; gCount++) {
    //This will obtain the layer that the selected graphic is on
    var tempMoveLayer = app.map.getLayer(app.map.graphicsLayerIds[gCount]);

    //Looping through all the graphics in the selected layer
    for (var graphicNum = 0; graphicNum < tempMoveLayer.graphics.length; graphicNum++) {

        //Comparing the selected graphic to the layer graphic currently in the loop
        if (tempMoveLayer.graphics[graphicNum].attributes.OBJECTID == results.features[j].attributes.OBJECTID) {

        //Activating the toolbar for the graphic. This is activated for every graphic in the loop
        moveToolbar.activate(esri.toolbars.Edit.MOVE | esri.toolbars.Edit.SCALE | esri.toolbars.Edit.ROTATE, tempMoveLayer.graphics[graphicNum]);

                                                            }
                                                        }
                                                    }


texts: Calculate row and column for WMTS request from Lat/Long

texts:
I'm writing a request to a WMTS server that has an upper left corner at 90, -180 and spans the whole Earth. I want to get a tile containing some point (targetLat, targetLong).

Suppose we have some tile matrix with a given matrixWidth, matrixHeight, tileWidth, and tileHeight, and some upper-left corner matrxLat, matrixLon.

Naively, one could calculate the N-S distance from the upper-left corner of the matrix to the target point by taking the haversine distance from (targetLat, targetLong) to (matrixLat, targetLong), and the E-W distance from (targetLat, targetLong) to (targetLat, matrixLong), then use simple algebra to find the correct tile from there.

The problem is that the E-W distance is going to change depending on which latitude I measure it from, so I don't think this will work. How can I translate my target lat-long point into a meter offset from the matrix corner in the correct coordinates?


texts:
I have a collection of Landsat images filtered over forest regions of India over a time interval of 6 years(2013 - 2018). I have again filtered this collection such that I have one image per 4 months, that is I have filtered on a quarterly basis so that I will have 3 images per year.

Now, for a particular point, there will be 18 pixels representing it over 6 years. I have to find, out of these 18 pixels, how many satisfy the condition - NDVI < 0.25. Suppose, 10 pixels satisfy this condition, then in my output image, the pixel representing that particular point should have its value = 10. Like this I have to do for all the points in my selected region.

I have no idea as to how to do this pixel wise computation. Can somebody help me on this?

In my code, I have done till filtering quarterly wise.

// FILTERING OUT FOREST COVER IN INDIA

var treeCanopyCover = hansen.clip(india).select('treecover2000');
var greater25 = treeCanopyCover.gte(25);

// NDVI FUNCTION

function addNDVI(image) {
var ndvi = image.normalizedDifference(['B5', 'B4']).rename('NDVI');
return image.addBands(ndvi);
}

// FILTERING THE COLLECTION

var startDate = ee.Date('2013-01-01');
var endDate = ee.Date('2019-06-30');
var filtered = l8.filterBounds(india)
      .filterDate(startDate, endDate)
      .filterMetadata('CLOUD_COVER','less_than',25)
      .map(addNDVI);

var maskedCollection = filtered.map(function(image){
  var x = image.updateMask(greater25);      // filtering over forest regions
  return x;
});

// FILTERING QUARTERLY WISE

var map_m = function(y)
{
  var months = [1, 5, 9];
  var filtered_col = months.map(function(m)
{
var start = ee.Date.fromYMD(y, m, 1)
var end = start.advance(3, 'month');
var filt = maskedCollection.filterDate(start,end).mean();

return filt.set('year', y)
          .set('month', m)
          .set('system:time_start',ee.Date.fromYMD(y,m,1)) ;
  });
  return filtered_col;
};

var years = ee.List.sequence(2013, 2019);
var quarterlyimages = years.map(map_m).flatten();
var z = ee.ImageCollection(quarterlyimages);


texts: Line layer not displaying correctly in QGIS

texts:
I am having issues displaying a .shp file correctly in QGIS 2.12.3.

Sometimes when I add it to a map, it displays like the red lines in the photo below (which is how it should be). However, sometimes it displays as the gold lines in the other photo below (notice there are several lines missing in the top-center section).

I am afraid there is an issue with the geometry of the layer, but I haven't figured out how to diagnose it or fix it.

Note that the layer was created by joining together three different line layers. Originally these layers were not in the same projection, so I couldn't join them together. I thought I was able to change the layers to the same projection and join them together and thus successfully display the correct layer, but I fear I may have done something wrong.

Note that the red lines below only represents one of the original layers that was joined. So the remaining two component layers aren't always appearing as part of the combined layer as they should.


texts:
I am trying to export mapping and boundary data from QGIS into AutoCAD to use for concept drawings in South Australia. At my previous company, I would just set up the projection as WGS84, then set the project CRS to MGA Zone 53 or 54, and it would all look perfect, and it would export easily.
When I try to do that now, the maps either look utterly bizarre, or they are in totally the wrong coordinates (not just 'off' by a bit).
I would always tick the 'On The Fly' option in previous versions (using V3.10 now), and in fact I temporarily rolled back to version 2.18 purely to tick that box, but the result was the same.
I should be seeing coordinates around E:280,000 - N:6,120,000. Currently the best I can get them is to E:55,500,000 - N:-4,100,000, which is clearly wrong.
Is anyone able to shed some light on how to fix this, and to tell me how badly I understand projections?


texts: How to edit multiple layer simultaneously with Python in QGIS 3.0?

texts:
I have many layers but only want to edit several of them in QGIS 3.0

layer1 = QgsVectorLayer('...','layer1','ogr')
layer2 = QgsVectorLayer('...','layer2','ogr')
layer3 = QgsVectorLayer('...','layer3','ogr')
layer4 = QgsVectorLayer('...','layer4','ogr')
layer5 = QgsVectorLayer('...','layer5','ogr')
layer6 = QgsVectorLayer('...','layer6','ogr')
layers_edit = QgsProject.instance().mapLayersByName('layer1' and 'layer2' and 'layer3')

In every layer an area is given ('area_calc')

sa = layers.dataProvider().fieldNameIndex('area_calc')

An error occurs:

Traceback (most recent call last):
  File "C:\PROGRA~1\QGIS3~1.0\apps\Python36\lib\code.py", line 91, in runcode
    exec(code, self.locals)
  File "<input>", line 1, in <module>
AttributeError: 'list' object has no attribute 'dataProvider'

Making a list of layers and startEditing() also doesn't work.

layers.startEditing()
Traceback (most recent call last):
  File "C:\PROGRA~1\QGIS3~1.0\apps\Python36\lib\code.py", line 91, in runcode
    exec(code, self.locals)
  File "<input>", line 1, in <module>
AttributeError: 'list' object has no attribute 'startEditing'

The next step I want to do is:

layers.dataProvider().addAttributes([QgsField ('eur1_km2',Variant.Double),QgsField('eur1_pol', QVariant.Double)])
layers.updateFields()

idx=layers.dataProvider().fieldNameIndex('eur1_km2')
idx1=layers.dataProvider().fieldNameIndex('eur1_pol')

for f in layers.getFeatures():
    f[idx] = 6
    attr_sa = f.attributes()[sa]
    attr_idx = f.attributes()[idx]
    f[idx1] = (attr_sa/1000000)*attr_idx
    layers.updateFeature( f )

layers.commitChanges()

But the same error occurs. Has anyone an idea how to edit these three layers simultanously with python in QGIS 3.0?


texts: I have a binary raster (values - 0 & 1) - I want to create a shapefile (KML for use in Google Earth) of the white area only (value = 1). Is there a simple method to do this?


texts: Issue on Tracing Geometric Network with ArcPy

texts:
I am trying to run this ArcPy example from ESRI TraceGeometricNetwork example 2 as:

# Import arcpy module
import arcpy

# Local variables:
Water_Net = "C:/GeometricNetworks/Montgomery.gdb/Water/Water_Net"
Flags = "C:/GeometricNetworks/Montgomery.gdb/Water/Flags"
Water_Net_2 = "Water_Net"

# Process: Trace Geometric Network
arcpy.TraceGeometricNetwork_management(Water_Net, Water_Net_2, Flags, "FIND_PATH", "", "", "Friction_Factor", "Friction_Factor", "", "NO_TRACE_ENDS", "", "", "", "AS_IS", "Friction_Factor", "Friction_Factor", "1-10", "AS_IS")

but I am getting this error on running the code

  Message File Name Line Position  Traceback
  call C:\Program Files (x86)\PyScripter\Lib\rpyc.zip\rpyc\core\netref.py 196  syncreq
  C:\Program Files (x86)\PyScripter\Lib\rpyc.zip\rpyc\core\netref.py 71
  sync_request C:\Program Files
  (x86)\PyScripter\Lib\rpyc.zip\rpyc\core\protocol.py 431  serve
  C:\Program Files (x86)\PyScripter\Lib\rpyc.zip\rpyc\core\protocol.py
  379
  _recv C:\Program Files (x86)\PyScripter\Lib\rpyc.zip\rpyc\core\protocol.py 337  recv
  C:\Program Files (x86)\PyScripter\Lib\rpyc.zip\rpyc\core\channel.py 50
  read C:\Program Files
  (x86)\PyScripter\Lib\rpyc.zip\rpyc\core\stream.py 166
  exceptions.EOFError: [Errno 10054] An existing connection was forcibly
  closed by the remote host

Can you please let me know what is happening? and How I can fix this?


texts:
I try to run a simple code in pyqgis with no gui and I get a segmentation fault error by only invoking QGSApplication.

The code is the following one (under Ubuntu16.04 and with python 3.5)

from qgis.core import *

# supply path to qgis install location
QgsApplication.setPrefixPath(QgsApplication.prefixPath(), True)

# create a reference to the QgsApplication, setting the
# second argument to False disables the GUI
qgs = QgsApplication([], False)

# load providers
qgs.initQgis()

qgs.exitQgis()

In order to reproduce the error, I made a simple Github repository (https://github.com/mbrasebin/test_pyqgis_crash) with the code and a Travis configuration file (to set the environnement).
The log error can be viewed on travis pages.

If you have an idea to correct the error (or at least to catch it) do not hesitate to answer and/or to fork the repository.

Mickaël

------------------- EDIT

By launching the script with gdb, the error seems to be related to a QT class used through qgis. If you have an idea to avoid or catch this error ...

Thread 1 "python3" received signal SIGSEGV, Segmentation fault.
0x00007fffeaf23261 in QSqlDatabase::close() ()
   from /usr/lib/x86_64-linux-gnu/libQt5Sql.so.5


texts: Creating vector layer in memory following the QGIS Python Programming Cookbook

texts:
I got the QGIS Python Programming Cookbook Second Edition, in hopes of improving my abilities around QGIS as well as making me confident in using Python. Keep in mind that my programming abilities are virtually nil, but I figured if I followed the directions in the book I could at least begin to grasp the basics.
Well, the first real exercise has me landing straight on my face! In it, you are supposed to create a vector layer in memory to display a point. The following is the code provided:
layer =  QgsVectorLayer('Point?crs=epsg:4326', 'MyPoint' , "memory")
pr = layer.dataProvider()
pt = QgsFeature()
point1 = QgsPoint(20,20)
pt.setGeometry(QgsGeometry.fromPoint(point1))
pr.addFeatures([pt])
layer.updateExtents()
QgsMapLayerRegistry.instance().addMapLayers([layer])

So when I enter it one line at a time, when I get to the pt.setGeometry(QgsGeometry.fromPoint(point1)) line I get the following error:

Traceback (most recent call last):
File "C:\OSGEO4~1\apps\Python37\lib\code.py", line 90, in runcode
exec(code, self.locals)
File "", line 1, in
AttributeError: type object 'QgsGeometry' has no attribute 'fromPoint'

If I continue on, after inputting the last line I get the following error:

Traceback (most recent call last):
File "C:\OSGEO4~1\apps\Python37\lib\code.py", line 90, in runcode
exec(code, self.locals)
File "", line 1, in
NameError: name 'QgsMapLayerRegistry' is not defined

If I run it altogether through the script editor, the first error appears. I also considered that they meant fromPointXY, but I get the following error:

Traceback (most recent call last):
File "C:\OSGEO4~1\apps\Python37\lib\code.py", line 90, in runcode
exec(code, self.locals)
File "", line 1, in
TypeError: QgsGeometry.fromPointXY(): argument 1 has unexpected type 'QgsPoint'

I am out of ideas here. Since this is a book to learn PyQGIS, I that I should not be in this predicament in the first place! Either there is a major error on the FIRST exercise, or I am doing something incredibly stupid. The only other thing I can think of is that I am not using QGIS 2.18 as the book was written for, but 3.0.3 and 3.4.3 (I also tried upgrading to see if that would help, it didn't). I really don't want to downgrade just to learn PyQGIS, but if you think it would work I may just do so.
Here is a link to the page in question via Google Books

I will be downgrading to 2.18 in order to complete the book. I hope that learning on an older version will still be worthwhile.


texts: Updating mapnik stylesheet to render WGS84?

texts:
We have a production tile server loaded using osm2pgsql serving osm-bright, osm-carto and osm-mapquest style tiles and would like to be able to server those in WGS84 projection, as well.

I've tried in vain to update the XML files properly but Mapnik doesn't seem to recognize any of the Map.SRS definitions.

Most recently, per instructions at the Mapnik XMLConfigReference page I've tried using a map srs of +proj=latlong+datum=WGS84 which resulted in this error:

  An error occurred while loading the map layer 'osmcartoproj': failed
  to initialize projection with: '+proj=latlong+datum=WGS84 +no_defs' in
  Map at line 4 of '/usr/local/share/maps/style/osmreproj/osmreproj.xml'

Can anyone shed any light on the correct SRS for generating tiles in WGS84?


texts: JTS "points must form a closed linestring"

texts:
I'm trying to create a polygon and calculate it's area using

    final GeometryFactory gf = new GeometryFactory();

    final ArrayList<Coordinate> points = new ArrayList<Coordinate>();
    points.add(new Coordinate(0,0));
    points.add(new Coordinate(0,1));
    points.add(new Coordinate(1,0));
    points.add(new Coordinate(1,1));
    points.add(new Coordinate(0,0));

    final Polygon polygon = gf.createPolygon(new LinearRing(new CoordinateArraySequence(points
        .toArray(new Coordinate[points.size()])), gf), null);

    System.out.println(polygon.getArea());

And I get the error:

 points must form a closed linestring

When I tried to use the points :

points.add(new Coordinate(0,0));
points.add(new Coordinate(0,1));
points.add(new Coordinate(1,0));
points.add(new Coordinate(1,1));
points.add(new Coordinate(0,0));

I've received area = 0 .

I'd like to know what I've got to do in order to calculate the area correctly without getting any errors.


texts: ArcGIS Desktop tool that does for records in a table what "Dissolve" does for features?

texts:
I have a data table that I need to trim down, merging records which have the same value in an ID field and performing a couple math operations on various other fields (summing, averaging). The Dissolve tool in ArcGIS does this for feature classes, but it doesn't work on tables. I realize there's probably a tool in ArcGIS that does this (it's basically a Dissolve with no geometry) but I haven't found it and am in a bit of a hurry.

Is there one? I have an ArcInfo license available if necessary.


texts:
I am using QGIS 3.4 (Madeira).When attemptng to export my layout as jpeg or pdf, I get the message that the layout has been successfully exported, along with the path, but nothing shows up in the destination specified.
I have tried adjusting the export resolution, the destination and restarted the program. (After each attempt the asterisk appears after the project title which shows that the project hasn't been saved.)

Any ideas?


texts:
I'm trying to build a graph or a map which would show the below example route as a straight line and so would require the geometries of the countries to be reprojected about each leg of the route. It would be equivalent to making the below route into the X axis of a graph where you could easily see and measure the proximity of a coastline.
I have all the data in a PostGIS database and I'm mainly working in Python.
What I would like to know is how to reproject geometries so the they are rotated about the geometry of a linestring. I have seen a lot of promising information of Proj but no way to make custom projections on the fly which I believe is what I need here.
How can I do this?


texts: Converting multipoint layer to point using PostGIS?

texts:
I have a Spatialite layer which is multipoint (following a QGIS Save As from shapefile). I need it to be PostGIS and POINT. I have tried ogr2ogr -nlt POINT and got 'new row for relation "xxx" violates check constraint "enforce_geotype_geometry"' error on insert.

I've also tried converting within PostGIS - i.e. bringing the layer in as multipoint and converting to point - no joy.

I'm using PostGIS 1.5.3 PostGreSQL 9.1.8. (considering upping to PostGIS2 after reading this: How to change the geometry type from Point to Multipoint within an existing table in PostGIS?) but would prefer to stay at 1.5 for now if there's a way round this. Happy to do this in either spatialite or PostGIS.


texts: Custom function works in "Zoom to Features" but not "Select Features"

texts:
I'm trying to update the Using Custom Functions tutorial to QGis 3 and Python 3 and I've run into an interesting issue. With the following code I can create an expression GetUtmZone()='33N' and it works in the Output preview, it works to zoom to the features but if I try to select features it says no features are found!

Update

There seems to be no issue with the function code as it works fine as a label and a map tip.

import math
from qgis.core import *
from qgis.gui import *

@qgsfunction(args='auto', group='Custom', usesGeometry=True)
def GetUtmZone(feature, parent):
    centroid = feature.geometry()
    longitude = centroid.asPoint().x()
    latitude = centroid.asPoint().y()
    zone_number = math.floor(((longitude + 180) / 6) % 60) + 1

    if latitude >= 0:
        zone_letter = 'N'
    else:
        zone_letter = 'S'

    return '%d%s' % (int(zone_number), zone_letter)

The zoom to's are clearly working but each select returns the warning: No matching features found.

versions etc:

  QGIS version
  3.4.4-Madeira

  QGIS code revision
  f6ddc62

  Compiled against Qt
  5.9.5

  Running against Qt
  5.9.5

  Compiled against GDAL/OGR
  2.2.3

  Running against GDAL/OGR
  2.2.3

  Compiled against GEOS
  3.6.2-CAPI-1.10.2

  Running against GEOS
  3.6.2-CAPI-1.10.2 4d2925d6

  PostgreSQL Client Version
  10.6 (Ubuntu 10.6-0ubuntu0.18.04.1)

  SpatiaLite Version
  4.3.0a

  QWT Version
  6.1.3

  QScintilla2 Version
  2.10.2
  PROJ.4 Version
  493

  Linux astun-desktop 4.15.0-45-generic #48-Ubuntu SMP Tue Jan 29 16:28:13 UTC 2019 x86_64 x86_64 x86_64 GNU/Linux


texts: Processing GPS data before map matching using R?

texts:
I have GPS data collected with a smartphone app for all travel modes (e.g., car, bike, walking, transit). Trips were detected by the app (so I have trip_id in the data). In theory, each point is collected every 5 seconds; it varies in reality (5-15 seconds). The accuracy is pretty low (from 4 to around 70, some points are completely off with accuracy of ~150). This accuracy was measured in meters, I think 4 means that the point has 68% chance to fall inside a buffer of 4m.

My goals are to clean this dataset and use it to determine speed, travel distance, and evaluate the routes. With that, I'll snap it into the road network (with a map matching tool such as nearest_osm in R or similar).

What accuracy threshold should I use to filter out "bad" points before I can snap these points to the network?


texts: Creating full-resolution footprints in ArcGIS mosaic dataset?

texts:
In ArcGIS Desktop 10.3, I want to create full resolution footprints in a mosaic dataset to avoid having black areas in the mosaic. When I create footprints ("Build Footprints" tool) using the default 2000x2000 sampling (the "Request Size" parameter), the mosaic will have black areas along the edges of images:

These pixels are RGB=0,0,0. To avoid the black, I want the footprints to have pixel-level precision, and I don't mind waiting for the tool to calculate them. So I set the "Request Size" parameter to -1, which "will compute the footprint at the original resolution", according to the documentation. The result is indeed a new set of footprints with pixel-level precision. But unfortunately, the footprints don't actually follow the edge of the image data. The result is pretty perplexing:

Some 0,0,0 cells are included in the footprint, and some aren't. And it seems almost random where the line is drawn. These jagged boundaries are all over the mosaic. This could very well be a bug, but it's also entirely possible that I'm doing something wrong. What might I be missing?

In case it helps, here is the entire command, with all parameters (most were defaults):

arcpy.BuildFootprints_management(in_mosaic_dataset="mos", where_clause="", reset_footprint="RADIOMETRY", min_data_value="1", max_data_value="255", approx_num_vertices="-1", shrink_distance="0", maintain_edges="NO_MAINTAIN_EDGES", skip_derived_images="SKIP_DERIVED_IMAGES", update_boundary="UPDATE_BOUNDARY", request_size="-1", min_region_size="100", simplification_method="NONE", edge_tolerance="", max_sliver_size="20", min_thinness_ratio="0.05")


texts: Export current user location from locate control using button in Leaflet

texts:
I am preparing a mobile application, to determine the location of natural and cultural goods along a route that I have previously uploaded on Leaflet maps. I use the "leaflet locate control" plugin to establish my position and follow the route.

To capture the points of interest, and save them in a file along the route, I tried to use the code of TomazicM (Export user location in Leaflet) to export lat and lng, but "mousedown" doesn't work so well on mobiles. It works perfect in the computer and sometimes in the mobil.

What I've tried to do is change the "mousedown" option to a button that performs the same function, but I haven't found the solution.

I want to use the button to save the coordinates of points of interest during a route. Just one point (my location) each time I press the button.
How can I achieve this?

This is what I have tried so far but does not work:

lc = L.control.locate({
    locateOptions: {
              enableHighAccuracy: true,
    showCompass: true,
    maxZoom: 15,
    strings: {
            title: "position ??"
    }}
}).addTo(map);
var cnt;
L.easyButton('far fa-save', function(btn, map) {
  map.on('locationfound', function (evt) {
  setTimeout(function() {
    cnt = 0;
      cnt++;
      if (cnt > 1) {
        cnt = 0;
        return;
      }

      saveAsLocalFile(JSON.stringify(evt.latlng), 'location.txt', 'text/plain;charset=utf-8');
  }, 100);
});
}).addTo (map);

function saveAsLocalFile(data, filename, type) {
  var file = new Blob([data], {type: type});
  if (window.navigator.msSaveOrOpenBlob) // IE10+
    window.navigator.msSaveOrOpenBlob(file, filename);
  else { // Others
    var a = document.createElement("a");
    var url = URL.createObjectURL(file);
    a.href = url;
    a.download = filename;
    document.body.appendChild(a);
    a.click();
    setTimeout(function() {
      document.body.removeChild(a);
      window.URL.revokeObjectURL(url);
    }, 0);
  }
}


texts:
I am trying to use r.reclass function on a Land Cover raster in QGIS 3.4.
The operation is quite basic stuff and the settings are attached below.

However, I get an error that prevents the creation of my output raster.
The error should be "FileNotFoundError: [Errno 2]", as printed below.

I am running QGIS 3.4 on a MacBook. I do not get any error when performing the same operation on QGIS 3.2 on Windows.

I attach the screenshot of the Log Messages pane for more information.


texts: Creating spatial data on-the-fly for timeline application?

texts:
Facts:

Data are stored in a MySql db with no spatial data. instead there is a uri of gazetteer service which will return spatial data in json format.

Project:

build a timeline webgis application (openlayers ??) with these data.

My problem is how to proceed ..

Is there a framework which can produce on the fly GeoJSON data out of the box ?


texts: Creating a raster with point count as cell value?

texts:
I have a PostGIS table with Point geometries in it. I want to create a raster object where each cell value is the number of points within a given raster cell. Cells with zero points in it can be 0, Null or anything else that lets me differentiate it from a plausible count output.

What I did so far is to create a new raster (adapted from this tutorial):

CREATE TABLE public.sampleraster(
  rid SERIAL primary key, rast raster
);

INSERT INTO public.sampleraster(rast)
SELECT ST_AddBand(
  -- Make empty raster
  ST_MakeEmptyRaster(
      5000, -- Raster width x (in pixels)
      3000, -- Raster width y (in pixels)
      -75, -- Upper left X coordinate
      45, -- Upper left Y coordinate
      0.008333333333333, -- X Cell Resolution (in degrees)
      0.008333333333222, -- Y Cell Resolution (in degrees)
      0, -- X skew
      0, -- Y skew
      4326 -- SRID (WGS 84)
    ),
    ARRAY [
      ROW(
        1,  -- Band index: sets this as the first band
        '32BUI'::text,  -- Pixel Type (string rep of ST_BandPixelType types)
        NULL,  -- Initialized pixel value
        NULL   -- Nodata Value
      )
    ]::addbandarg[]
);

The Problem is now to fill the raster with count values. I'm counting the Points per cell like this

SELECT
round (st_x (geom) / 0.008333333333333) * 0.008333333333333 lon,
round (st_y (geom) / 0.008333333333222) * 0.008333333333222 lat,
count(*) ct
FROM sampledata
WHERE
st_x (geom) BETWEEN -75 AND -33.33333
AND
st_y (geom) BETWEEN 10 AND 35

GROUP BY round (st_x (geom) / 0.008333333333333), round (st_y (geom) / 0.008333333333222)
ORDER BY round (st_x (geom) / 0.008333333333333), round (st_y (geom) / 0.008333333333222)

Now I'm struggling to get the values into the raster. I think I need to get the values in a 2D array (going from the tutorial), but I don't know how to go about this. Any ideas here?

If somebody is interested in recreating the issue, you can use this query to generate a bunch of random points in the study area:

CREATE TABLE sampledata AS
WITH x AS
    (select generate_series (1,1000),
     st_setSRID(st_makepoint(
        (random()*50)::float - 120,
        (random()*25)::float + 30),
     4326) geom,
     random()*100::int num_a)
SELECT geom FROM x;


texts: Reduce extent of raster processing in ArcPy or ArcGIS Desktop?

texts:
I need to reduce the processing extent of my raster by some sort of automatic procedure in python. Currently the borders of the dataset I've been using has anomalies (due to the classification of the landsat scanning lines).

I would like to reduce the extent of the raster processing so that when I retrieve the information above by a SetNull function the borders are not processed. Something like

Now this can easily be done with an extract by mask function (i.e. figure2) although that would require manual selection of the masking region. Is there a menthod I can utilize to calculate the extent region (ignoring nodata values) minus the border width in the x and y direction of each row/column and then utilize the arcpy.env.extent in Python?


texts: Attribute Based Hierarchical Clustering Method took so long to process

texts: I am currently using the Attribute-based clustering plugin in QGIS. The K-Mean clustering method works and only took just about a minute to complete clustering. But when I ran the Hierarchical clustering, QGIS became "Not Responding" and even if I waited for almost 5 hours, the process still wasn't complete. Any advice on what to do. I am currently using the monthly precipitation of a region to determine which area has the same Precipitation pattern in a year.


texts: Inherited values between polygons in QGIS

texts:
In this project, I have three polygons, see the image below

'Polygons_From', includes input attributes, i.g. "Value_1", "Value_2", "Value_3", and $area for each quadrate (on the image there are two types of them).
'Polygons_To' polygons for which new values "Value_1_new", "Value_2_new", "Value_3_new" need to be defined
'Border' defines an area where the calculation should be done.

My task is to transfer values from 'Polygons_From' into 'Polygons_To' based on a certain aspect/ratio. Of course, I can simply achieve it with the 'Union' function as was suggested in this thread.
It uses $area-factor, but perhaps there are better approaches that will lead me to a much more efficient output?
So, I am asking for a recommendation regarding the technique/workflow that can be applied in QGIS 2.18. Shall I consider only $area or there are some other criteria?
Moreover, I looking forward to finding some papers that may improve my general understanding of the methodology.

References:

Calculating the area of one layer inside area of another layer QGIS
Calculating polygon area inside another polygon in QGIS
Combining data from overlapping polygons and assign
How to combine data from overlapping polygons?
Calculating area of polygon inside region?
Calculate proportional area of polygon within another layer's polygon


texts: Getting rid of cells of 0 values in raster using QGIS?

texts:
I have a DEM like below image. The black area shows 0-values.

I want to remove cells of 0-values in the raster.

In raster calculator, I wrote con("demclip_re@1">0,"demclip_re@1") but it gives "Incorrect expression".

What is wrong with my expression?


texts: PostGIS: Is it interesting to create an index on geom column and an alpha numeric field?

texts:
I would like to know if it's interesting to create an index on an other column than the geom one.

In some cases, I would like to do spatial analysis. If not, I would like to analyse only attributes.

So, should I make a second index on an alpha/numeric field?

EDIT: The practical case is few inserts (adn feature by feature) and a lot of selects on spatial or attributes (one main field as an ID (int) ) or both


texts: Can't display Mapbox tileset with layer control

texts:
I think the problem is in the mbURL. Simplified script:

    osmUrl    = 'https://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png',
    esriUrl   = 'https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServe\
r/tile/{z}/{y}/{x}',
    mbUrl = 'https://api.tiles.mapbox.com/v4/mtnbiker.d7jfhf8u.json?access_token=pk.eyJ1IjoibXRuYmlrZXIiLCJhIjoiNmI5ZmZjMzAyNzJhY2Q0N2ZlN2E1ZTdkZjBiM2I1MTUifQ.6R3ptz9ejWpxcdZetLLRqg';

var mbMap   = L.tileLayer(mbUrl),
    osmMap  = L.tileLayer(osmUrl),
    esriMap = L.tileLayer(esriUrl);

var map = L.map('map', {layers: [osmMap]})
           .setView([<%= @location.latitude %>, <%= @location.longitude %>], 17);

var baseLayers = {
  "OSM Street" : osmMap,
  "1908 LA [?]": mbMap,
  "Satellite"  : esriMap
};

L.control.layers(baseLayers).addTo(map);

You can see this in (in)action at alpha site and then click on stack on the right and select 1908 LA. No errors, but no display. I can display the tiles—click on Map on the nav bar, it's the same map I want to display. Something wrong in the syntax of mbURL, but I've looked searched and tried many alternatives. (The lat and long are Rails calls to a PostgreSQL database and they work, because they also display the marker. The ? mark in 1908 LA is to suggest the user that there may be a problem.)

Code is at https://git.heroku.com/secure-shore-68966.git.
Tileset is at https://www.mapbox.com/studio/tilesets/mtnbiker.d7jfhf8u/.


texts: Existing solutions to automate generalisation/simplification of geographical data?

texts:
I am looking for existing solutions (softwares, libraries, methods, etc.) to automate the generalisation/simplification of geographical data as presented there, and also there.

Do you know and use some tools to simplify your geographical data?

I already had a look at these ones:

maptimize to simplify the representation of KML marker layers,
Modelbuilder by ESRI
Radius clarity by laserscan
mapshaper to compute online filtering algorithms on geometries,
Axpand by axes-systems

The ICA commission on generalisation and multiple representation website is a huge information source on this topic. Many papers presents advanced methods and use cases related to this topic.

The very famous douglass & peucker filtering algorithm is often used for geometrical simplification and is present in almost all GIS softwares. However, its effect is to make geometries very bony/angular without really simplifying the geometry shape. For example, it is not very adapted for lines representing roads or rivers, that have to be smooth and not angular.


texts: qgis's georeferencer renders my white map with a translucent blue tint, can I prevent this?

texts:
I'm running QGIS 2.0.1-Dufour on MacOSX and when loading a PNG or TIFF into qgis as a raster its always getting rendered with a blue tint.

Here's what it looks like before being loaded into QGis's Geoferencer.

and after:

Is there any way I could save this file in photoshop or gimp that would allow it to maintain its original colorization?

GDALINFO (png)

    davodavo$ gdalinfo color\ key.png
    Warning 1: libpng: iCCP: Not recognizing known sRGB profile that has been edited
    Driver: PNG/Portable Network Graphics
    Files: color key.png
    Size is 308, 588
    Coordinate System is `'
    Metadata:
      XML_com.adobe.xmp=<x:xmpmeta xmlns:x="adobe:ns:meta/" x:xmptk="XMP Core 5.4.0">
       <rdf:RDF xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#">
          <rdf:Description rdf:about=""
                xmlns:exif="http://ns.adobe.com/exif/1.0/">
             <exif:PixelXDimension>308</exif:PixelXDimension>
             <exif:PixelYDimension>588</exif:PixelYDimension>
          </rdf:Description>
       </rdf:RDF>
    </x:xmpmeta>

    Image Structure Metadata:
      INTERLEAVE=PIXEL
    Corner Coordinates:
    Upper Left  (    0.0,    0.0)
    Lower Left  (    0.0,  588.0)
    Upper Right (  308.0,    0.0)
    Lower Right (  308.0,  588.0)
    Center      (  154.0,  294.0)
    Band 1 Block=308x1 Type=Byte, ColorInterp=Red
      Mask Flags: PER_DATASET ALPHA
    Band 2 Block=308x1 Type=Byte, ColorInterp=Green
      Mask Flags: PER_DATASET ALPHA
    Band 3 Block=308x1 Type=Byte, ColorInterp=Blue
      Mask Flags: PER_DATASET ALPHA
    Band 4 Block=308x1 Type=Byte, ColorInterp=Alpha

gdalinfo (tiff) -- same blue tint

    Driver: GTiff/GeoTIFF
    Files: colorkey.tif
    Size is 308, 588
    Coordinate System is `'
    Metadata:
      TIFFTAG_DATETIME=2014:07:09 17:27:31
      TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
      TIFFTAG_SOFTWARE=Adobe Photoshop CS6 (Macintosh)
      TIFFTAG_XRESOLUTION=143.9926
      TIFFTAG_YRESOLUTION=143.9926
    Image Structure Metadata:
      INTERLEAVE=PIXEL
    Corner Coordinates:
    Upper Left  (    0.0,    0.0)
    Lower Left  (    0.0,  588.0)
    Upper Right (  308.0,    0.0)
    Lower Right (  308.0,  588.0)
    Center      (  154.0,  294.0)
    Band 1 Block=308x8 Type=Byte, ColorInterp=Red
    Band 2 Block=308x8 Type=Byte, ColorInterp=Green
    Band 3 Block=308x8 Type=Byte, ColorInterp=Blue


texts: How do I know what parameters I can call for WMS layers?

texts:
I'm trying to learn OpenLayers but the documentation is hard to understand.

For WMS layers, ex:

 var wms = new OpenLayers.Layer.WMS( "OpenLayers WMS",
            "http://vmap0.tiles.osgeo.org/wms/vmap0", {
                "transparent":"true",
                "layers":"basic",
                "format":"image/png" },
              { "reproject":"true" }
              );

How do I know what all the possible parameters are?

Is there a list somewhere of all the parameters?


texts: gdal_translate creates images that are mirrored

texts:
I'm having a problem with gdal_translate - when translating from an XYZ file to AAIGrid or geotiff the raster is mirrored and "moved" south on one axis.

The top part of the attached image is the original file, when displayed in QGis, the bottom file has been translated to AAIGrid format with gdal_translate using

gdal_translate -of AAIGrid -a_srs EPSG:3301 grid.xyz grid.grd

The same happens when I try to translate the XYZ file to Geotiff.

Edit: added file sample.
The file looks like this, the first 2 colums are coordinates in EPSG:3301 local coordinate system and the third column is height in metres

719935 6575005 30.709999
719945 6575005 31.08
719955 6575005 30.805
719965 6575005 30.772499
719975 6575005 30.2775
719985 6575005 30.1175
719995 6575005 30.012501
715005 6575015 28.525
715015 6575015 28.715
715025 6575015 28.834999
715035 6575015 28.6875
715045 6575015 28.452499
715055 6575015 29.147499

Output from gdalinfo

gdalinfo 65714_dem_10m.xyz.grid
Driver: AAIGrid/Arc/Info ASCII Grid
Files: 65714_dem_10m.xyz.grid
       65714_dem_10m.xyz.prj
Size is 500, 500
Coordinate System is:
PROJCS["Estonian_Coordinate_System_of_1997",
    GEOGCS["GCS_EST97",
        DATUM["Estonia_1997",
            SPHEROID["GRS_1980",6378137,298.257222101]],
        PRIMEM["Greenwich",0],
        UNIT["Degree",0.017453292519943295]],
    PROJECTION["Lambert_Conformal_Conic_2SP"],
    PARAMETER["standard_parallel_1",59.33333333333334],
    PARAMETER["standard_parallel_2",58],
    PARAMETER["latitude_of_origin",57.51755393055556],
    PARAMETER["central_meridian",24],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",6375000],
    UNIT["Meter",1]]
Origin = (715000.000000000000000,6575000.000000000000000)
Pixel Size = (10.000000000000000,-10.000000000000000)
Corner Coordinates:
Upper Left  (  715000.000, 6575000.000) ( 27d46'16.00"E, 59d15'32.02"N)
Lower Left  (  715000.000, 6570000.000) ( 27d45'58.29"E, 59d12'50.69"N)
Upper Right (  720000.000, 6575000.000) ( 27d51'31.03"E, 59d15'22.83"N)
Lower Right (  720000.000, 6570000.000) ( 27d51'12.91"E, 59d12'41.51"N)
Center      (  717500.000, 6572500.000) ( 27d48'44.56"E, 59d14' 6.79"N)
Band 1 Block=500x1 Type=Float32, ColorInterp=Undefined