texts: pgrouting pgr_dijkstra forcing route along a selected edge?

texts:
I have routing working where a user clicks 2 points on the map, I find the closest vertex to those 2 clicked locations and pass that to a pgr_dijkstra function.

What I would like is give the edge a user clicks a dynamically lower cost, which would force the route to go that way, even if there is another route that is the same cost, or maybe even less.

So given this image, a user clicks A & B, the route could go either way, as the length/cost is the same. But since the user clicked the "B" on the left edge, that should be given priority.

Is this possible?

This is my query now:

SELECT id, geom, trailid, cost, reverse_cost, source, target, ST_Reverse(geom) AS flip_geom
FROM pgr_dijkstra(
   'SELECT id, source::int, target::int, length::float AS cost, reverse_cost::float AS reverse_cost
   FROM edge_table
   WHERE ST_Contains(bbox, 4326), geom)',
source, target, false, true), edge_table
WHERE id2 = id ORDER BY seq;

This is a prime example, I could do either way, but I would prefer to route up the green climbing trail.
http://www.trailforks.com/region/vedder-mountain/planner/?waypoints=49.076580,-121.986912:49.070167,-122.000207


texts: MySQL primary key column not visible in QGIS 2.6

texts:
I've tried adding a vector from Mysql table (Mysql 5.1, MyISAM) to QGIS 2.6.
All columns are visible except Primary key column(s).

Could anyone confirm this behaviour? Could we consider it a bug?


texts: Suitability Analysis with Raster Calculator and Weighted Values

texts:
I am trying to perform a suitability analysis based on certain criteria which I have rasters of. I would like to use the Raster Calculator function in ArcMap to determine suitability scores in a given area, with some criteria weighted more than others.

For example, one criterion is income, another is tree canopy percentage, and a third is percentage of vacancies. I'd like to give a equal ranking/weight to tree canopy and vacancy and half weight to income.

How should I input this into the raster calculator?

I do understand that if I wanted to give all equal weight I would just input:
 "treeCanopy" + "vacancy" + "income"
and get my resultant scores in the output raster.


texts: Read Pixel Values from WMTS RasterLayer in PyQGIS

texts:
I'm using a WMTS service to serve image tiles, and I want to read the pixel values of these tiles over some extent.

I am able to do this with a QgsRasterLayer of a local file (.tif) without an issue, using the following:

import numpy as np
from qgis.core import *

layer = QgsProject.instance().mapLayersByName('MyLocalTiff')[0]
provider = layer.dataProvider()
extent = provider.extent()
height, width = layer.height(), layer.width()

img = np.zeros((height, width, 3), np.uint8)

for band in range(3):
    block = provider.block(band, extent, width, height)

    for y in range(0, height):
        for x in range(0, width):
            img[y, x, band] = block.value(y, x)

However, when I use the same technique with a WMTS RasterLayer, every pixel value is NaN. QGIS is able to display the layer just fine and all the image data is displayed on the canvas without issue.

How do I access these pixel values?


texts: where do I download POIs from Wikipedia? placeopedia used to have them but is now gone

texts:
I am trying to improve OpenStreetMap in Siem Reap, Cambodia.

A way I am doing this is to add temple locations. I hear that there is a dataset of temple locations available from Wikipedia (presumably amongst a broader dataset of POIs on Wikipedia).

How to extract these POIs from Wikipedia? Placeopedia used to have them, apparently, but now it's down.

One way I tried to do this was to download Wikipedia POIs in the Android app called OSMAND+. This, however, is in .OBF format.

Does anyone know how to convert from .OBF format to another one like .KML so I can view the file in JOSM and add missing temples?


texts: Extracting values from raster under line using GDAL?

texts:
How can I find the cells which lie on a straight line between two points?
I.e. how can I find the indices of the cells (in grey) which are intersected by the line defined by the two points (in blue):

I want to do this in order to read off values from the underlying DTM and thus extract the height profile along the line....

Using the geotransform of the raster, it is easy to convert the two points to cell offsets, so my first thought was to create many points along the line and convert each of these. However the problem here is creating points at the right distance to successfully capture each intersected cell (while also not creating a ridiculous number of points so that it would be incredibly inefficient)

Looking for ideas that are either programmatic (python/C++) or theory (e.g. no arcgis tools! gdal ok...)


texts: Might my file be too big to be written with 'writeRaster' in R?

texts:
I am trying to write a GeoTiff of a rasterstack that I created after some calculations. I am not quite sure how big the output-file might be. But I assume it has more than 4 GB. I can plot the stacked file but everytime I try to write the file, R ist busy for several hours and I cannot spot anything in my output folder. I tried the code with some smaller images and it seems to work fine.

Is there a size-limit for 'writeRaster' in R and, if so, is there a workaround?

I read there is BigTiff, but it seems not to be supported within the formats of writeRaster.

 rs <- writeRaster(rasterstack, filename=
                paste(getwd(),sprintf("/Raster/%s/rasterstack_%s.tif",
                                      year,year),sep = ''),format="GTiff",overwrite=TRUE)


texts: Creating Polygons from black/white map or points using QGIS?

texts:
I just started working with qgis and I have two maps of a town. On the one map buildings are outlined by Points, like this.

And the other map is just black with White Buildings. like this.

My task is now to create a layer in which the buildings are polygons. I looked up some questions here and found some interesting stuff, but it never really worked for me. I tried some tools in qgis, like lines to Polygons, but it never worked.

Can you think of a way to auto-create a layer like this or should I just do it by hand?

Of course I'm asking this because it's a very big amount of data and I would rather not do it by hand.

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

This is what the identifier gives me. "Kanal" is german for "channel" or "band" I think.


texts: Distance between centriod of the polygon and each and every point of the polygon

texts:
I'm trying to calculate the distance between centriod of the polygon and each and every point of the polygon using GeoPandas Python

I need to find the distance between centroid(18.9568383649  72.831281) and geometry POLYGON ((72.83164071429999 18.9568102857,
                                    72.8314213636 18.9569958182, 72.831512
                                   18.9569766667,72.83164071429999
                                   18.9568102857))

How I can find the distance? I'm not tied to using GeoPandas.


texts: Has Export to AI (Adobe Illustrator) changed at ArcGIS for Desktop 10.2?

texts:
In the past when I would export vector layers from ArcGIS each layer would be nested in a group with a clipping path over the top.

Now when I Export to AI from 10.2 that layer's features are 5-6 layers deep.

Did this software behaviour change at 10.2?


texts: Z-Index in GeoJSON features using leaflet

texts:
I have a watershed polygon with sub-watersheds contained within the larger watershed. I would like the sub-watersheds drawn on top.

Can i control the drawing order of features within a geojson object or do i have to add the features separately and then control the layer drawing order?

I'm using leaflet 0.7.2


texts: R: Error in checking point inside polygon

texts:
I am trying to check if location point is inside the polygon and it works for certain location point. But for certain location point output says point is outside polygon even if it is inside.

Rmap <- readShapePoly('MAP.shp')#rome shape file
mappoints <- fortify(Rmap)

data_raw<-data.frame(Longitude=12.477,Latitude=41.91020)
ind<-point.in.polygon( data_raw$Longitude, data_raw$Latitude,    mappoints$long, mappoints$lat)

ggplot(mappoints) + aes(long, lat) + geom_polygon(fill="grey80",aes   (group=group))  + scale_fill_manual(values = pretty_palette)+ geom_point(data=data_raw,aes(x=Longitude,y=Latitude))

Output:
ind is always 0 which means point is outside polygon. Original plot is


texts:
Insert Cursor fails to write results to table. The Search Cursor is returning the results, but errors when trying to insert into field.

# Insert Cursor for GDB table to GDB table
# Import modules, environment settings
import arcpy, traceback, string
arcpy.env.overwriteOutput = True

# Database connection variables on SPGISPROC1
OutFGDB = r"P:\batch\APO\District_overlay.gdb\pa_fp_dat1"
InFGDB = r"P:\batch\APO\District_overlay.gdb\pa_fp_rel1"

#Create cursor for input table items
InCursor = arcpy.SearchCursor(InFGDB)

# Create a variable that stores the value for column in a given row of the input table

for InRow in InCursor:

    parcel      = InRow.getValue('pafp_frq_PARCEL')
    fp_zone1    = InRow.getValue('pafp_frq_ZONE2')

    # Copy the stored value from the input table into the output table
    gdbCursor = arcpy.InsertCursor(OutFGDB)
    row = gdbCursor.newRow()

    if fp_zone1:
        fp_zone1 = fp_zone1.upper

    #Fill out fields with variables from input table
    row.pafp_frq_PARCEL = parcel
    row.pafp_frq_ZONE2  = fp_zone1

    #Insert new row into output file
    gdbCursor.insertRow(row)

#Delete cursor and row objects to remove locks on the data
del row
del cursor

Error Message:

Traceback (most recent call last):
  File "C:\Python27\ArcGISx6410.1\Lib\site-packages\pythonwin\pywin\framework\scriptutils.py", line 325, in RunScript
    exec codeObject in __main__.__dict__
  File "C:\gisdata\Python\Insert_table_Cursor.py", line 28, in <module>
    row.pafp_frq_PARCEL = parcel
  File "C:\Program Files (x86)\ArcGIS\Desktop10.1\arcpy\arcpy\arcobjects\_base.py", line 35, in __setattr__
    return setattr(self._arc_object, attr, ao)
RuntimeError: ERROR 999999: Error executing function.


texts: Does qgis support postgresql range types?

texts: I added a postgis layer with date range column in Qgis but that field is not listed in Qgis.


texts: How to perform Spatial Query to remove features that intersect/touch/overlap in QGIS?

texts:
I have two polygon layers and many of the features are identical (or should be identical). I want to isolate only the features that are in one layer but not the other.
I'm using QGIS 3 and I tried using the Vector>Geoprocessing>Difference tool but it left me with slivers of features because they weren't perfectly identical.

I used this layer with slivers by eliminating all features with an area ~0, but I'd still like to know if there's a Spatial Query I should be running.


texts: How to export .mdb file to shapefile in MapInfo Professional?

texts:
I have received some GIS files, in .mdb format.  Is it possible to project this .mdb file onto my map; say by converting it into a shapefile etc?

In MapInfo Professional, I could open this .mdb file but I don't know if I could translate it to another format.


texts: R Freezes when Plotting Spatial Data Frames

texts:
I have run into a problem where anytime I try and plot a Spatial DataFrame in R, the console freezes, regardless of the size of the file.  I imagined this must be due to an issue with the path to one of the packages or the GDAL framework, but I have uninstalled and re-installed all of them and am still having an issue.  Any ideas?

Some example code, although I doubt that's the issue, as it is happening with all SpatialDataFrames:

library(raster)
india <- getData("GADM", country="IND", level=1)
plot(india)

I don't get an error, it just runs indefinitely until I manually force quit R.


texts:
I am trying to create tiles given a collection of features.

UPDATE

Let's say I got two tiles per zoom level 1 and my features are about a small area near Berlin. I expect one of the two tiles empty and the other one filled with all the features with proper dimension (I expect to see something like a point because features are about a small area). What I get is an empty tile (correct) and a image 256 x 256 representing the features like it was a thumbnail, not a tile.

This is what I did so far to render features in a tile:

...
MapContent content = new MapContent();
GTRenderer renderer = new StreamingRenderer();
BufferedImage image = new BufferedImage(width, height, BufferedImage.TYPE_INT_ARGB);
Graphics2D graphics = image.createGraphics();
ReferencedEnvelope mapArea = featureCollection.getBounds();
Style shpStyle = SLD.createSimpleStyle(featureCollection.getSchema());
Layer layer = new FeatureLayer(featureCollection, shpStyle);
content.addLayer(layer);
renderer.setMapContent(content);
Rectangle tileArea = new Rectangle(0, 0, 256, 256);
renderer.paint(graphics, tileArea, mapArea);
File fileToSave = File.createTempFile("mytesttile.png");
ImageIO.write(image, "png", fileToSave);
...

but this code produces a small (256 x 256) image of the given collection, not a tile. I am in short of time and I can't find some good docs or a working example. Can you point me to the right direction? I am using geotools version 17.2.


texts: Storing information from a Combobox in layer in QGIS/PyQt

texts:
I am building a plugin for QGIS and want to add a string to a created polygon. In my dialog I was able to create a ComboBox filled with 10 different strings which can be added while creating the polygon. But when I identify this object this information doesn't seem to be stored. Is there a way to connect this information to the layer? Here is what I created until now:

Objects = ['Trees on grass', 'Trees on avenue', 'Grass field', 'Green roofs', 'Cool roofs', 'Facade greening', 'Solar Panels', 'Pond/river', 'Fountain', 'Surface albedo change']

#insert this string in Combobox
self.ObjectType = QComboBox()
object = self.ObjectType.insertItems(0, Objects)

# Add this object feature to layer
for layer in QgsProject.instance().mapLayers().values():
   if 'field=' + self.tr('Drawings') + ':string(255,0)' in layer.dataProvider().dataSourceUri()[-28:]:  # must have its first field named Drawings, string type
     self.layers.append(layer)
     self.layerBox.addItem(layer.object)

Does anyone have a clue to enhance the last part?


texts: ImportError: DLL load failed after QGIS Upgrade from 3.14.1 to 3.14.16

texts:
I have upgrade my QGIS from 3.14.1 to 3.14.16 and I get this error on startup.
Traceback (most recent call last):
  File "", line 1, in
  File "C:/OSGEO4~1/apps/qgis/./python\qgis\__init__.py", line 70, in
    from qgis.PyQt import QtCore
  File "C:/OSGEO4~1/apps/qgis/./python\qgis\PyQt\QtCore.py", line 24, in
    from PyQt5.QtCore import *
ImportError: DLL load failed: Kan opgegeven procedure niet vinden.

Python versie:
3.7.0 (v3.7.0:1bf9cc5093, Jun 27 2018, 04:59:51) [MSC v.1914 64 bit (AMD64)]

QGIS versie:
3.14.16-Pi 'Pi', a235a149f0

Python pad:
['C:/OSGEO4~1/apps/qgis/./python', 'C:/Users/joeri/AppData/Roaming/QGIS/QGIS3\\profiles\\default/python', 'C:/Users/joeri/AppData/Roaming/QGIS/QGIS3\\profiles\\default/python/plugins', 'C:/OSGEO4~1/apps/qgis/./python/plugins', 'C:\\OSGEO4~1\\apps\\Python37', 'C:\\OSGEO4~1\\apps\\Python37\\Scripts', 'C:\\OSGeo4W64\\bin\\python37.zip', 'C:\\OSGEO4~1\\apps\\Python37\\DLLs', 'C:\\OSGEO4~1\\apps\\Python37\\lib', 'C:\\OSGeo4W64\\bin', 'C:\\Users\\joeri\\AppData\\Roaming\\Python\\Python37\\site-packages', 'C:\\OSGEO4~1\\apps\\Python37\\lib\\site-packages', 'C:\\OSGEO4~1\\apps\\Python37\\lib\\site-packages\\win32', 'C:\\OSGEO4~1\\apps\\Python37\\lib\\site-packages\\win32\\lib', 'C:\\OSGEO4~1\\apps\\Python37\\lib\\site-packages\\Pythonwin']

My toolbox processing is gone, I even cannot install anymore plugins, please help I need to get a certain QGIS model back.


texts: Migrating ArcSDE data to Oracle Spatial?

texts:
We've got all the geospatial information of our company within an ArcSDE Database (over oracle  Server):

points
lines
polygons

We consume these data with ArcGIS Desktop and with ArcGIS Server.

We are planning to migrate all the information to Oracle Spatial.

Is using Spatial ETL a way to do that ?


texts:
As part of work in statistical mapping for school, I need to use QGIS to find the spatial extent of commercial nodes.  A node is defined for this exercise as a zone of 75m in diameter where the density of shops is greater than or equal to 5.

The problem is less one of finding the spatial extent of the nodes than to find and map the nodes themselves.  To find the nodes, I have 5 shapefiles, each corresponding to a type of commerce in Villeurbanne (the study area for this exercise).  I have merged these five shapefiles into a single one, in the hope of simplifying my life (but perhaps that's a mistake...). I therefore have available a vector file where 2018 points represent the various shops in Villeurbanne.

What I would like to do is to analyze a circle of 75m in diameter (or a radius of 37.5m) around each point representing a shop, in order to know how many shops are found in this circle.  I don't know whether it's possible with vector software... I know that a similar operation exists in raster software like ILWIS for analyzing data in continuous grids, because I have done that before, but I have to work in QGIS for this excercise.

Once the shops have been found within a 75m diameter zone, I will have to be able to create the regions of intersection, which would give me the nodes.

I'm a beginner with QGIS, an application that I have used only once or twice last year by following instructions, but this year no instructions are indicated.

Dans le cadre d'un travail à réaliser en carto statistique pour la fac, je dois déterminer l'étendue spatiale de nodules commerciaux sur Q-Gis. Un nodule est défini dans l'exercice comme une zone de 75m de diamètre où la densité de commerces est supérieure ou égale à 5.

Mon problème n'est pas tant de déterminer l'étendue spatiale des nodules que de déterminer et cartographier les nodules eux-mêmes. Pour déterminer les nodules, je dispose de 5 shapefile correspondant chacun à un type de commerce à Villeurbanne (la zone d'étude de l'exercice). J'ai rassemblé ces cinq shapefile en un seul, dans l'espoir de me simplifier la vie (mais c'est peut-être une erreur...). J'ai donc comme base de travail un fichier vectoriel où 2018 points représentent les différents commerces de Villeurbanne.

Ce que je souhaiterais faire, c'est analyser un cercle de 75m de diamètre (ou 37,5m de rayon) autour de chaque point représentant un commerce, pour savoir combien de commerces se trouvent dans ce cercle. Je ne sais pas si c'est possible sur un logiciel vectoriel... Je sais qu'une opération similaire existe sur les logiciels raster comme Ilwis en analysant l'information contenue dans les mailles contigues pour l'avoir déjà fait, mais je dois travailler sur Q-Gis pour cet exercice.

Une fois les commerces trouvés dans une zone de 75 m de diamètre, il faudrait que je puisse regrouper les zones d'intersection, ce qui me donnerait mes nodules.

Je suis novice sur Q-Gis, logiciel que je n'ai utilisé qu'une ou deux fois l'an dernier en suivant des fiches des consignes, mais cette année, aucune procédure n'est indiquée...

J'espère que ma question est claire

Merci d'avance à tous.


texts:
How does one go about creating a convex hull from a set of points in SAGA GIS? I read that this is possible in QGIS, but I prefer SAGA and would like to see if it is possible.

Furthermore, if it is possible to create said convex hull, is it possible to create multiple convex hulls around points with the same value of a classification variable with one command?

I have a few hundred thousand observations with fairly accurate latitude and longitude values that also have their respective FSA from their postal code. I would like to create a polygon map of said city with each FSA as a separate polygon so I can do some aggregate analysis on neighborhoods created from manually unioning these polygons.

I know Statistics Canada provides an FSA boundary file (I used it fairly extensively when I worked for them and do have access to it for free through my University), but I'd prefer to learn how to do this on my own.


texts: JoinField_management Appears to Append Field Values

texts:
In the snippet of code below, I'm attempting to create variables, then pass them to the arcpy.JoinField_managment tool:

#Set join parameter variables
inFeatures = arcpy.GetParameterAsText(0)
inField = arcpy.GetParameterAsText(1)
joinTable = arcpy.GetParameterAsText(2)
joinField = arcpy.GetParameterAsText(3)
arcpy.MakeFeatureLayer_management(inFeatures,"inFeatures_lyr") #make feature layer from feature class
#Join road feature class to buffer distance table
arcpy.JoinField_management("inFeatures_lyr", inField, joinTable, joinField)

The result appears to actually append the field values from joinTable to the attribute table in inFeatures_lyr.  I've tried replacing the variables with hard-coded values and get the same result.  I can't remove the join and and can only get the fields out by deleting them as I would any other field in the feature class.  How do I structure this to create the join correctly, do my geoprocessing, then remove the join programmatically?


texts: Where is the ArcGIS Workflow Manager REST API documentation?

texts:
I hope it is online and I'm just missing with my searches.

I'm looking for this API documentation for 10.0 (but I will +1 answers relating to the 10.1 version):


texts: cartodb.js map doesn't display custom legend or attribution

texts:
I have developed a map at https://afterfibre.nsrc.org/af_dev.html which should display a custom attribution in the lower right corner:

 L.tileLayer('https://dnv9my2eseobd.cloudfront.net/v3/cartodb.map-4xtxp73f/{z}/{x}/{y}.png', {
    attribution: 'Undersea cable data courtesy Greg Malknecht
Mapbox Terms &amp; Feedback'
    }).addTo(map_object);

and should have a custom legend:

      var aftLegend = new cdb.geo.ui.Legend.Custom({
        title: "Fibre Status",
        data: [
          { name: "Live",               value: "#58A062" },
          { name: "Under Construction", value: "#F07971" }
        ]
      });
      $('#map_object').append(aftLegend.render().el);

but no legend appears and the attribution remains the default.  The map  seems to work well in every other respect.  No javascript errors detected.  The full code is at

https://github.com/stevesong/nsrc-afterfibre/blob/master/af_dev.html

I have been through it several times but am at a loss to work out what I am doing wrong.


texts: Georeferencing in ArcGIS - Zero-Order (Shift)

texts:
Zero-Order (shift) transformation is meant for rasters that only require shift, which as I understand it means moving or translating the image. Two questions:

Does it mean that two or more control points always yield residuals? How does this transformation handle two or more control points?
When I do use more than two points it results in an extremely enlarged raster (each pixel kilometers across) which looks like a bug but I assume could be the result of how the transformation works. Any way to explain this?


texts: Problem with integration of GRASS-Gis in creating a process for PyWPS

texts:
Hey everyone:) It's possible that I have some kind of miss-conception in my head and just wanted to ask it here. I have a space-time-datacube in GRASS-Gis and now I'd like to create a process in PyWPS that lets me calculate some statistics for an index (e.g. the ndvi).
I do have a script and I can create the process and see it on localhost. However, when I send a request to the server I always get:
PS&identifier=sen2_index&version=1.0.0 HTTP/1.1" 200 -
127.0.0.1 - - [22/Jun/2020 16:58:14] "GET /favicon.ico HTTP/1.1" 404 -
Traceback (most recent call last):
  File "/usr/local/lib/python3.6/dist-packages/pywps/app/Process.py", line 236, in _run_process
    self._set_grass(wps_request)
  File "/usr/local/lib/python3.6/dist-packages/pywps/app/Process.py", line 362, in _set_grass
    from grass.script import core as grass
ModuleNotFoundError: No module named 'grass'
127.0.0.1 - - [22/Jun/2020 16:59:10] "GET /wps?request=Execute&service=WPS&identifier=sen2_index&version=1.0.0&datainputs=start=2019-07-01;end=2019-12-01 HTTP/1.1" 200 -

So it apparently doesn't find Grass. I'd like to give some more information about my installation etc., but I'm really sure how to access these informations. I do have miniconda installed and I always do a conda deactivate prior to start GRASS. At least to start the gui this is for some reason necessary (I think for matching the python-version).
If someone has any comment on what I could do or how to proceed, would be awesome:)
Update
I configured my pywps.cfg and modified the line gisbase=/usr/bin/grass78
to
gisbase=/usr/lib/grass78
What I get now in my console is this here:
 "GET /wps?request=DescribeProcess&service=WPS&identifier=sen2_index&version=1.0.0 HTTP/1.1" 200 -
Traceback (most recent call last):
  File "/usr/lib/grass78/etc/python/grass/lib/ctypes_loader.py", line 76, in load
    return ctypes.cdll.LoadLibrary(path)
  File "/usr/lib/python3.6/ctypes/__init__.py", line 426, in LoadLibrary
    return self._dlltype(name)
  File "/usr/lib/python3.6/ctypes/__init__.py", line 348, in __init__
    self._handle = _dlopen(self._name, mode)
OSError: libgrass_datetime.7.8.so: cannot open shared object file: No such file or directory

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/usr/local/lib/python3.6/dist-packages/pywps/app/Process.py", line 236, in _run_process
    self._set_grass(wps_request)
  File "/usr/local/lib/python3.6/dist-packages/pywps/app/Process.py", line 414, in _set_grass
    from grass.pygrass.gis import make_mapset
  File "/usr/lib/grass78/etc/python/grass/pygrass/gis/__init__.py", line 13, in <module>
    import grass.lib.gis as libgis
  File "/usr/lib/grass78/etc/python/grass/lib/gis.py", line 23, in <module>
    _libs["grass_gis.7.8"] = load_library("grass_gis.7.8")
  File "/usr/lib/grass78/etc/python/grass/lib/ctypes_loader.py", line 62, in load_library
    return self.load(path)
  File "/usr/lib/grass78/etc/python/grass/lib/ctypes_loader.py", line 78, in load
    raise ImportError(e)
ImportError: libgrass_datetime.7.8.so: cannot open shared object file: No such file or directory
127.0.0.1 - - [23/Jun/2020 08:24:54] "GET /wps?request=Execute&service=WPS&identifier=sen2_index&version=1.0.0&datainputs=start=2019-07-01;end=2019-12-01 HTTP/1.1" 200 -


texts: Determining search radius used in Kernel Density analysis?

texts:
When conducting kernel density analysis in ArcGIS (10.2) if you leave the search radius field blank I realise it uses a default based on an algorithm.

Is there any easy/quick way to find out what search radius was applied after a new surface is created? In the geoprocessing results it's left blank.


texts: Why is Majority so much slower than mean in Focal statistics?

texts:
I am trying to use this solution to clean up small holes in my rasters. It uses focal statistics on the no data cells to replace them with a local average. The raster have a small number of unique values so using MEAN in focal statistics increases the number of values and gives dodgy results (I want a hole in a section where value = 100 to be 100 and not 98.9). The shouldn't be a problem because I can use MAJORITY in focal statistics.

What is stange though is that this python code takes 9 minutes to run:

filled = sa.Con(sa.IsNull(in_raster),sa.FocalStatistics(in_raster,
                     sa.NbrCircle(50, "CELL"),'MEAN'), in_raster)

still a little slow but manageable.

Then this code took 4 hours and 26 minutes to run (on the exact same in_raster):

filled = sa.Con(sa.IsNull(in_raster),sa.FocalStatistics(in_raster,
                    sa.NbrCircle(50, "CELL"),'MAJORITY'), in_raster

The only difference is mean vs Majority. Could anyone explain to me why there is such a difference in performance and if there is anything I can do about it?

I have a slightly laborious work-around of reclassifying the output to return to my small set of unique values, but I'd like to understand what seems like a simple solution is unworkably slow


texts: ol v.2.13 - Selecting features with draw or zoombox

texts: I'm attempting to add a feature to a legacy mapping project. I need the user to be able to draw a rectangle and have the features/markers inside of it selected and pushed into the selected array. Is there an easy way to go about doing this? I would be fine simply adding a select feature to the already existing zoomBox function to make it easier. I've had difficulty implementing the examples I've seen and am new to OpenLayers. I'm working in ol v. 2.13.


texts: Using layer name with accent on letter in ArcPy?

texts:
I have a lot of MXDs (more than 200) and all of them use a polygon shapefile "study_area_v00.shp".
In the data frame, its name is just "Study Area".

Every time we change the boundary of the polygon, we save a new copy to "study_area_v01.shp" an so on.
All the attributes in table are the same, nothing changes, but the boundary.

Every time we open all the mxds and change the data source, one by one because none of pythons I found, work fine.

This script works fine, and change all the layerfiles in all MXD in a folder.
But there is a little detail in my case: my layer name in the table of contents have a word with acccent (Área de Estudo) and the script doesn't recognize it. I translated all the script for all of you to understand, but I'm from Brazil and we have a lot of words with accent, and the arcpy doesn't understand it.

Can you help me?

import arcpy, os
inputPath = os.curdir
outputPath = os.curdir

#layerfile name we want to change, in table of contents
layerfile = "Área de Estudo"

#Folder path of the new shapefile
FolderPath = r"T:\.......\Utm"

#Name of the new shapefile, in folder path (without '.shp')
NewShapefile = "AreasEstudo_PDETrevo_v00"

#Loop through each MXD file
for filename in os.listdir(inputPath):
    fullpath = os.path.join(inputPath, filename)
    if os.path.isfile(fullpath):
        if filename.lower().endswith(".mxd"):

            mxd = arcpy.mapping.MapDocument(fullpath)
            df = arcpy.mapping.ListDataFrames(mxd)[0]
            for lyr in arcpy.mapping.ListLayers(mxd, "*", df):
                if lyr.name == layerfile:
                        print layerfile
                        lyr.replaceDataSource(FolderPath, "SHAPEFILE_WORKSPACE", NewShapefile)
            print fullpath
            mxd.save()


texts: Plotting a route from a list of Points without access to ArcGIS Online?

texts:
On ArcMap 10.5.1 with a basic license. I have several excel files I joined with a municipal address layer to create a map to visualize the area for workers.

I am now trying to figure out how best to create a "recommended" route between these points without having any access to the tools I imagine I should be using for this. I cannot access any of the online portions of ArcMap's routing capabilities.

Within the map I'm using I have one way streets labelled, and am using a municipal parcel for the street visualization. Would there be a workflow to get this done given my limitations?


texts: Updating bounding box of layer using geoserver API

texts:
Is there a way to set the Bounding Box of a geoserver layer using the REST Api? I get the BB extent using a postgis query and then I could pass it to a json or xml and post it. But I haven't found any examples online.

I read that gsconfig doesn't support this functionality.
My geoserver version is rather old (2.5).


texts: NDVI values obtained with different sensors

texts: In a research paper with NDVI values, differences were found for the same pixel between Landsat 5, Landsat 8, MODIS or Sentinel made after the corresponding atmospheric corrections. Can this be right? If so, is there any way to correct it?


texts:
I have Climate Data grids_germany_multi_annual_air_temp_mean_1981-2010_17.asc.gz from the ftp-server of the Deutsche Wetterdienst as raster files of multi-annual mean temperature.
Source:
https://opendata.dwd.de/climate_environment/CDC/grids_germany/multi_annual/air_temperature_mean/

After unzipping, I want to change the size of the raster in QGIS. I'm aware that the smaller rasters will have the same value as the original ones. But I'll need a higher resolution for clipping later on.
I tried several things, but I can't figure out how to change the raster size.


texts: QGIS: Is there a way to make the "Python Console" a separate window?

texts:
QGIS 3.10; Ubuntu 18.04
When Plugins->PythonConsole is accessed, it opens as a panel, changing the shape of the map canvas.
Can it be a separate window instead?
Alternatively, is there a way to run a python command; e.g., size() and resize(), without having to open the console?


texts: Modify Service Layer Credits in ArcGIS Pro

texts:
Is there a way to modify (edit, move, resize, delete) Service Layer credits in ArcGIS Pro?

From How to remove Esri basemap credits? I know that it is possible in ArcMap but I could not find a way to do it in ArcGIS Pro.

The default placement is often sub-optimal and modification would alleviate most issues.


texts: Multi-labeling with multiple colours in QGIS

texts: I would like to label data entries with multiple fields (that is easy, using the expression editor) in QGIS 2.18. But I fail to alter the styling of the labels, so that say the colour of the label from one data field (e.g., in the first line) was different to the second label value (second line). An older answer had the standard solution (duplicating the layer, Two different colours in label) I am wondering whether QGIS 2.18 or even the up and coming 3.0 will have a solution.


texts: Spatial adjustment for vector shp files in QGIS

texts:
How to do spatial adjustment for vector shp files in QGIS?

I have a set of 4 cadastral maps of adjacent locations. I have converted them form DWG to DXF and then to vector shapefiles. They are not geo-referenced/ geo-rectified? I want them to be linked to a co-ordinate system (geo-referenced/ spatial adjustment) so that I can join the four parts together. I also have a shapefile of gps points corresponding to one of the four points.

Can you please tell me how to do this in QGIS or GRASS, not ArcMap?

Thanks,
New_QGIS_user


texts: QGIS 3.0.3 takes a lot of time to save selected features

texts:
I'm editing a large .xml layer with over million points and 2GB of data.

I'm trying to save selected features (about 100 000 points).

After an hour of saving, I have just 77 kb of data saved. With this pace, it'd take weeks to have a full file.

Any idea why it takes so long?


texts: Python Script to Identify and Flip Lines Digitized in the Wrong Direction?

texts: I have a river dataset that was digitized both in the direction of flow and opposite the direction of flow. I need to find an automated method to determine which lines were digitized from downstream to upstream and flip them automatically. I do not have spatial, 3D or network analyst. I am looking for a python code. These lines were created in MicroStation using a bare earth LiDAR dataset. They were hand digitized and it is not a stream centerline but a line following the edges of the river.


texts:
In my buffer plugin, i tried the following code

# Get buffer distance and selected city
dist = float( myLineEdit.text() )
city = myComboBo.currentText()

# Get the feature corresponding to the selected city
expr = QgsExpression('"Cities" = \'' + city + '\'' )
it = myLayer.getFeatures( QgsFeatureRequest( expr ) )
feature = it.next()

# Perform a buffer on the feature geometry
myBufferPolygon = feature.geometry().buffer(dist, 50)

And i had also attached my code which may helps in debugging the error

    # Create the dialog (after translation) and keep reference
    self.dlg = city_searcher_classDialog()

    # Declare instance attributes
    self.actions = []
    self.menu = self.tr(u'&City_searcher_plugin')
    # TODO: We are going to let the user set this up in a future iteration
    self.toolbar = self.iface.addToolBar(u'city_searcher_class')
    self.toolbar.setObjectName(u'city_searcher_class')
    QObject.connect(self.dlg.Statecombo, SIGNAL("currentIndexChanged(int)"), self.findTheCityNames)
    QObject.connect(self.dlg.Citycombo, SIGNAL("currentIndexChanged(int)"), self.calcDistance)
    QObject.connect(self.button_box, QtCore.SIGNAL(_fromUtf8("accepted()")), self.buffer)
    #QObject.connect(self.dlg.browsepushButton, SIGNAL("clicked()"), self.saveTheFile)

 result = self.dlg.exec_()
    # See if OK was pressed
    #if result:
    if result:
        QObject.connect(self.button_box, QtCore.SIGNAL(_fromUtf8("accepted()")), self.buffer)

And my buffer function is

 def buffer(self):
    distance = float( self.dlg.DistancelineEdit.text() )
    #new_distance = distance
    #number = self.dlg.number_buf.value()
    #QgsGeometry.buffer(new_distance , 2)
    city = Citycombo.currentText()
    # Get the feature corresponding to the selected city
    expr = QgsExpression('"name_1" = \'' + city + '\'' )
    it = myLayer.getFeatures( QgsFeatureRequest( expr ) )
    feature = it.next()

    # Perform a buffer on the feature geometry
    myBufferPolygon = feature.geometry().buffer(dist, 50)`

And the error i am getting is attached below

And FYI button_box is the name of the OK button in QT Designer.Sorry for the lengthy post.


texts: Select features in statement CASE ELSE or CASE in QGIS 2.10.2

texts:
I need to label only US states and canadian provincies

Here is my attribute table

and my expression

Program show that is no error here, but it still don't work. What I'm doing wrong?


texts: Using "right()" and "strpos()" to get rightpart of string for QGIS Field Calculator?

texts:
The strpos() function searches from the left side of the string. But instead of always trimming the string, I would like to get the position from the right side. Is there anyway to do this in the field calculator?

Here is the specific case of application:

I wanted to get the part of the string in "FolderPath" up the first '/' from the right.


texts:
What is the expected QGIS 2.2 release date?

It would be helpful to schedule my custom development in QGIS for "Relations" management.


texts: Select Feature Based on Iterate Row Value in ModelBuilder

texts:
I have a table showing commuting trips from home to work. home_id represents home and work_id represents work. Now I want to step through the trips table using Iterate Row Selection in ModelBuilder, and get the work_id value. Then use the work_id value to select the corresponding feature in the work location feature class at each iteration.

I used Get Field Value tool in ModelBuilder to get the work_id value. But how can I pass the work_id value to the Select Layer by Attribute tool to use in the SQL expression?

I want something like work_location_id = %work_id%, but this one gives the "ERROR 00358: Invalid expression" error. What is the correct way to do this?


texts: Creating index field from number of items separated by "|" in field

texts:
I am trying to create an index field from a intersection layer I have created using the SAGA, polygon self intersection tool.
The SAGA tool outputs a field which shows which buffers intersect each other separated by a "|". I would like just a simple count of the number of intersections - so I need a way to count the number of numbers between the "|" in that field.


texts: What are the licensing rules for using cartodb's public REST API with leafletjs

texts:
I am currently working on a project which needs some offline mapping. I am using Leaflet.js as my base library. I have also seen that leaflet can be used with the following tile providers:

https://leaflet-extras.github.io/leaflet-providers/preview/

The idea would be to download a minimal set of tiles from one of the sources given in the previous link in order to use them in "offline mode" (no internet connection).

I was initially planning on using cartodb because their base tiles are "clutter-free" and a good starting point to add my own annotations.

Depending on the zoom level, I would need between 1000-10000 tiles that, I will have to update every 6 months.

I am trying to figure out if I am not breaking any terms of services by accessing their public REST API (http://{s}.basemaps.cartocdn.com/light_nolabels/{z}/{x}/{y}.png) in order to download tiles for offline viewing.

There are no API keys or any other form of identification in order to access the afore-mentioned service. Can assume that this is therefore public data (ie free)?

The cartodb terms of service are given here: https://carto.com/terms/#acceptable-use

I am not lawyer myself, but from my understanding it seems to me that it is OK to use their public REST API in order to cache tiles offline as long as I provide attribution.

Can anyone advise on this or tell how they solved a similar issue?

Edit:
OK answering my own question (I can't add a separate answer because my question has been put on hold). Using the carto "base maps/tiles" is free as long as you are not building any commercial apps.

If you are building a commercial app which embeds the base maps, or if you are using the maps in order to track assets (my case), you need a commercial API key. All clearly explained here:

https://carto.com/location-data-services/basemaps/


texts: Where to get osm history extracts in 2019?

texts:
I was trying to make an animation of my city growing, but stumbled upon lack of data. In 2018 geofabrik stopped providing history extracts (osh.pbf), and that was the last source of history files, as I was able to find, all other links on the internet were dead long before.

So, where does one get cropped osm data with history in 2019?

Is processing a full history dump to extract local area a viable option for a modern desktop PC?


texts: I want to calculate the distance between two points using network analysis

texts:
function drawroute(allcordinate) {
        var Coordinate = new Array();
        Coordinate = allcordinate;
        var long = Coordinate.split(",")[0];
        var lat = Coordinate.split(",")[1];
        var long1 = Coordinate.split(",")[2];
        var lat1 = Coordinate.split(",")[3];
        routeTask = new esri.tasks.RouteTask("http://111.93.174.107:6080/arcgis/rest/services/network/KendraparaNetwork/NAServer/Route");
        //set up the route parameters
        routeParams = new esri.tasks.RouteParameters();
        routeParams.stops = new esri.tasks.FeatureSet();
        routeParams.outSpatialReference = { "wkid": 102100 };
        dojo.connect(routeTask, "onSolveComplete", showRoute);
        dojo.connect(routeTask, "onError", errorHandler);
        //define the symbology used to display the route
        stopSymbol = new esri.symbol.SimpleMarkerSymbol().setStyle(esri.symbol.SimpleMarkerSymbol.STYLE_CROSS).setSize(15);
        stopSymbol.outline.setWidth(4);
        routeSymbol = new esri.symbol.SimpleLineSymbol().setColor(new dojo.Color([246, 26, 26, 0.7])).setWidth(5);
        var stop;
        var start = new esri.geometry.Point(lat, long, new esri.SpatialReference({ wkid: 4326 }));
        var end = new esri.geometry.Point(lat1, long1, new esri.SpatialReference({ wkid: 4326 }));
        for (var i = 0; i < 2; i++) {
            if (i == 0) {
                stop = map.graphics.add(new esri.Graphic(start, stopSymbol));
                routeParams.stops.features.push(stop);
            }
            else {
                stop = map.graphics.add(new esri.Graphic(end, stopSymbol));
                routeParams.stops.features.push(stop);
            }
        }
        if (routeParams.stops.features.length >= 2) {

            routeTask.solve(routeParams);
            lastStop = routeParams.stops.features.splice(0, 1)[0];
        }
    }


texts:
I have created a Mapbook with each page containing a section of the main map and an index map.  They both seem to be georeferenced when exported in a PDF.

When a user with an Apple product uses the map with Avenza PDF, the current position, tracks, and place marks are directed to the index map, not the main map.

Is there a way to georeference ONE and NOT both of the maps?


texts: How to install missing resampling methods in GDAL Warp (QGIS)?

texts: I am trying to resample my data based on "sum" resampling method. But unfortunately, I do not see this in the drop down menu in Warp.I am running a GDAl version 3.1, but it does not have this function. Can someone please help me figure out how to "install" this?


texts: Where can I find the nightly build release number?

texts: As far as I know every QGIS nightly build has a release number such as 81be559? Where can I find this number?


texts: Adding hyperlinks to Leaflet polylines created via Leaflet.Draw?

texts:
[Cross posting from StackOverflow as I realized this may be a more specific community for this]

I'm using Leaflet.Draw, and I want to allow the lines that I draw to have hyperlinks in them so that someone can click on the link and open a new page.

How do I add a link to a polyline created via Leaflet.Draw? When I addLayer to drawnItems, how can I give that layer another attribute that could be a URL?

var map = L.map('map').setView([-41.2858, 174.78682], 14);
mapLink =
    'OpenStreetMap';
L.tileLayer(
    'http://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png', {
    attribution: '&copy; ' + mapLink + ' Contributors',
    maxZoom: 18,
    }).addTo(map);

var drawnItems = new L.FeatureGroup();
map.addLayer(drawnItems);

var drawControl = new L.Control.Draw({
    edit: {
        featureGroup: drawnItems
    }
});
map.addControl(drawControl);

map.on('draw:created', function (e) {
    var type = e.layerType,
        layer = e.layer;
    drawnItems.addLayer(layer);
});


texts:
I wonder if it is possible to load German NAS-Data into QGIS or to convert them to Shapefiles. I'm using QGIS 1.8 on Ubuntu 12.04.

NAS is the new Geodata-format, based on ALKIS-Data:
http://de.wikipedia.org/wiki/Normbasierte_Austauschschnittstelle

I installed the current Master on my Ubuntu 12.04 machine today and typed "ogr2org --formats" but NAS didn't appear on the list.

If I check my ogr-version it says 1.7.3 - as I understood NAS is supported from 1.8.0 onwards.

Did you do anything special to install that version?


texts:
I have a .shp file with includes a column of format Double that is the Arrival time.
The time values entered are decimal time. For instance, 1 record has a decimal time of 7.11194, which converts to 07:06:43
What is the easiest method to convert this value so I can use it in QGIS?
I've tried using the totime value in QGIS, but no luck.


texts: How to fix white tile issue when exporting Sentinel-2 Satellite images as video to drive

texts:
I'm trying to retrieve Sentinel-2 satellite images from Google Earth Engine using Python, and export them as a video to google drive. But the problem I'm facing is that regardless of the export settings, the video ends up as completely black with a white tile that moves every second.
I've also tried changing the coordinates, filters and bands, but it still gives me the same result. Any idea why this might be the case?
Images below for reference:

This is the code:

import ee
import time
from ee import batch

# Authenticating and connecting to Drive
# Fetching the Sentinel Image Collection

## Authenticating
ee.Authenticate()

## Initialize
ee.Initialize()

## define your collection
collectionSentinel = ee.ImageCollection('COPERNICUS/S2').filterDate('2015-06-23', '2021-03-18')
print('collection done')

coodinates = [[66.79667713688107,24.885210852891888], [67.31749775455685,24.885210852891888], [67.31749775455685,25.185855686990557],  [66.79667713688107,25.185855686990557],  [66.79667713688107,24.885210852891888]]
print("coordinates received")

geoRegion = ee.Geometry.Polygon([coordinates])
pathrowSentinel = collectionSentinel.filterBounds(geoRegion)
print('region done')

##Filter cloudy scenes.
cloudsSentinel = pathrowSentinel.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 30))
print('clouds filtered')

## select the bands, we are going for true colour... but could be any!
bandsSentinel = cloudsSentinel.select(['B4', 'B3', 'B2'])
print('bands selected')

##make the data 8-bit.
def convertBit(image):
    return image.multiply(512).uint8()

## call the conversion
outputVideoSentinel = bandsSentinel.map(convertBit)
print('converted')

#Export to video.
print("about to build Sentinel video")
outSentinel = batch.Export.video.toDrive(outputVideoSentinel, description='khi_video_sentinel', dimensions = 720, framesPerSecond = 12, region=geoRegion, maxFrames=10000)
outSentinel.start()
print("process sentinel sent to cloud")


texts:
I'm stuck with a problem:
I got three raster's layers: "Sand", "Clay" and "BD" (Bulk density) of the same region.
How can I obtain a "table" or a "Data frame" with information for the values of the raster's for each pixel?
I prefer to use QGIS or R.
My main goal is to obtain a lookuptable of this features of the soil:

So i will merge pixels with the same information for Clay, Sand and BD into a unique MU_GLOBAL index.


texts: How to convert post gis database into shapefile, I found that I need to use pgsql2shp command on postgis database query , But how to use it and where to use it. Do I have to use it in the postgre sql database command line or in the normal command line?


texts: How to add map.addControl, WMSLegend in OpenLayers 2

texts:
I am trying to get the categorized symbology for a single WMS layer in the layer switcher.

See my code below (I am getting error OpenLayers.Control.WMSLegend is not a constructor)

 My OpenLayers Map

//OpenLayers.ProxyHost = "/cgi-bin/proxy.cgi?url=";

var map;

function init() {
    map = new OpenLayers.Map('ch05_wmsfeatureinfo', {});

// Setup our two layer objects
            var wms_layer_map = new OpenLayers.Layer.WMS(
                'Base layer',
                'http://vmap0.tiles.osgeo.org/wms/vmap0',
                {layers: 'basic'},
                {isBaseLayer: true}
            );

        var okombahe_town_and_townlands = new OpenLayers.Layer.WMS(
            'Townland_Boundary',
            'http://192.168.0.112:8080/geoserver/Erongo/wms',
            {layers: 'Erongo:okombahe_town_and_townlands',
            transparent: true,
            isBaseLayer: false},
            {visibility: false, opacity:0.5}
        );

        var okombahe_extensions = new OpenLayers.Layer.WMS(
            'Extension_Boundaries',
            'http://192.168.0.112:8080/geoserver/Erongo/wms',
            {layers: 'Erongo:okombahe_extensions',
            transparent: true,
            isBaseLayer: false},
            {visibility: false, opacity:0.5}
        );

        var okombahe_erven = new OpenLayers.Layer.WMS(
            'Okombahe Erven',
            'http://192.168.0.112:8080/geoserver/Erongo/wms',
            {layers: 'Erongo:okombahe_erven',
            transparent: true,
            isBaseLayer: false},
            {visibility: false, opacity:0.5}
            );

                    map.addLayers([
                    wms_layer_map,
                    okombahe_town_and_townlands,
                    okombahe_extensions,
                    okombahe_erven,
                    ]);

map.addControl(new OpenLayers.Control.LayerSwitcher({}));
map.addControl(new OpenLayers.Control.WMSLegend({
class: 'ol_legend'

}));

map.setCenter(new OpenLayers.LonLat(30, -30), 4);

OpenLayers.Control.WMSLegend = OpenLayers.Class(OpenLayers.Control, {
    class: null,
    wmsVersion: '1.3.0',
    format: 'image/png',
    draw: function () {
        var legendP = document.createElement('p');
        legendP.innerHTML = 'Legend:';
        this.div = document.createElement('div');
        this.div.className = this.class + ' olControlNoSelect';
        this.div.appendChild(legendP);
        var layers = this.map.layers;
        for(var i=0;i

info = new OpenLayers.Control.WMSGetFeatureInfo({
        url: 'http://192.168.0.112:8080/geoserver/Erongo/wms',
        title: 'Identify features by clicking',
        queryVisible: true,
        eventListeners: {
            getfeatureinfo: function(event) {
                if (event.text.length>658) // <<=== use this, empty feature return 658 character or just check the length using alert(event.text.length);

            {

            var pops = map.popups;
             for (var a = 0; a < pops.length; a++) {
                 if (pops.length == 1) {
                     map.removePopup(map.popups[0]);
                 }
             };

            map.addPopup(new OpenLayers.Popup.FramedCloud(
                "ayam",
                map.getLonLatFromPixel(event.xy),
                null,
                event.text,
                null,
                true
            ));
            }
        }
    }
});

 map.addControl(info);
info.activate();
}

Activated


texts: About geodesic or azimuth and UTM and Google Mercator projections

texts:
I need to calculate to azimuth between two points.
I'm using Projnet and GeographicLib to get the azimuth. So I have a UTM projection map then I need Projnet to get Lat/Lon coordinates because GeographicLib uses Lat/Lon to calculate this azimuth.
Ok I get the angle but now my question is if I need to reproject this azimuth to my UTM projection map???


texts: find number of points closer to lines inside a polygon ArcGIS

texts: I have a map that contains 3 shapefiles representing projects (points), studies performed (lines) and sectors (polygons). I would like to know how can I find the sector (polygon) that contains the higher number of projects (points) closer to the studies performed (lines). See image.


texts: QgsRasterBandStats returns nan for mean and sum of squares in QGIS3.4

texts:
I have a problem with certain raster files in QGIS 3.4.
When requesting statistics via QgsRasterBandStats the average, the sum returns nan.
The other statistics minimum, maximum, sum of squares and standard deviation if they return expected values, float.
I think, that the error could be that the method considers the nan in the calculation, in the mean and sum.
How can I solve it?
I want to avoid the alternative of reading the raster with Gdal again directly. In that case it would be better to redo all my code and work directly with GDal and numpy
An example of the code:
capa=iface.activeLayer()
provider=capa.dataProvider() `
stat=capa.`dataProvider().bandStatistics(6,QgsRasterBandStats.All)
print(stat.sum)
   nan
print(stat.mean)
   nan
print(stat.elementCount)
   70848
print(stat.stdDev)
   0.02975916513582063


texts:
I always get ERROR: Unable To Seek after I run r.resamp.rst

r.resamp.rst input=map30 ew_res=10m ns_res=10m elevation=map10

It actually starts to work but after it finish I get the error and no output.

Fresh install of Grass GIS 7

I'm suspecting this has something to do with invalid input files, and resolutions.
So the question is:

What projections and input file parameters are supported by r.resamp.rst?


texts:
I'm storing a large list of lnt/lat points in a MySQL database. At the moment these are projected points for the UK in M, but in the long term I would like to make sure that I can store the coordinates of points for all over the world.  What datatype should I use?

I started using decimal(18,12), but was unsure whether this precision is needed or whether I could just use a float. I included my code, in case there is anything else I should be considering:

CREATE TABLE UKTest
(
lat FLOAT,
lng FLOAT
)

I'm quite new to SQL so I'm unsure whether the data type is significant when running a query. I assumed that the overall memory would vary based on the variable type. Is there a benefit of using a spatial database over a basic MySQL database for this type of work?


texts: Should I start a new project using MapGuide OS?

texts:
I'm wary of starting a debate based on opinion, so I lifted my title from this question which seems to be doing fine. I'm also worried about offending MapGuide enthusiasts so please understand I am simply trying to avoid unwise strategic decisions.

When I try to install, use, develop with, ask questions about, or otherwise entertain MapGuide OS I get the feeling that starting a new project with MapGuide is too risky:

API documentation has not been updated for years and developers are asked to use 2.0 docs for the 2.2 and 2.4 product
It took several months from the time 2.2.0 was released for http://mapguide.osgeo.org/ to recognise that 2.1.0 was not the latest stable release, suggesting that no one is actively maintaining the site
Install documentation is spread wide and thin. I recently spent the best part of a week trying to install on 64-bit linux for use with Oracle and ultimately had to give up (then installed GeoServer in less than one hour)
Forum threads regularly end unanswered or stumble along with hopeful suggestions only to end with an unanswered question from the OP
I personally have never managed to get a .NET map viewer up and running due to compile errors that others have reported but had no luck fixing
It seems like a lot of the effort that goes into MapGuide comes from a very small number of people (especially Jackie Ng) which doesn't bode well should they ever get too busy or lose interest in the project.

Someone recently asked a similar question on the MapGuide forum but it took a decidedly ambiguous tone and ended with no clear answer.

I appreciate that a lot of work has gone into the project, and also that it provides functionality I can't find elsewhere (this is a big part of the problem I have), as it can provide a complete viewer experience in short order. However in addition to the list above I have other issues:

Dependency on stateful sessions and associated load balancing considerations
Requirement for a configuration tool or web service to create opaque configuration files / databases that can't be edited or viewed easily
PHP web viewer's need for a specific version of PHP which is unlikely to be compatible / ideal for other PHP applications on a shared server
Distinction between IE and other browsers in the viewers due to ActiveX

Yes, I know that this is an Open Source project and so if I think it needs improving I should get involved, but honestly I would much rather donate my time to other projects like GeoServer or OpenLayers as these products have both made my life much easier and have well defined roadmaps.


texts: How to transform a shapefile's projection in ArcGIS?

texts:
I am working with a series of shapefiles of France, which were created privately by someone (as opposed to an official institution). We can call them the P-maps. There is no .prj file associated with them. When I open them in QGIS or ArcMap, and set their coordinate system to be WGS84, they "look alright". However, they don't overlay with a correct shapefile of France in WGS84 (a map in WGS84 from the IGN, let's call it IGN) , so another projection was used (you can see the picture, the IGN shapefile is circled in blue; one of the P-maps is pointed in red).

I know I could potentially try to identify the correct projection by trying each individual projection one by one, by I am weary of using this method, as I am not sure that this person correctly georeferenced these P-maps. I think that they establish them by creating the shapefile from an image that was not georeferenced.

Is there a way to define a transformation to each file (I have around 50 P-maps, and they all align to each other) so that  GIS would define the right coordinates systems, having the comparison point of the IGN shapefile in the WGS84 projection? Something like when you georeference a raster map. Since all the  maps align, I could apply this transformation to each one of them.

The layer extent:
top 590.9218000
right 620.0
bottom 30.0
left 20.02860

And the file you can download


texts:
I'm trying to implement this function to create a hex grid over a layer.

However, when I implement it and export it via pgsql2shp the results ends up empty. The query I used is below

CREATE TABLE hex_grid (gid serial not null primary key);
SELECT addgeometrycolumn('hex_grid','the_geom', 0, 'POLYGON', 2);

CREATE OR REPLACE FUNCTION genhexagons(width float, xmin float,ymin  float,xmax float,ymax float  )
RETURNS float AS $total$
declare
    b float :=width/2;
    a float :=b/2; --sin(30)=.5
    c float :=2*a;
    height float := 2*a+c;  --1.1547*width;
    ncol float :=ceil(abs(xmax-xmin)/width);
    nrow float :=ceil(abs(ymax-ymin)/height);

    polygon_string varchar := 'POLYGON((' ||
                                        0 || ' ' || 0     || ' , ' ||
                                        b || ' ' || a     || ' , ' ||
                                        b || ' ' || a+c   || ' , ' ||
                                        0 || ' ' || a+c+a || ' , ' ||
                                     -1*b || ' ' || a+c   || ' , ' ||
                                     -1*b || ' ' || a     || ' , ' ||
                                        0 || ' ' || 0     ||
                                '))';
BEGIN
    INSERT INTO hex_grid (the_geom) SELECT st_translate(the_geom, x_series*(2*a+c)+xmin, y_series*(2*(c+a))+ymin)
    from generate_series(0, ncol::int , 1) as x_series,
    generate_series(0, nrow::int,1 ) as y_series,
    (
       SELECT polygon_string::geometry as the_geom
       UNION
       SELECT ST_Translate(polygon_string::geometry, b , a+c)  as the_geom
    ) as two_hex;
    ALTER TABLE hex_grid
    ALTER COLUMN the_geom TYPE geometry(Polygon, 4326)
    USING ST_SetSRID(the_geom,4326);
    RETURN NULL;
END;
$total$ LANGUAGE plpgsql;

--width in the units of the projection, xmin,ymin,xmax,ymax
SELECT genhexagons(1.0, -124.8, 25.10, -66.9, 49.6);

And then I export it using the pgsql2shp pgsql2shp -f test -h wideload -u minhmai -P password! test_nd_2016g "select * from hex_grid"

But when I open it in QGIS, the shapefile is blank. I'm not sure about what the function really does but I assume it should have something when exported, especially since the output indicated that there are rows.


texts: Graphical Modeler to rotate points based on line layer azimuth

texts:
I am trying to create a model using the QGIS 2.0 Processing Toolbox graphical model to replicate the process described in this thread: How to horizontally align irregularly shaped point symbol to vector layer in QGIS?

However, I get an error when trying to execute the Distance to Nearest Hub algorithm: Invalid name attribute: bearing. When I go through the process manually, using the same data, this doesn't happen.

UPDATE: running the field calculator using the python console doesn't populate my vector layer's bearing field as I expected it to. A field called "bearing" is created in a new output layer (i've corrected to link this to distance to nearest hub instead of how it is currently), but the field is blank.

I ran the script like so after import processing:

processing.runalg("qgis:fieldcalculator", resLayer, 'bearing', 0, '(atan((xat(-1)-xat(0))/(yat(-1)-yat(0)))) * 180/3.14159 + (180 *(((yat(-1)-yat(0)) < 0) + (((xat(-1)-xat(0)) < 0 AND (yat(-1) - yat(0)) >0)*2)))', 'bearingdata.shp')


texts: Setting Layer in the Georeferencing ToolBar in arcMap 10?

texts:
Does anyone out there knows how to set the layer in the georeferencing toolbar. There is no information in the object models. Using
ICommandItem ic =   ArcMap.Application.Document.CommandBars.Find(uid, false, false) as ICommandItem;
I found that type is georeferencing_Layerslistcontrol. Any ideas?


texts: Displaying symbols correctly in legend in ArcGIS Desktop after PDF Export?

texts: I am trying to create some maps in which have some point symbols. when I export the layout to pdf, the symbols are not appear. I have read some similar questions in the group but were not helpful for my case. there is an example of the mentioned problem in below.


texts: GeoPandas/Geoplot losing first row when reprojecting

texts:
I'm trying to work with the London Borough borders file (https://data.london.gov.uk/dataset/statistical-gis-boundary-files-london) in geopandas/geoplot.
It comes in National Grid Easting/Northing coordinates (epsg="7405" I think) and I'm trying to convert to Lat Longs to be able to map, but the first row of any dataframe I convert keeps being changed to an "INF INF" polygon!
No idea why it's happening.  I've tried doing different rows, only one row, changing the CRS....same issue.
all_border = geopandas.read_file("boroughs\London_Borough_Excluding_MHW.tab")

lsoa_border.to_crs(epsg=4327)


texts: Permissions for QGIS 3 project files on Windows 10

texts:
When saving a project under QGIS 3.x on a local Windows PC under e.g. c:\tmp another user logged on to the PC cannot open the project due to lacking permissions. The project file saved from QGIS 3 permissiones are restricted to groups SYSTEM and Admin and for the current user:

Project files saved under QGIS 2.x have in addition full access for groups User and Everyone:

Changing permissions for the target folder (c:\tmp) has no affect. Remarkable: This behavior does not apply to i.e. Layer definitions saved from QGIS 3, these have the same permissions as project files saved from QGIS 2 (cp.above screenshot).
Is this a bug or a feature or am I missing any settings?


texts: Field Calculator: Prefixed, padded IDs from a number field

texts:
Using the Field Calculator, I want to create prefixed, padded IDs from an existing number field. By "padded", I mean I want the IDs to have a specific number of characters, regardless of how many digits the numbers have.

For example, I want to use the values in OBJECTID to create IDs like:

OBJECTID 8 >>> INSP000008
OBJECTID 9  >>> INSP000009
OBJECTID 10 >> INSP000010
OBJECTID 11 >> INSP000011

Note: All the resulting IDs have the same number of characters (10).

How can I do this?


texts: Using Google Maps to track work per building and allowing users to manipulate that data

texts:
I work part time for a political campaign office, who are going door-to-door for several months at a time, all over a city. We need a system to track where we have been, and how long ago we were last there - to be able to know where to visit each week, etc.

Now, I've never done any work with the Google Maps API or MyMaps before, but it does seem like a fairly easy platform for that kind of project. However, my research hasn't given me much to work with. Is it possible to create a system like this using Google Maps, or would I need to build my own application?

The features I'd need are:

Placing permanent polygons on top of a map surface, covering individual buildings or blocks.
Assigning those polygons with a value, indicating when those buildings/blocks were last visited. Perhaps this could be accomplished using color marking?
Allowing users to modify the "last visited"-variable for each polygon, without allowing them to move or modify the polygons themselves.

This will be an internal system used only by staff, so I'm not worried about security. The access limits are only necessary for making sure that no-one accidentally removes polygons or "last visited"-data.

Any advice to get me in the right direction would be great. Thank you!


texts: Access is Denied in cmd prompt shp2pgsql

texts:
I was wondering if anyone has come across the error access is denied in command prompt when wanting to load a shapefile using shp2pgsql for postgresql.

I have set my user access security level to bare minimums,  given every single user/admin etc.. full rights to the shapefile folder and command prompt.  I am puzzled to what it could be.  I have fairly heavy firewalls but I don’t think that should be preventing access.


texts: ArcScene: creating 3D polygon features from a DEM

texts:
I have a DEM with a lot of topographic variance and a separate layer of polygon features (w/o Z values)throughout. When I use 3D Analyst to convert features(polygon) to 3D, only a portion of the polygons are overlaid. Many of the polygons seem to "be inside the DEM" or underneath, specifically the polygons that appear to be on a hillside. I've tried several extrusion methods, but can not find a clear visualization of the data.
Any recommendations are welcomed.


texts: TypeError: 'tuple' object is not callable (QGIS 3)

texts:
I'm getting this error when trying to run the following code:

line = iface.activeLayer()

a = processing.run("qgis:fieldcalculator", {'INPUT':line,'FIELD_NAME':'len','FIELD_TYPE':1,'FIELD_LENGTH':10,'FIELD_PRECISION':3,'NEW_FIELD':True,'FORMULA':' $length ','OUTPUT':'memory:'})
layer1 =a['OUTPUT']
featparede = layer1.getFeatures()
parede = []
for f in featparede:
    parede.append(f['len'])
print (parede)

Now I have shutted down the Pc and restarted it. The code worked once but in the second time it gave the error. Now it only gives the error.

This is the complete error:

Traceback (most recent call last):
  File "C:\PROGRA~1\QGIS3~1.6\apps\Python37\lib\code.py", line 90, in runcode
    exec(code, self.locals)
  File "<input>", line 1, in <module>
  File "<string>", line 8, in <module>
TypeError: 'tuple' object is not callable


texts: Is it possible to apply a gaussian filter to rasters in ArcMap without applying a calculation or otherwise program it yourself?


texts: Install postgis 2.3.0 for postgres 9.5.3.1 on Win7 in unattended mode?

texts:
Does any one know how to do this?

The PostgreSQL install is easy and works like a charm but I'm trying to also install PostGIS in unattended mode with no user interaction, and it keeps popping up the gui and asking me to agree with the terms.


texts:
I have a layer that is published by GeoServer (2.12.0) and then cached by MapProxy (1.12.0). The layer is a layergroup on the server consisting of a polygon and a GeoTIFF showing a proximity gradient how the layer looks like delivered by the GeoServer is like this:

Nice and smooth gradient but when I open the layer on my mapproxy it looks like this:

What I think is happening is that MapProxy upscales the image from a higher zoom level. Why this does happen I don't know.
The merge.yaml entry for the layer looks like this:
sources:
  layer_name_source:
    type: wms
    http:
      client_timeout: 5 # seconds
    req:
      << : *geoserver_url
      layers:   'layer_name_on_geoserver'
      transparent: true

caches:
  layer_name_cache:
    cache:
      type: file
      directory: /cache/dir/
      directory_layout: tms
    sources: [layer_name_source]
    grids: [GLOBAL_WEBMERCATOR]

layers:
  - name: layer_name_layer
    title: layer_name
    sources: [layer_name_cache]


texts: Finding format of data source used in official Carto example

texts:
I'm trying to apply my own data to this example, but I'm not sure what format the original data was in. I've checked through Carto's stock datasets in the web dashboard, but I can't find it.

Is there somewhere I could find the bird_journey data used in the Carto training guide, or is there a way I could somehow access it through javascript using sources?


texts: Chinese characters appear in attribute table

texts: I am running QGIS 2.4 on a mac, and all of the text in my attribute tables appears in Chinese.  I've tried uninstalling and re-installing, changing the locale to U.S. English, and nothing helps.  I've tried opening attribute tables in other applications, and they appear in English everywhere else, so the problem seems to be in importing them to QGIS.


texts: How to get an array of simple lines from one multipoint linestring in PostGIS?

texts: How can I extract some LINESTRING consisting of 3 or more POINTs to several LINESTRINGs each of 2 POINTs in PostGIS?


texts: Positing Labels around on buffers around a centroid

texts:
I'm using Qgis 2.16.2, I am having an issue positioning labels around a centroid.
To help you to built a mental picture of what i am doing: I have a centroid point with 400 metre buffers so 400, 800, 1200 etc.. I am trying to label each buffer with corresponding values

The issue i have is that whenever i create labels they default to being horizontal (next to each other) along the buffers, which is great but i want  them to be vertical (above each other) mainly for composer purposes. Is there anyway i can change the positioning of the labels without manually moving them and just using the Label positioning tools? or at least a way that when i manually move them to ensure they are all perfectly lined up?


texts:
I have a Python expression that I use in Field Calculator.
!PROPR_NAME!.split(",")[0]

I am writing a script in PyCharm and I am wanting to add this expression, however I am uncertain of how to re-write it in ArcPy code.
Does anyone know the code for split?


texts: Not all tiles loading in view (web) with Mapbox?

texts:
For some reason, my map does not load all the tiles when viewed on the web. Here is a link to the test data - http://bestersurveys.co.za/htmlsite/maptest.html

I have nothing funny in my code and really can't find a solution on the net, so can someone help me out?

Source Code:

<!DOCTYPE html>

    <meta charset='utf-8' />
    <title></title>
    <meta name='viewport' content='initial-scale=1,maximum-scale=1,user-scalable=no' />

    <link href='https://api.tiles.mapbox.com/mapbox-gl-js/v0.35.0/mapbox-gl.css' rel='stylesheet' />

<button id='fly'>Fly</button>


texts: Find NODATA values in TIF files using Python or QGIS

texts:
I have got about 7000 tifs in a directory. They all have the same defined NODATA value but not all of them contain cells with NODATA.
How to tell which ones have NODATA cells programmatically? Python is preferred but PyQGIS is also good.


texts: What reliable sources (paid or free) exist for fossil fuel infrastructure data?

texts:
I am trying to find GIS data for various types of petroleum and natural gas infrastructure (pipelines, refineries, gathering lines, pump stations, etc). I've found quite a few sources, mostly for-profit companies selling the data. However, without spending thousands of dollars on the data, it's hard to know what I'd be getting from any of these vendors. Most of them state the resolution of the data but I'm more concerned with the completeness of it. I've also found several government sources, which seem to offer more reliable data but the different layers I'm interested in are scattered across several agencies or just not available. So,

I'm interested in data in the US in this case.

The sources I've found so far include:

Paid: Platts, MAPSearch, HTSI, Rextag

Free: U.S. Energy Information Administration, National Pipeline Mapping System (restricted access)


texts:
I am not sure if it is a right place to ask this one but let me try.

Monitoring programmes are known for various targets and for various spatial levels.
They may include regions, countries or even continents. It is widely accepted that
monitoring systems are of increasing importance as economic activities and their
impacts on nature become more and more global.

I am looking for a various monitoring programmes. The best will be that related to forestry (but in fact i am seeking for anyone) like this one

I am interested
how do you evaluate the selected level I and level II indicators used in the ICP
forest ecosystem monitoring programme?
Do you consider them exhaustive to get a comprehensive view of the state of the
forest ecosystems under the influence of air pollution?

Also if you have some examples in mind it will be great to share it!


texts:
So trying to work through the setup of geoserver and using an example script see below.

<!DOCTYPE html>

    <title>Simple Map</title>
    <link rel="stylesheet" href="http://openlayers.org/en/v3.15.1/css/ol.css" type="text/css">

When opening the page though i cannot see the sample shapefiles 'topp:states'
Anyone know what i am doing wrong?

I assume i have something wrong with the URL for the source.


texts: Convert Ordnance Survey NTF data to QGIS compatible format

texts:
Is there a way of converting Ordnance Survey NTF data to QGIS that will correctly import the data and build polygons?

QGIS will read the data as is, but it imports the polygons only as centrepoints.


texts: Are there any sites that gather free archaeology datasets?


texts: How to make the text color changing within a label in QGIS?

texts:
I would like to keep in black a part of my label (title) and make the other one varying according to a value.

In this (simplified) example, as required, the value color is varying according to a concentration, but also the title.

I would like to keep the first line in black and make the second one varying.

Text formula:

'Drill ' ||  "id"    ||  '\n'  ||
CASE
    WHEN "Conc" >100
        THEN
            "Polutant"  || ': ' ||  "Conc"
        else
        ''
END

Color formula:

if  ("Conc"> 220, color_rgb( 255,0,0) , color_rgb( 0,0,0))

Can QGis handle such a situation and how, or have I to program a function?
In this last case, do you have an example to put me on the right way?


texts:
I'm trying to parse a shapefile with ogr and export the values to WKT. The geometries contain a m value but this is not showing in the WKT output.

Does ogr support 4D/m values?

# My code
from osgeo import ogr
shp = ogr.Open(shapefile.shp)
lyr = shp.GetLayer()
for feature in lyr:
   # output - 'LINESTRING (557729 4970394 1869,557024 4970372 1806, ...)'
   # It should be - 'LINESTRING (557729 4970394 1869 1,557024 4970372 1806 2, ...)'
   print(feature.GetGeometryRef().ExportToWkt())


texts: Can I load a WFS in a deegree Processlet?

texts:
I am very new to deegree and Processlets.
After working through the documentation and implementing both examples I am stuck. I cant find further examples or documentation on the possibilities of Processlets.
More concrete - is it possible to provide a Processlet with the id of a WFS as an input, perform calculations on the feature attributes is return a list of `feature results?
I am also very happy with further resources or examples.
EDIT: I did go through further examples and could not find an example of requesting features from an existing WFS.


texts: How to override QgsMapCanvas default zoom behavior?

texts:
I'm writing a standalone QGIS Python app. I've subclassed QgsMapTool and have successfully implemented some selection capabilities. Now I'm trying to implement my own zoom behavior, using the same map tool, but QgsMapCanvas seems to have some built-in zooming behavior. How can disable the built-in zooming and get complete control over the mouse wheel?

class MyMapTool(QgsMapTool):
    ...
    def wheelEvent(self, event):
        new_extent = QgsRectangle(...)
        self.canvas.setExtent(new_extent)
        # Something else is messing with the extent.
        # How can I disable it and use my new_extent?

# Set the tool to the canvas. This is the only Map Tool I use.
tool = MyMapTool(self.canvas)
self.canvas.setMapTool(tool)

[EDIT]: I'm using QGIS 2.18.


texts: Referencing text from Excel sheet in text box of ArcGIS Desktop?

texts:
I'm trying to reference data in an xlsx document in a text box (Insert --> Text) in ArcGIS.

Is it possible to make it so text in the text box is referenced directly from the xlsx document so it updates when the document is updated?

I'm looking into it, but, I'm having trouble finding this information online.


texts: Converting Lat Long CSV file to meters using ArcGIS Desktop?

texts:
I have a large file in a CSV or excel format and want to convert the Lat longs into Northings and Eastings in order to then put them into ArcGIS Desktop and convert those Northings and Eastings into a custom UTM coordinate system, here is the file.
https://data.edmonton.ca/City-Administration/Property-Information-Data-Current-Calendar-Year-/dkk9-cj3x

How do I change a few or all of them?

I am unable to convert it into excel and put it into ArcGIS Desktop, that doesn't seem to work as the scale is completely off.


texts:
I have never worked with the raster calculator before and I can't solve my problem.

There are two raster datasets. One contains roads ( "closed", "in use", "under construction") and the other dataset contains 4 countries ("1","2","3","4").
I just want to select from roads "in use" and all the countries that contain the roads ...

I am using ArcMap 10.2.2
So I tried something like this:
(in use has the Value 2)

("roads" == 2) & ("countr" == 1)

This works. But how can I get all the countries that contain that road(s)?

My goal is to compute a new output like a "Clip". Cells of roads "in use" and cells from countries should be "cliped".
Is this possible using the Raster Calculator? it looks as if i didnt get the concept of Rasterdata ;-(


texts: How to Import Metadata of feature classes with them into ArcSDE?

texts: I have a database that has numerous feature classes of different types.  None of these feature classes have metadata.  I understand that when you import feature classes to a database, the metadata does not come with.  I'm trying to import metadata to features that have never had metadata and I'm not seeing a way to do this.


texts: Assigning a classification to a shapefile to export to CAD

texts: I'm trying to export a number of polygon shapefiles from Arcmap 10.3 to AutoCAD. The polygons are divided into four different classes and I need to retain the information of which polygon belongs to which class. The symbology attached to the shapefiles does not export with the polygons so is there a method of assigning a classification to a shapefile to denote the class to which it belongs?


texts: Problem doing a coordinate transform with GDAL python

texts:
I'm trying to do a concordant transform from the format of the shapefile (3309) to google maps(4326).

This code works on one Centos 6.6 box with the SCL python 2.7 and a custom gdal 2.0 build. Another similarly configured box gives an error. I seem to have the two datums installed:

$ /opt/gdal-custom/bin/gdalsrsinfo epsg:3309

PROJ.4 : '+proj=aea +lat_1=34 +lat_2=40.5 +lat_0=0 +lon_0=-120 +x_0=0 +y_0=-4000000 +datum=NAD27 +units=m +no_defs '

OGC WKT :
PROJCS["NAD27 / California Albers",
    GEOGCS["NAD27",
        DATUM["North_American_Datum_1927",
            SPHEROID["Clarke 1866",6378206.4,294.9786982138982,
                AUTHORITY["EPSG","7008"]],
            AUTHORITY["EPSG","6267"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4267"]],
    PROJECTION["Albers_Conic_Equal_Area"],
    PARAMETER["standard_parallel_1",34],
    PARAMETER["standard_parallel_2",40.5],
    PARAMETER["latitude_of_center",0],
    PARAMETER["longitude_of_center",-120],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",-4000000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["X",EAST],
    AXIS["Y",NORTH],
    AUTHORITY["EPSG","3309"]]

$ /opt/gdal-custom/bin/gdalsrsinfo epsg:4326

PROJ.4 : '+proj=longlat +datum=WGS84 +no_defs '

OGC WKT :
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4326"]]

Snippet in question:

 def coordinateTransformation(self,EpsgFrom,EpsgTo):
    from_ogrSpatialReference = osr.SpatialReference()
    to_ogrSpatialReference = osr.SpatialReference()
    from_ogrSpatialReference.ImportFromEPSG(EpsgFrom)
    to_ogrSpatialReference.ImportFromEPSG(EpsgTo)
    retObj = osr.CoordinateTransformation(from_ogrSpatialReference,to_ogrSpatialReference)
    from pprint import pprint
    print("CoordinateTransformation: ")
    pprint(retObj)
    return retObj
. . . .

self.coordTransform = self.coordinateTransformation(4326,3309)
. . .
intersectgeometry = featuregeom.Intersection(poly)
from pprint import pprint
print ("Intersection: {0}".format(intersectgeometry.ExportToJson()))
pprint (self.coordTransform)
intersectgeometry.Transform(self.coordTransform)

Output:

<osgeo.osr.CoordinateTransformation; proxy of <Swig Object of type 'OSRCoordinateTransformationShadow *' at 0x7f9f64652120> >
Intersection: { "type": "LineString", "coordinates": [ [ 1286274.293049275642261, -89857.546843856573105, 0.0 ], [ 1284590.659783720737323, -90115.061295109800994, 0.0 ], [ 1279467.874491626396775, -92807.228362500201911, 0.0 ], [ 1273444.6429856531322, -94958.127677664160728, 0.0 ], [ 1270756.154574947431684, -98174.25142274517566, 0.0 ], [ 1264401.513631681911647, -104188.64924998767674, 0.0 ], [ 1255148.083878455683589, -111527.525126685388386, 0.0 ], [ 1247811.151507917791605, -119020.223048191517591, 0.0 ], [ 1243559.735155994538218, -124702.111474551726133, 0.0 ], [ 1235582.033464746084064, -130480.617263145744801, 0.0 ] ] }
<osgeo.osr.CoordinateTransformation; proxy of <Swig Object of type 'OSRCoordinateTransformationShadow *' at 0x7f9f64652120> >

Stack trace:

 File "/opt/rh/httpd24/root/var/www/MyFile.py", line 380, in MyFunction
    intersectgeometry.Transform(self.coordTransform)
  File "/opt/rh/httpd24/root/var/www/wsgi-virtualenv/lib/python2.7/site-packages/GDAL-2.0.0-py2.7-linux-x86_64.egg/osgeo/ogr.py", line 5236, in Transform
    return _ogr.Geometry_Transform(self, *args)
RuntimeError: OGR Error: General Error

What am I doing wrong?


texts: I have a PostGIS road network layer (all features are linestrings) which I am accessing in Python and would like to get the geometry for each feature as a list of points.  Something like [(lat1,lon1),(lat2,lon2),(lat3,lon3)...] would be ideal.  I am currently accessing the geometry by using the ST_AsText(geom) function, and then parsing the resulting string into points.  But this seems clunky, and since this seems like a pretty common/basic operation I would expect there to be a standard way to do it, but I am not seeing anything in the documentation.  Is there a function or good way to do it out there that I am just missing?


texts: Select attributes in a column that are similar to attributes from another column

texts:
I have two map layers. One layer contains names of roads, on which some road accident occures. Those are the roads from whole region. The road names on this layer also contains the km or the number of a house, on which accident happend.

On the other layer I have names of roads only from the regional road network.

You can see that the recording method is very different. Now what I need is to select roads from the layer 1, that have the same name as on layer 2. I was thinking about using LIKE function, but I only found how to use it ones per attribute.


texts:
I am new to GIS but slowly getting the learning curve. I am working with QGIS version 2.18.7.
I have been editing several layers (both Raster and vector shape files) but suddenly the 'Toggle Editing' feature is greyed out and therefore I can no longer edit them. Although I have previously been able to do this. I have tried removing Filter but this still did not help.


texts: How to return coordinates of my location using LocationDisplayManager (ESRI API)?

texts:
I'm trying to create an application which works with lat and long of a device position. I'm using LocationDisplayManager class to find and to map a location.

MapView mMapView = null;

public void getPos(MapView mMapView) {
    LocationDisplayManager loc = mMapView.getLocationDisplayManager();
    loc.setAutoPanMode(LocationDisplayManager.AutoPanMode.LOCATION);
    loc.start();
}

public void showPos(View view) {
    getPos(mMapView);
}

@Override
protected void onCreate(Bundle savedInstanceState) {
    super.onCreate(savedInstanceState);
    setContentView(R.layout.activity_main);

    mMapView = (MapView) findViewById(R.id.map);

    if(getLon(mMapView) < 55.00 && getLat(MapView) < 55.00){
        TextView itemLoc = (TextView) findViewById(R.id.textLoc);
        itemLoc.setText("Lon: " + getLon(mMapView) + " Lat: " + getLon(mMapView));
    // code...
    }
}

I cannot return numerical values of lat and long that I need for future manipulations
I tried to implement the following code:

public double getLon(MapView mMapView) {
    LocationDisplayManager loc = mMapView.getLocationDisplayManager();
    return loc.getLocation().getLongitude();
}
public double getLat(MapView mMapView) {
    LocationDisplayManager loc = mMapView.getLocationDisplayManager();
    return loc.getLocation().getLatitude();
}

It gives error:

  Caused by: java.lang.NullPointerException at this line: return loc.getLocation().getLongitude();

Do you have any idea how to solve it?


texts: Images missing in QGIS print composer template using mapped drive file path

texts:
My question relates to this previous question.
QGIS Composer Big Red X
My images are missing when I load a print composer from a template
and all I see is red X.  Investigating this further the problem only occurs when I used a mapped network location e.g. Z:\QGIS\Composers, if I look at the saved file path it incorrectly reads ../../GIS/Composers.  If I use the full path \\myserver\GISshare\QGIS\Composers then there is no problem and image loads correctly.

All users have the same mapped drive that they commonly use to work with and are trying to access templates and images from this same location.  All users are using QGIS 2.18 and Windows 7 or 10.  Does anyone know why QGIS won't save the mapped file path correctly in this case and if it can be resolved?


texts: Callout lines option is missing on QGIS 3.4.3

texts:
I have seen that callouts are available starting with Qgis 3.1 but for some reason I don't have it in 3.4.

I have a work around but now I think there may be other options that aren't installed.


texts:
I noticed within the "Project Properties" there's an option for, Measurements --> Ellipsoid. What should this setting be configured to?
From what I read online, if this is set incorrectly, it can skew your measurements and provide inaccurate results.
Note: I'm currently using EPSG:26917. QGIS 3.14


texts: How to calculate the form (or shape) factor of a polygon?

texts:
I am using QGIS to calculate the form (or shape) factor of polygons following what Louf and Barthelemy did here: http://rsif.royalsocietypublishing.org/content/11/101/20140924

The form factor is defined as the ratio between the area of the block and the area of the circumscribed circle (something similar as a convex hull but with a circular buffer). The form factor is always smaller than 1 (since it is the area of the polygon divided by the area of the circumscribed circle).

To do so, I have a shapefile with the polygons of each block in several cities. I calculated the area of each polygon, but I need to find a way to calculate the area of the circumscribed circle. At the same time, it would be very useful if the calculation method can be repeated through a script or a function in the field calculator so it can be repeated for a large dataset of city blocks.


texts: Read specific columns shapefile with PySAL

texts:
I'm using PySAL for work and I have encounter a problem.
I want to calculate Moran’s I for serval fields, I can open and create a weight matrix from my shapefile (w) however I don't know how to select a specific field (y).

mir = pysal.Moran(y, w, permutations = 9999)

http://pysal.readthedocs.io/en/latest/users/tutorials/autocorrelation.html#moran-s-i

I just want something like :

f = "Path to my shapefile"
y = array from f

I don't have arcpy to read my shapefile, which I guess would make think easier.


texts:
I have installed up to 2.6.1 of QGIS on my computer.

Is it recommended that i remove all older versions or will it cause conflicts if they remain?


texts:
I have 2 spatial tables:
Streets: linestrings of streets from a country (more than 3 million)
Portals: points of adress from a country (around 18 million)
I need to move each point to the closest linestring matching 2 conditions, the query to do that is like this:
update portals pt
     set geom = (select st_closestpoint(s.streets_geom, pt.portal_geom)
                from streets s
                where pt.field1 = s.field2
                and pt.field3 = s.field3
                order by pt.portal_geom <-> s.streets_geom
                limit 1)

So far the query works well in tests but when I do it for all the data, it takes too long (several days) I have spatial index on both geometries.
How can I improve the query to do it faster? I can't figure it out
EDIT: Explain of the query


texts: PDF or HTML versions of Spatial Query Lab (FOSS4G Academy) lectures by Rick Smith

texts:
I am Spatial Science Teacher at TAFE NSW teaching Spatial Information Systems online.  In 2019 our department started including links to many of the Spatial Query Lab SQL lectures in our resources. We found them to be excellent as they were software-neutral and pitched lots of GIS concepts at just the right level for our students who need theoretical information at a simpler level than many university level texts.  The Labs for these classes are available here: https://github.com/FOSS4GAcademy/GST102FOSS4GLabs but not the lectures.
Does anyone know how to contact Rick Smith (previously at Richard.Smith@tamucc.edu) or anyone who might have a copy of these lectures?


texts: How to collect the graphics element drawn on a graphics layer?

texts:
The graphics element may consist of point, circles, polygon, free hand line drawing or marker symbols.

I m using arcobjects.


texts: Move the position of shapefile objects automatically

texts:
I have a shapefile with a big number of polygon objects. The problem is that these objects have a systematic shift.

Is there any possibility to place them all at the right position automatically using Python- or GDAL- scripts?

I'm using QGIS 2.10 on a Ubuntu 15.04 System.


texts: Add Column and Populate using the Python Console

texts:
I am trying to add a new column to the attribute table and populate the column with a string, using the python console in QGIS. I have done the same but populating the column with an integer (FIELD_TYPE=1), which works fine, but as soon as I change the column to a string (FIELD_TYPE=2), it no longer works. The code I am trying to run:

out3 = "%s\\roth_vh_%s_3.shp"%(shpPath,year)
processing.runalg('qgis:fieldcalculator',\
out2,\
"polar",\
2,\
4,\
0,\
True,\
'vh',\
out3)

The error message I get is "Evaluation error: Column 'vh' not found". So for some reason it is looking for a column vh rather than taking it as the string to populate all the fields in the column with. It's confusing me though, because if I open the attribute table and run a field calculation to populate a column with a string, I would simply select the column, type the string into the field calculation bar (e.g. 'vh') and click update all, and it would work. Any ideas why the field calculation of 'vh' is not working in pyqgis please?


texts: How to render values onClick() event in componentDidMount()

texts:
I am using OpenLayer library and I need change values in ComponentDidMount() but when I am clicking it does not affect.  I have used forceUpdate().

  handleDisable = () => {
  this.setState({ status: true });
  this.opacity_value = this.state.text;

  console.log(this.opacity_value);
  this.forceUpdate();
};

initMap = () => {

  console.log(this.opacity_value);
  wms_layer.setOpacity(this.opacity_value);

};
componentDidMount = () => {
  this.initMap();
};

render();
{
  return (
    <div>
      <Button
        content="Save"
        primary
        icon="save"
        onClick={() => this.handleDisable()}
      ></Button>
    </div>
  );
}


texts: Ordinary kriging problem after exporting to raster (Geostatistical Wizard)

texts:
I'm not sure what I'm making wrong. I used the Geostatistical Wizard to perform an ordinary kriging. The resulting Geostatistical Layer looks like this:

After using the GA Layer To Grid Tool by right clicking - Data - Export Raster, the resulting raster looks like this:

I think it's not a problem generated by the GA Layer to Grid Tool because if I go to Layer Prop - Symbology in the Layer Properties of the GA Layer, and change to Grid, the resulting grid also shows the breaks.

Can someone tell me, what is happening? Why do you see the regular borders in the second figure? Here you can see my input for the kriging:


texts: How to read Greek fonts (ISO-8859-7) in shapefile attributes within QGIS 1.8.0?

texts: I'd like to ask whether there is a way to read greek fonts in a shapefile's  attribute table within QGIS 1.8.0.


texts: Adding raster layer from WMTS URL using PyQGIS

texts:
I am trying to add a WMTS raster layer in QGIS 3.8 using Python (for my plugin).

I have found this Loading Layers - Raster in the documentation, but I cannot make it work.

I tried this:

QgsRasterLayer(f'url={url}&crs=EPSG:4326&format=image/png&tileMatrixSet=4326&username={name}&password={pwd}', 'some layer name', 'wmts').isValid()

But it returns False.

I see the same question asked a lot on this site, but no answer solves my problem.

How can I deal with this?

The URL works fine when I add the layer through the UI.

Other post : Loading a WMS layer to QGIS using python


texts: In a Geoserver Layer how do you make the color ramp smaller?

texts:
Is there a way to resize a color ramp decoration in a geoserver layer? As you can see in the WMS jpg request below the ramp overlaps my map.

I know how to anchor and move the color ramp, but ultimately I need to make it smaller.

Here is my layout.xml I see the size="auto" but when I try to change it the ramp simply doesn't display.

<layout>
  <decoration type="text" affinity="bottom,left" offset="36,6">
    <option name="message" value="2016 Leaf Out"/>
    <option name="font-size" value="22"/>
    <option name="font-family" value="Arial"/>
    <option name="halo-radius" value="2"/>
  </decoration>
  <decoration type="legend" affinity="bottom,right" offset="6,-55" size="auto"/>
</layout>


texts: On 4 of 74 Data Driven Pages dynamic text box prints out the codes instead of implementing them?

texts:
I am using Data Driven Pages in ArcGIS Desktop.

Here's the text box entry:

<UND>PWS Well Identification</UND>
<BOL>PWS Name:</BOL> <dyn type="page" property="attribute" field="OwnerID" domainlookup="true"/>
<BOL>Primary Station Code:</BOL> <dyn type="page" property="attribute" field="PS_CODE" domainlookup="true"/>
<BOL>Owner Well Name: </BOL> <dyn type="page" property="attribute" field="SOURCE_NAM" domainlookup="true"/>

On almost all of the DDP pages, the underlines and bolds are implemented properly, where the underlined bits are underlined and bold bits are bolded, but on 4 of them, it looks like this:

<UND>PWS Well Identification</UND>
<BOL>PWS Name:</BOL> Turlock Golf & Country Club
<BOL>Primary Station Code:</BOL> 2400078-001
<BOL>Owner Well Name: </BOL> WELL #1 - CLUBHOUSE/POOL AREA

Its printing out the codes, but its going and getting the information from the attribute table like its supposed to.

That's on 4 of the 74 pages. the rest work normal.

Any ideas why this is happening?