texts:
I am trying tu use David Rumsey' Georeferencer tiles  in R leaflet package without success (Only the OSM basemap is displayed).

According to the package vignette I tried :

   library(leaflet)
            # XYZ direct link
            leaflet() %>%
              setView(0.65, 0, zoom = 5) %>%
              addTiles() %>%
              addTiles("https://maps.georeferencer.com/georeferences/700526190853/2017-12-30T11:48:27.589686Z/map.json?key=D7AwmpRP1H6pUic6DIK3")

and

 library(leaflet)
            # WMS tiles
    leaflet() %>%
      addTiles() %>% setView(0.65, 0, zoom = 5) %>%
      addWMSTiles(
        "https://maps.georeferencer.com/georeferences/700526190853/2017-12-30T11:48:27.589686Z/wmts?key=D7AwmpRP1H6pUic6DIK3&SERVICE=WMTS&REQUEST=GetCapabilities",
    layers = "1", # I assume
    options = WMSTileOptions(format = "image/png", transparent = FALSE),
    attribution = "")

Note that the access at the map links need registration.
Note that I used the 2.1 package release.


texts: Data defined symbol rotation for a categorized renderer

texts:
I have a CSV File with coordinates and a certain column called azimuth. Azimuth field can have any value from 0 to 360 degrees.  I want to plot the points as triangles, but rotate them depending upon the azimuth column.

I understand that , there is a option to categorize the style. However, I have to go to each category and  change the rotation manually.

How can I make the rotation value read directly from the Azimuth


texts: Grouping a dataframe whilst making a pie chart

texts:
I have a data layer which consists of an attribute table of ~30 columns and ~6000 rows. This table contains information about buildings that are located close to parks in Amsterdam. I want to create a map in which I create for each park in Amsterdam a pie chart which shows the amount of buildings near that park. I know how to make the pie charts and how to display these on a map, the issue here is that I need to group the data to fall into parks. As of now, when I make the pie charts, the pie charts are made for EACH building (as I have 6000 buildings) and hence I need to group these buildings to each park. I have around 6000 observations for just 20 parks. So ideally I want to have 20 pie charts, one for each park.
I used a few tools (e.g. summary statistics) but these all give me back a table and hence I cannot make pie charts out of this anymore. I have attached three screenshots, 1 of the table, 2 of the output I got (and which I want to group to parks). I have also added a example of how I want my data to look like (note: the counties here substitute the parks in my data)


texts: Changing format of column from character to decimal or float in Mapinfo

texts:
We have some Mapbasic software that updates character fields in existing mapinfo tables to decimal or float. The character fields are mostly numeric but some values are N/A to indicate not applicable which is more correct for us that 0.

This is working fine for values under 1,000 (note that the format has a comma in it) but for values greater than 1,000 it is treating the comma as a decimal point I think. Any numbers greater than this are altered to 1 2 etc
WE are using commands as follows:

    Alter Table "tmp_astext" ( modify LV Float )
    Alter Table "tmp_astext" ( modify LV Decimal(16,4) )

Is there any way we can stop this behaviour without having to make the users remember to change the formats when creating the files initially (they are created from excel). I've tried using the command

        Set Format Number "9,999.9"

in the Mapbasic window but it doesn't seem to help


texts: How to stop QGIS auto snapping new poly lines as they are being created?

texts: Im creating a polyline based data set which involves digistising around a large amount of street outlines. Is there a way to stop QGIS from auto snapping the new polyline section to a straight ahead bearing? See attached image for what I mean. The red line being the existing line section.


texts:
I have a shapefile which I read in via
import geopandas as gpd
shapes_df = gpd.read_file(shp_f)
print(shapes_df.crs)

The output was {'init': 'epsg:3857'}
I then reprojected via
shapes_df = shapes_df.to_crs("epsg:32629")
print(shapes_df.crs)

The output was epsg:32629
What is the meaning of the 'init' value for the CRS in the original CRS? Should I be concerned that there is no 'init' value in my reprojected data? Did the projection not work properly?


texts: Displaying rasters in a folder | QGIS 3.0.3

texts:
I have a satellite image. And in the folder where this picture is stored, it is well displayed.

But when I divide it with the help of script into parts and save in a separate folder, then get such a display of parts in the folder.

Here is the script:

# grid - grid lying on the satellite image

for i in range(grid.featureCount()):
        grid.select(i)
        OUTPUT = 'LC08_L1TP_155016_20180203_20180203_01_RT_B1_{}'.format(i+1)
        parameters = {'INPUT': 'D:/Study/PROJECT/Storage scene/LC08_L1TP_155016_20180203_20180203_01_RT_B1/LC08_L1TP_155016_20180203_20180203_01_RT_B1.TIF',
                        'MASK': QgsProcessingFeatureSourceDefinition(grid.name(), True),
                        'NODATA' : 0,
                        'ALPHA_BAND' : False,
                        'CROP_TO_CUTLINE' : True,
                        'KEEP_RESOLUTION' : False,
                        'OPTIONS' : '',
                        'DATA_TYPE' : 5,
                        'OUTPUT' : 'D:/Study/PROJECT/Storage clip/LC08_L1TP_155016_20180203_20180203_01_RT/LC08_L1TP_155016_20180203_20180203_01_RT_B1/' + OUTPUT + '.TIF'}
        processing.run('gdal:cliprasterbymasklayer', parameters)
        grid.removeSelection()

Maybe it depends on the parameters?

How to fix it?

In QGIS this is shown normally.

But in the folder it looks like this:

How to make a part in the folder displayed normally, as in QGIS?

In the future, OpenCV will use these parts of the satellite image to determine the images. ("В будущем OpenCV будет использовать эти части изображения спутника для определения образов." in Russian)


texts:
I'm working on a stormtracking project, and I realized that the longitude sign of some values in my GeoPandas dataframe are wrong; i.e. I need negative longitudes instead of positive ones.
When I trace my storm tracks, I have an ugly horizontal line (the red one) which can be explained by a change in the sign of longitude in geometry.linestring.
For example:

My GeoDataFrame looks like:
    DEBUT    MEMBER FIN        ID   VORT    PRES    CIRC           geometry
11  2020112415  0   2020112418  49  78.0    952.88  0.0 LINESTRING (-175 60, -174.75 61.5)
12  2020112418  0   2020112421  49  80.0    950.16  0.0 LINESTRING (-174.75 61.5, -175.25 62.5)
13  2020112421  0   2020112500  49  81.0    948.5   0.0 LINESTRING (-175.25 62.5, -176.5 63.5)
14  2020112500  0   2020112503  49  81.0    946.03  0.0 LINESTRING (-176.5 63.5, -177.75 64)
15  2020112503  0   2020112506  49  79.0    945.21  0.0 LINESTRING (-177.75 64, -179.5 64.25)
16  2020112506  0   2020112509  49  77.0    945.48  0.0 LINESTRING (-179.5 64.25, 178.75 64.5)

Then I'm trying to change the last row with spent many hours on it:
16  2020112506  0   2020112509  49  77.0    945.48  0.0 LINESTRING (-179.5 64.25, 178.75 64.5)

Is there any function in GeoPandas to easily manipulate the coordinates?
I'm turning around ... here's my last attempt to modify longitudes:
def swap2(x):
    coords = list(x.coords)
    coords = [Point(-t[0], t[1]) for t in coords]
    return LineString(coords)

tempete=tempete.geometry.map(swap2)


texts: pgRouting problem: creating topology where segments touch

texts:
I am working my way through the example in the pgRouting manual (2.2). After running pgr_createTopology() and pgr_analyzeGraph():

pgr_test=# SELECT pgr_createTopology('edge_table',0.001);
NOTICE:  PROCESSING:
NOTICE:  pgr_createTopology('edge_table', 0.001, 'the_geom', 'id', 'source', 'target', rows_where := 'true', clean := f)
NOTICE:  Performing checks, please wait .....
NOTICE:  Creating Topology, Please wait...
NOTICE:  -------------> TOPOLOGY CREATED FOR  18 edges
NOTICE:  Rows with NULL geometry or NULL id: 0
NOTICE:  Vertices table for table public.edge_table is: public.edge_table_vertices_pgr
NOTICE:  ----------------------------------------------
 pgr_createtopology
--------------------
 OK
(1 row)

pgr_test=# select pgr_analyzeGraph('edge_table', 0.001);
NOTICE:  PROCESSING:
NOTICE:  pgr_analyzeGraph('edge_table',0.001,'the_geom','id','source','target','true')
NOTICE:  Performing checks, please wait ...
NOTICE:  Analyzing for dead ends. Please wait...
NOTICE:  Analyzing for gaps. Please wait...
NOTICE:  Analyzing for isolated edges. Please wait...
NOTICE:  Analyzing for ring geometries. Please wait...
NOTICE:  Analyzing for intersections. Please wait...
NOTICE:              ANALYSIS RESULTS FOR SELECTED EDGES:
NOTICE:                    Isolated segments: 2
NOTICE:                            Dead ends: 7
NOTICE:  Potential gaps found near dead ends: 1
NOTICE:               Intersections detected: 1
NOTICE:                      Ring geometries: 0
 pgr_analyzegraph
------------------
 OK
(1 row)

So there are 2 isolated segments, which turn out to have id = 17, 18:

pgr_test=# SELECT a.* FROM edge_table a, edge_table_vertices_pgr b, edge_table_vertices_pgr c WHERE a.source=b.id AND b.cnt=1 AND a.target=c.id AND c.cnt=1;

 id | dir | source | target | cost | reverse_cost | x1  | y1  |       x2       | y2  |                                      the_geom
----+-----+--------+--------+------+--------------+-----+-----+----------------+-----+------------------------------------------------------------------------------------

 17 | B   |     14 |     15 |    1 |            1 | 0.5 | 3.5 | 1.999999999999 | 3.5 | 010200000002000000000000000000E03F0000000000000C4068EEFFFFFFFFFF3F0000000000000C40

 18 | B   |     16 |     17 |    1 |            1 | 3.5 | 2.3 |            3.5 |   4 | 0102000000020000000000000000000C4066666666666602400000000000000C400000000000001040

The id=18 case is due to crossing and can be fixed later by pgr_nodeNetwork(), but the id=17 case is due to touching and I thought it should be handled by createTopology(), and it persisted even after pgr_nodeNetwork(). I tried bigger tolerance values to no avail. The manual says tolerance is in projection unit, but I don't know whether this is relevant for the sample data in the manual?

I am using pgRouting-2.2.0 on Windows with Postgresql-9.5, PostGIS-2.2.


texts: Performing double classification for symbols in QGIS?

texts:
I am using QGIS version 2.14.1 (Essen). I have symbols that I have classified according to certain value from one of the columns in my shapefile.

Is it possible to do a second classification according to a second column?

Basically, I am using data with road works and the deviations linked to roadworks. My first classification is according to the size of the roadwork and I would like the second one to be according to the type of deviation (for cars, pedestrians, buses, etc).


texts:
Really don't know how to ask this right, but here it goes. How to add condition to ST_ClusterWithin to cluster data depending on distance AND some other attribute, like type of something?

  ST_ClusterWithin — Aggregate. Returns an array of
  GeometryCollections, where each GeometryCollection represents a set of
  geometries separated by no more than the specified distance.

Sample code from documentation:

WITH testdata AS
  (SELECT unnest(ARRAY['LINESTRING (0 0, 1 1)'::geometry,
               'LINESTRING (5 5, 4 4)'::geometry,
               'LINESTRING (6 6, 7 7)'::geometry,
               'LINESTRING (0 0, -1 -1)'::geometry,
               'POLYGON ((0 0, 4 0, 4 4, 0 4, 0 0))'::geometry]) AS geom)

SELECT ST_AsText(unnest(ST_ClusterWithin(geom, 1.4))) FROM testdata;

Result:

st_astext
---------
GEOMETRYCOLLECTION(LINESTRING(0 0,1 1),LINESTRING(5 5,4 4),LINESTRING(0 0,-1 -1),POLYGON((0 0,4 0,4 4,0 4,0 0)))
GEOMETRYCOLLECTION(LINESTRING(6 6,7 7))


texts: How to disable double-click zoom in Arcgis Javascript API

texts:
I am using ArcGIS java script API for creating map. I am trying to disable the double-click zoom in the map. Here is link of my map: http://urbanaccessregulations.eu/mindcrew/mobile-popup.html

I used both methods which are given in documentation of Arcgis JS:

map.isDoubleClickZoom=false;
map.disableDoubleClickZoom();

But still the double click zoom is not disabled - what am I doing wrong?


texts: How to conduct logistic regression between two rasters

texts:
I have two rasters. The first raster is a region showing sinkholes(depressions) of the region(cell values are 1 for sinkholes and 0 for non-sinkhole cells). The second raster is the TPI(Topographic Positioning System) of the same region. I want to conduct a logistic regression analysis between the two rasters(Sinkhole~TPI). can anybody kindly tell me how I can do this?(preferably R)

(The purpose is to build a sinkhole detection model using TPI index so that it can use the model to detect sinkholes in other parts of the study area. )


texts: Selecting features in an ArcGIS attribute table based on values of other features in the same table using SQL

texts:
I have a line feature class (saved in a feature dataset within a File GDB) which is comprised of a set of intersecting lines.

I have used the Add Geometry tool to add Start_X, Start_Y, End_X and End_Y fields to the attribute table.  The lines are not all orientated in the same direction, so for certain pairs of features in the table, End_X = End_X and End_Y = End_Y.
I want to use the Select Layer By Attributes tool with an SQL query (if possible) to select only lines in the table with matching end coordinates and then run the Flip Line tool on the selection in order to reverse the direction of the selected lines.

It's also worth noting that the OIDs are in a random order, e.g. line 1 may be connected to line 10.

Is it possible to use an SQL statement to select features in an attribute table based on the values of other features?

I tried

SELECT * FROM table_name WHERE:
END_X = END_X AND END_Y = END_Y AND OBJECTID <> OBJECTID

I'm sure it will be obvious to most that all that does is search each row individually and come up with nothing, though hopefully it gives an idea as to what I'm trying to do.

If it is not possible using SQL, is it something that could be done with Arcpy Search Cursors?  (I am an extreme novice with Arcpy/Python scripting. I just know Search Cursors exist and can be used to find records in tables.)


texts: Hiding password of PostGIS layer in QGIS Plugin

texts:
How can I hide the password in QGIS Vector Layer (PostGIS) in a Plugin (via PyQGIS)?
If the layer is created.
uri = QgsDataSourceURI()

uri.setConnection(self.uri.host(), self.uri.port(), self.uri.database(), self.uri.username(), self.uri.password())

uri.setDataSource(schema,vista, geom, '', pk_field)

vlayer = QgsVectorLayer(uri.uri(False), vista,'postgres')

I just opened the Python console in QGIS, select the layer to discover the database password and type.
layer = iface.activeLayer()
layer.source()

And it displays the user and password of the PostGIS database I am connecting.
I want to avoid that so important information, like the user and password of the PostGIS database I am connecting to, could be displayed so simply, without a way to hide it.


texts: Features show on data view, but not in Attribute table?

texts:
There are 966 points for tree data that was collected in the field. When I try to open the attribute table only 2 points' data show up. I've tried to create layer from selected, because when I choose select all, it says:

  966 out of 2 selected

But that hasn't worked, and I have also tried to export the data to Excel as well says transfer the file to a new map completely, nothing has worked. When I try to sort the points an error shows up. Is the layer just completely corrupt?

The original data was downloaded from ArcGIS Online as a file geodatabase. The data was collected with Trimble Catalyst receivers into a hosted feature layer on ArcGIS Online. All we are trying to do at this point is get the attribute table to copy over to Excel.


texts:
This seems like it ought to be straightforward:

from osgeo import gdal
gpkg = gdal.OpenEx('path_to_gpk',gdal.OF_VERBOSE_ERROR,allowed_drivers=['GPKG'])
error = gpkg.DeleteLayer(0)

After running this error is always a OGRERR_FAILURE. The geopackage has been created by QGIS so I would have thought should be standard.

Can't see any way to debug further into the GDAL calls apart from a GDB session, which I might try later.


texts: Access the array of a loaded layer in QGIS3 with pyQGIS

texts:
The code below is a clunky attempt to get a raster numpy array loaded up for basic machine learning from the QGIS Python console.  It usually works but has stability issues with the GetRasterBand() and I have asked a separate question.  As an improvement I am trying to figure how to get the array values directly from the loaded rasters in QGIS.  I.e. you open the raster layers in QGIS, do whatever is needed (e.g. training polygons) and then run the code.

#Part 1. Importing Python modules
import numpy as np
import gdal
from sklearn.cluster import KMeans
import os

#Part 2. Specify the image file to be classified and the class raster with the training data.

MyFolder = 'D:\\RemoteSensing\\Practical_7\\' #Change this to match the folder on your drive.  Note that we use \\ instead of \
MyImageRoot = 'S2_mini' #Specify the root image name. Do not nclude the extension, e.g. 'glacier_'.  If a color image, just write the file name without extension.
Nclasses = 4 #how many classes are in the image
MyClassRaster = 'KmeansTEST_4cls'
BandList = ['B02'] #Indicate the bands you want to load.  B

#Part 3 Load the image
fileName = os.path.join(MyFolder + MyImageRoot + BandList[0] + '.tif')
ImageFile = gdal.Open(fileName)
band1 = ImageFile.GetRasterBand(1)
array1 = band1.ReadAsArray()

I am trying to get array1 by direct querry to a loaded band in the QGIS layers.  I have found QgisProject:

layer = QgsProject.instance().mapLayersByName(BandList[0])

But I can't find how to access the raster values for the layer and have them as a numpy array for the next processing steps.  I'm working with QGIS 3.2 and 3.4.


texts: Python Window Crash and import arcpy issue

texts:
Trying to do anything in the Python shell within ArcGIS 10.2 and it just crashes once I hit enter after writing code, no matter what I do. I've tried to uninstall, deauthorize, reinstall, re-authorize, hoping it would repair, but it still crashes.

On a similar note (not sure if they're connected), I can not import arcpy outside of ArcGIS, whether it's in the 2.7.3 shell or WING IDE. It prompts that:

  raise TypeError("This object does not support enumeration")

  TypeError: This object does not support enumeration

I don't understand any of this, I've always had an issue trying to import arcpy, I don't get it. I thought it was as simple as updating the Python paths to the arcpy library. Any help or suggestions?

Here's the result of sys.path:

  ['', 'c:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy', 'C:\',
  'c:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy\arcpy\gdal',
  'C:\WINDOWS\SYSTEM32\python27.zip',
  'c:\Python27\ArcGIS10.2\DLLs', 'c:\Python27\ArcGIS10.2\lib',
  'c:\Python27\ArcGIS10.2\lib\plat-win',
  'c:\Python27\ArcGIS10.2\lib\lib-tk', 'c:\Python27\ArcGIS10.2',
  'c:\Python27\ArcGIS10.2\lib\site-packages', 'C:\Program Files
  (x86)\ArcGIS\Desktop10.2\bin', 'C:\Program Files
  (x86)\ArcGIS\Desktop10.2\ArcToolbox\Scripts',
  'c:\Python27\ArcGIS10.2\lib\site-packages\win32',
  'c:\Python27\ArcGIS10.2\lib\site-packages\win32\lib',
  'c:\Python27\ArcGIS10.2\lib\site-packages\Pythonwin']


texts: Downloading OpenStreetMap data by osmosis?

texts:
I need to download a geographic data and unzip the download file using osmosis command. I tried this:

osmosis --rx test.osm --bb top=48.78 left=6.02 bottom=48.59 right=6.38

but this command filters OSM files by BBox and it does not import the data from the OSM server.


texts: Pulling images into QTWidget (plugin dialog box) from web server?

texts:
I'm trying to show images in a dialogue box (built using QT Creator and interacted with through my plugin).  I have the images as a list of URLs, but can't figure out how to get them to display in the dialogue box. I thought perhaps I could just include an HTML image tag, but that just shows a file icon.

I wonder if I need to download them into a temp filestore first.


texts: Is there a way to automate the preprocessing of LANDSAT data?

texts:
I'm currently enrolled in a remotesensing course, using LANDSAT data.
Part of my course is about standardizing the scenes based on radiance values. I assume that because this is such a tedious task, someone has already figured out a way to automate it.

Are there any tools, ESRI or otherwise that will standardize LANDSAT scenes?


texts: Accessing File Geodatabase to convert to shapefiles using ArcObjects?

texts:
The end result of what I'm trying to do is convert a FGDB into a folder of shapefiles. To do this, I need access to all the data I'd expect to see in a shapefile. Once I get all the information I need from the file geodatabase, it should be really easy.

I am stuck on trying to read the file geodatabase spatial data (Xs and Ys)

UPDATE: Using the snippet of code posted by Rich Wawrzonek below, I can read the names of the datasets. Now I need to know how to read the Xs and Ys of each dataset!

    IEnumDataset enumDataset = iwrk.Datasets[esriDatasetType.esriDTFeatureClass];
    IDataset dataset;
    while ((dataset = enumDataset.Next()) != null)
    {
         //Do something with the feature class here
         Console.WriteLine(dataset.Name);
    }


texts: How to do a Correlation of various Raster Data and show the result in a Correlation Matrix

texts:
I am doing spatial correlation of soil Data sets.
I did various Kiriging Interpolations with the krige function. The output was cut to a target raster to have the same raster size for all Kriging Interpolations of my datasets.
Usually to do a correlation between two kriging outputs I use the cor.test function method = "pearson". Later I copy and paste the Consule result to Excel and calculate the r-squared and creat a Correlation Matrix.
But now I have so much Kriging output that it would take forever to do everything manually!
Is there a possibility to do a correlation of all Kriging Output (Formal Class SpatialGridDataFrame) with one function and demonstrate the result in a Correlation Matrix?


texts: Exporting DDP to multiple multi-page pdf files

texts:
In ArcMap 10.4 I have an index layer in DDP enabled mxd. The index layer has 50 unique Name values (used as the PageName in DDP) but there are more rows with the same Name value providing different scales, spatial extents etc. I'm looking for an ArcPy script to export pdf files with all pages of the same Name.

Example:

Name   Scale
A      5000
A      600
B      25000
B      5000
B      2500
...

Result:

A.pdf (2 pages)
B.pdf (3 pages)
...

Would you be so kind to provide some hint? Apparently, I'm stuck in the process of correctly defining the 'for' loop.

So far I tried to use the modified code from Cindy Jayakumar https://gis.stackexchange.com/a/71702/93026, which was closest to what I found:

import arcpy
import os

mxd = arcpy.mapping.MapDocument("CURRENT")
output_folder = r"C:\temp"
ddp_lyr = mxd.dataDrivenPages.indexLayer
pdf_field = "Name"

#Get a list of unique values to group by
pdf_list = set([row[0] for row in arcpy.da.SearchCursor(ddp_lyr,pdf_field)])

for pdfName in pdf_list:
    arcpy.mapping.PDFDocumentCreate(os.path.join(output_folder,pdfName + "_inv.pdf"))
    for i in range(1, mxd.dataDrivenPages.pageCount + 1):
        mxd.dataDrivenPages.currentPageID = i
        cur = str(mxd.dataDrivenPages.pageRow.getValue(pdf_field))
        if cur == pdfName:
            tempPdf = os.path.join(output_folder,"temp_" + cur + str(i) + ".pdf")
            arcpy.mapping.ExportToPDF(mxd, tempPdf)
            finalPdf.appendPages(tempPdf)
            os.remove(tempPdf)
    finalPdf.saveAndClose()
    del finalPdf, tempPdf

The above code does what I need to some extent as from some reason does not export the pdf files for all unique values in the "Name" field. I was not able to identify why some have been exported and some not. I know there is much easier way to do this as 2 years ago I already had the code working (though not saving it).


texts: QGIS only reads first polygon in KML file

texts:
I am new to QGIS and I am struggling to read in a KML file of polygons.

I have a KML file that contains polygon information for a number of squares formatted as follows:

<?xml version="1.0" encoding="UTF-8"?>
    <kml xmlns="http://www.opengis.net/kml/2.2">
<Placemark>
<name>Place@001@001@010</name>
<Polygon><outerBoundaryIs><LinearRing>
<coordinates>
25.271,-13.416
25.271,-13.415
25.272,-13.415
25.272,-13.416
25.271,-13.416
</coordinates>
</LinearRing></outerBoundaryIs></Polygon>
</Placemark>

<Placemark>
<name>Place@001@001@011</name>
<Polygon><outerBoundaryIs><LinearRing>
<coordinates>
25.371,-13.426
25.371,-13.425
25.372,-13.425
25.372,-13.426
25.371,-13.426
</coordinates>
</LinearRing></outerBoundaryIs></Polygon>
</Placemark>

<Placemark>
<name>Place@001@001@011</name>
<Polygon><outerBoundaryIs><LinearRing>
<coordinates>
25.391,-13.456
25.391,-13.455
25.392,-13.455
25.392,-13.456
25.391,-13.456
</coordinates>
</LinearRing></outerBoundaryIs></Polygon>
</Placemark>

</kml>

When I do: Layer -> Add Layer -> Add Vector Layer and choose my file, it does load into QGIS but only a single square is plotted. I imagine I am making an elementary mistake (likely in the format of the KML file),  but I cannot find a reference to figure it out.

QGIS Version: 3.8 Zanzibar


texts: How to create a dropdown menu in QGIS Toolbar with Python

texts:
I Have a set of some QAction in my QGIS-Plugin and they are created and added to the toolbar like this

icon = QIcon(":/action1.png")
self.action1 = QAction(icon, u"Action 1", self.iface.mainWindow())
QObject.connect(self.action1, SIGNAL("triggered()"), self.someMethod)
self.iface.addToolBarIcon(self.action1)
self.iface.addPluginToMenu(u"&My Plugin", self.action1)

But now I'd like to group my Actions like in this "Select Feature" Screenshot:
http://img195.imageshack.us/img195/6375/selectfeaturedropdown.png

Is this possible with Python in QGIS?


texts: Creating noise maps using ArcGIS for Desktop or Open Source Desktop GIS?

texts:
In a previous role I helped create Noise Maps for Northern England.

The predicted noise levels was done using the IMMI package (which I can recommend) and the cartographic output was done in ArcMap.  There are also other packages such as LimaArc which have tighter integration with ArcMap.

I want to take on a personal project to try and create noise maps for large urban areas in Australia.  I am interested in both mapping noise from traffic (roads, trams, rail), maybe having a stab at aircraft noise and also looking at how the noise from level crossing bells can be heard by surrounding residents.  As its for self development purposes, I have no budget for buying dedicated software like LimaArc or IMMI.

I was wondering if anyone in the GIS-SE community has found a way to create noise maps in a similar approach to how IMMI/LimaArc does?

--EDIT--
I had a play with the noise module for OrbisGIS (user16714 answer and comments below), but had some problems getting it to work - looks like it has a lot of potential but still has to be installed against a beta version and not entirely sure how well it will handle larger datasets. (Will post back if I get further with this).

Ideally I really want to do this all in ArcGIS or a more established open source package (QGIS, GRASS, SAGA, etc.).  Have seen some papers that mention tools to do this in ArcGIS, but can't source any scripts/tools as of yet.

Has anyone seen any scripts (Python) or free extensions where this can be achieved directly in ArcGIS?


texts:
How can I apply rasterizing process in QGIS 3.0. I know it is very simple quesiton. However I cannot find the rasterize tool.

EDIT

I rasterized a layer with GDAL\rasterize tool by using script which is:

parameters={"INPUT": layer,
                           "FIELD": "PolygonId",
                           "DIMENSIONS": 0,
                           "WIDTH": rps,
                           "HEIGHT": rps,
                           "RAST_EXT": "%f,%f,%f,%f" % (xmin, xmax, ymin, ymax),
                           "TFW": 1,
                           "RTYPE": 6,
                           "NO_DATA": 0,
                           "COMPRESS": 0,
                           "JPEGCOMPRESSION": 1,
                           "ZLEVEL": 1,
                           "PREDICTOR": 1,
                           "TILED": False,
                           "BIGTIFF": 2,
                           "EXTRA": '',
                           "OUTPUT": dosya}
processing.runAndLoadResults("gdal:rasterize",parameters)

However it does not work now.
I found the rasterize.py under this:

C:\OSGeo4W64\apps\qgis-dev\python\plugins\processing\algs\gdal\rasterize.py

I cannot use it again. What is the solution for rasterize by using at least these parameters;

Input, Field, Width, Height, Rasterize extension, Output


texts: Joining separate parcels based on group ID into one large group

texts:
What tool would I use to group many parcels into its own respective group. For example, if parcels A, B , C have group ID of 1 then I would like to have parcels A, B, and C combined into one larger parcel. Is there a tool that does this without me manually digitizing it?

Also, I have another layer called street center lines, and would like the larger parcel to have the street center line as its boundary if it is lying next to a street center line.

Please refer to the images below

https://imgur.com/a/FpA2fQ7

The first photo shows many parcels that have the same ID group.

The second photo shows a close up view of street center lines next to the parcels.

The third photo shows the desired result. Each larger grouped parcel has a boundary that is next to each other.


texts: Merging two segments with overlapping nodes in QGIS?

texts: I have two line segments that overlap nodes. I want to select them, and merge them into one complete segment, so that the overlapping nodes are combined to one node.


texts: Zonal statistics of polygon overlays in PostGIS

texts:
I am trying to do a seemingly simple task -- extract mean pixel values from a raster based on a polygon overlay in PostGIS. I consulted several sources, including this blog, this blog, and this post). I loaded a Landsat tile into my database (in-db) and a polygon overlay consisting of three polygons (ID 1,2,3). The Landsat tile is a pansharpened image with a 15m spatial resolution and 4 bands.

Based on help files and the aforementioned posts/questions, I arrived at the following PostGIS query:

SELECT id, (SUM((ST_SummaryStats(a.rast, 1, true)).sum)/SUM((ST_SummaryStats(a.rast, 1, true)).count)) as mean
FROM imagery.l8_2015_09_12 AS a, analysis_results.zonal_stats_test as b
WHERE ST_Intersects(b.geom,a.rast)
GROUP BY id;

The results are:

1: 9902.49
2: 10079.68
3: 12355.90

However, when I test the results against ArcGIS's "zonal statistics to table" tool, I get the following results:

1: 8089.61
2: 7527.62
3: 12290.05

Close, but not identical. I am almost there but I was wondering if someone could help me figure out the source of the difference between both techniques. I suspect that there is something wrong with my ST_summaryStats method.


texts: How to add source in GeoJSON layer using OpenLayers

texts:
I am trying to add GeoJSON layer using OpenLayers which is in the PostgreSQL db. When I am running query as

select ST_AsGeoJson(ST_Transform(geom,900913)) as farmcode from farmdatadd where areatbl = 26

I am getting multiple row of MultiPolygon on which I want to add as a layer on map. But only the first array of layer is added and remaining all is not appearing.
My code is:
function addGEOJSONLayer(geojson) {
            var vectorSource = new ol.source.Vector({
                  features: new ol.format.GeoJSON().readFeatures(JSON.parse(geojson[0])),
             });
            if(geojson.length > 1) {
            console.log("here");
                for(var k = 1; k < geojson.length; k++) {
                console.log(k);
                    vectorSource.addFeature(new ol.Feature(new ol.format.GeoJSON().readFeatures(JSON.parse(geojson[k]))));
                }
            }

            var vectorLayer = new ol.layer.Vector({
                  source: vectorSource,
                  style: styleFeature
            });
            var ext = vectorLayer.getSource().getExtent();
            map.getView().fit(ext, map.getSize());
            map.addLayer(vectorLayer);
        }

I have done the same thing in old version of Openlayer as :
if (jsonString != "") {
        //var holdingFeature = JSON.parse(jsonString);
        var holdingFeature = jsonString;
        console.log(holdingFeature);
        console.log(holdingFeature.length);
        for (var i = 0; i < holdingFeature.length; i++) {
            console.log(holdingFeature[i]);
            vector_layer_feature.addFeatures(geojson_format.read(holdingFeature[i]));
        }
    }
    for (var j = 0; j < vector_layer_feature.features.length; j++) {
        vector_layer_feature.features[j].style = { fillColor: "#4169E1", fillOpacity: 0.5 };
    }
    map.zoomToExtent(vector_layer_feature.getDataExtent());
    vector_layer_feature.refresh({ force: true });

Any suggestions?


texts: Can I add a spreadsheet with addresses to QGIS and then geocode it?


texts: Polygon layer with only one feature and lots of separated areas

texts: using qgis, I have a polygon layer. It has only one feature in the attribue table. But it has lots of different areas, clearly separated from each other. Now I want to have this separate areas also in separate features, to be able to sort according to size etc. Trying different ways with the Geoprocessing tools did not help (singlepart to multipart, dissolve).


texts:
I am geocoding some addresses and it works just fine. However, all the points are rearranged afterwords. How can I stop it from rearranging the points? Is this a locator issue?

Geoprocessor GP = new Geoprocessor();
            GeocodeAddresses geocodeAddresses = new GeocodeAddresses();
            geocodeAddresses.in_table = workSpaceLocation + "\\" + addressTable;
            geocodeAddresses.address_locator = addressLocator;
            geocodeAddresses.in_address_fields = "Street Address;" +
               "City <None>;" + "State <None>;" +
               "ZIP <None>";
            geocodeAddresses.out_feature_class = workSpaceLocation + "\\" + outputFeatureClass;
            geocodeAddresses.out_relationship_type = "STATIC";

                IGeoProcessorResult result = GP.Execute(geocodeAddresses, null) as
                    IGeoProcessorResult;
                if (result != null)
                {
                    if (result.Status != esriJobStatus.esriJobSucceeded)
                        MessageBox.Show("Failed to Geocode.", "Error", MessageBoxButtons.OK, MessageBoxIcon.Error);

                }
                else
                {
                    MessageBox.Show("Failed to Geocode.", "Error", MessageBoxButtons.OK, MessageBoxIcon.Error);
                }


texts: Backwards compatibility for QGIS 3.2.1 Processing Modeler

texts:
I have created multiple models using the Processing Modeler in QGIS 2.18.15. However, I could not open these .model files in the user interface for the Graphical/Processing Modeler in QGIS 3.2.1, which seems to only accept .model3 or .MODEL3 files. Also, there is an option to save the (.model) file as a Python (.py) file in QGIS 2.18.15.

Is there a way that I can make either the (.model) or (.py) files from QGIS 2.18.15 compatible for QGIS 3.2.1 Graphical/Processing modeler? Since the alternative to that would be to re-create the Processing Modeler manually from scratch.


texts: Python code for locking elements in Layout view

texts: Is it possible to lock elements in layout view in Python? I know it is possible in the production cartography extension from ESRI but all I need is the funcionality to lock elements. I am wondering if there is a way to do this in Python.


texts: Performing rigid translation of polygon using GeoTools Java API?

texts:
I've been spending some days trying to find the best way to do a "rigid translate" of a polygon in a certain direction and distance, without much luck. To provide an example, consider translate like moving a polygon on 30 degrees south, for 500 meters.

Is there any API in Java GeoTools that you have used successfully to do such moving of a polygon?

If not for the polygon as a whole in one operation, I can do a single point translation for each point of the polygon, if no better solution can be found. I don't need to do this real time, so I have processing time to spend to do this.

The approach I'm following so far is of calculating the shifting on the single point as for this link to the geodetic calculator

In the specific the part called:

Generate location away from a point

    `GeodeticCalculator calc = new GeodeticCalculator();
    // mind, this is lon/lat
    calc.setStartingGeographicPoint(45.4644, 9.1908);
    calc.setDirection(90 /* azimuth */, 200 /* distance */);
    Point2D dest = calc.getDestinationGeographicPoint();
    System.out.println("Longitude: " + dest.getX() + " Latitude: " + dest.getY());`


texts: Trying to convert NetCDF file to shapefile in QGIS

texts:
I have a NetCDF file containing surface sea temperatures that I want to convert into a ESRI shapefile.
I have tried following this link: https://freegistutorial.com/how-to-export-layer-to-shapefile-on-qgis/
but I cannot see the ESRI Shapefile format in the drop down list.

I loaded the CDF file -> right click -> Export -> save as, and I come to this screen where I can choose the format type.

There is no shapefile option and I am unsure what I'm doing wrong here.


texts: Keeping dev scripts in sync with projects

texts:
How do you keep your development scripts in sync with your source code management (SCM) and your projects?

For instance right now I have things like D:\code\reusable-stuff that goes to an SCM (e.g. GitHub, Bitbucket, ...). Then in my projects I copy the particular scripts I want to use. However this is a pain and inefficient because it's very rare the script is just used as is. So I tend to avoid keeping things in sync, and rely on occasional big all-in-one lump commits that wave away many benefits of using revision control and inevitably leading to more headaches later.

To illustrate with a single concrete example:

There is a script D:\code\reusable-stuff\conversion\clip_all_layers.py (along with a dozen other scripts in the same folder, none of which I care about right now).

Clip_all is replicated in Project_2015_01\scripts\clip_all_layers.py, which is used in Project_2015_01\Toolbox.tbx\Awesome_Hicky_Model.

I discover an off-by-one bug which meant in some circumstances every 2nd layer is skipped.

I change and test in situ with the project data where the bug was encountered. It takes 3 hours to come to a clean working copy, with a dozen intermediate changes and tests. The first 3 changes are dutifully replicated to D:\code and committed, but the ones between r3 and r12 go bye-bye; too much work. Oh yeah, when committing r12-final from d:\code, another fix was discovered, so we had to replicate back to the project.

2 weeks later, I'm looking at the undocumented jump between r3 and r13 and scratching my head. Important logic is missing and I no longer remember enough to know why solution X was used, and today I need to know that.

There has to a smoother way.


texts: How does PostGIS project from 4978 to 4326?

texts:
I have data in EPSG 4978 (3D Cartesian XYZ using WGS84) and just for kicks, I tried projecting it to 4326 (lon, lat 2D using WGS84). I expected an error but lo and behold, I got a point back.

Here's my query to verify that my point is EPSG 4978:

=>select st_astext(pos_geom) from pos.pos_ts limit 1;
                      st_astext
------------------------------------------------------
 POINT Z (1427635.07972 -5672506.90998 2534091.68628)

Attempt at transformation:

=> select st_astext(st_transform(pos_geom, 4326)) from pos.pos_ts limit 1;
                           st_astext
----------------------------------------------------------------
 POINT Z (-75.8734015076523 23.5640490271681 -20.0468848655)

Like what?! My first reaction was that the output was probably garbage (-20.04688?) but when I use an external library to do the same transformation, it appears that the lon lat values are correct. So my question is twofold:

What is PostGIS doing? AFAIK this behavior, while useful, violates the EPSG standard for 4326. 4326 is a 2D CRS with no Z coordinate in the spec.
Please oh please tell me that the Z coordinate is elevation relative to the ellipsoid! I want it to be true but can't find it documented and my external projection method doesn't return elevation.


texts:
First off: While I have a background in GIS, I only did geocoding once, over two years ago, and the tutorials you find online seem highly specific to their own examples. I'm sorry if what I lay out below is excessively vague but that goes to show how in deep I am I suppose.

I have a CSV that has addresses (i.e. 1234 Wall Street) in one field, and the zip in another. I have clipped my roads down to the geography I am using - made of a few zip codes in NYC.

In ArcGIS 10.3, I tried building out an address locator using both dual, and single range locators. I think that the dual range one is correct, based on the fact that the street data has the "left from" and "left to" fields, etc.

The following few things have happened in my trials:

The geocoding runs, but 100% of results are unmatched.
The geocoding says it has completed, but the status bar never fills.
The geocoding apparently completes, and a shapefile of the points is created, but it is not projected or displayed. Have tried reprojecting but to no avail.


texts: How to avoid rasterized vector to cover minimum area?

texts:
I rasterized vector data(Lines) with QGIS.
The resulting raster covers minimum area, but actually I want it to cover more area for my purpose(See image).

How can I make the resulting raster data cover more area specifically ?


texts:
I'm trying to convert a DXF 2010 file to KML using ogr2ogr:

ogr2ogr -f "KML" outfile.kml infile.dxf

In the file, $INSUNITS is 6, that's meters, if I'm correct?
The problem is that DXF is not georeferenced, and I get errors like:

ERROR 1: Latitude -1113854.310000 is invalid. Valid range is [-90,90]. This warning will not be issued any more
Warning 1: Longitude -684136.440000 has been modified to fit into range [-180,180]. This warning will not be issued any more

I'm pretty new to this stuff, and in the last weeks I haven't found documentation that's explaining what I'm doing wrong. Can I add a suggested lat/lon as a parameter to ogr2ogr? Is it actually possible to do a conversion DXF to KML?


texts: Leaflet.Draw options - dynamically changing polyline color

texts:
I am playing around with the Leaflet.Draw plugin and its options. What I try to do is dynamically setting the color option for the polyline to a variable like this:
let options = {
    position: 'topright',
    draw: {
        polyline: {
            shapeOptions: {
                color: penColor,
                weight: 15
            }
        }
    },
    edit: {
        featureGroup: drawnItems,
        remove: true
    }
};

Unfortunately when changing the penColor variable after the page loads the color does not. Is there a way to change  polyline color after the plugin has loaded, like rereading the options on a variable?


texts: Join dataframe with spatial object and get its spatial boundaries - R

texts:
I am trying to do Spatial + Non-Spatial object join here. However, I am trying to get the spatial coordinates of the matching ids. I tried to use sp::merge, but that is doing the other way around.

My dataframe has ids, I am trying to get the spatial boundaries of those ids alone from spatial polygons data frame.

df looks like,

id
1
2
3
2
1

spatialpolygondf will also have a column called id in its attributes. How can I join them together to get the boundaries ?


texts:
I've wrote below code for visualize precipitation changes over a year with ui.chart command in two different provinces.

var persiann = ee.ImageCollection('NOAA/PERSIANN-CDR')
.filterDate('2017-01-01', '2017-01-30');

var province1 = table2;

var province2 = table3;

var regions = ee.FeatureCollection([
  ee.Feature(province1,{label : 'province1'}),
  ee.Feature(province2,{label : 'province2'})]
  );

var chart = ui.Chart.image.seriesByRegion(
  persiann,regions,ee.Reducer.mean(),'precipitation',1000,
  'system:time_start','label'
  );

print(chart);

however it shows ERROR :

  Error generating chart: Collection.geometry: Unable to use a collection in an algorithm that requires a feature or image. This may happen when trying to use a collection of collections where a collection of features is expected; use flatten, or map a function to convert inner collections to features. Use clipToCollection (instead of clip) to clip an image to a collection.


texts: GEE. : Exported bands must have compatible data types; found inconsistent types: UInt16 and UInt32

texts:
What are the possible reasons for this error?
This is the code:

Map.addLayer(PuntoQuellaveco1);
Map.addLayer(operaciones);
Map.addLayer(Abast);

 var image = ee.ImageCollection('COPERNICUS/S2')
 .filterBounds(PuntoQuellaveco1)
 .filterDate('2016-07-01','2016-08-31')
 .sort('CLOUD_COVER', false);

var conteo = image.size();
print('conteo de imágenes', conteo);

var mejorImagen = ee.Image(image.sort('CLOUD_COVER').first());

print('La primera con menos nubosidad', mejorImagen);

Map.addLayer(mejorImagen, {bands: ['B4', 'B3', 'B2'], max: 3000}, 'image');

var fechaAdquirida = mejorImagen.get('DATE_ACQUIRED');

print('Fecha adquirida', fechaAdquirida);

 Export.image.toDrive({
  image: mejorImagen,
  description: 'imageToDriveExample',
  scale: 10,
  region: Abast
   });

The variables are already created, They just can't be copied because of this:


texts: IndexError: list index out of range on APP STARTUP v 3.16.1, 3.16.0.2, 3.14.0

texts:
During start up of the app, the following error occurs.
An error has occurred while executing Python code:
IndexError: list index out of range
Traceback (most recent call last):
  File "n:/QGIS3~1.16/apps/qgis/./python/plugins\processing\script\ScriptAlgorithmProvider.py", line 112, in loadAlgorithms
    folder = ScriptUtils.resetScriptFolder(folder)
  File "n:/QGIS3~1.16/apps/qgis/./python/plugins\processing\script\ScriptUtils.py", line 114, in resetScriptFolder
    commonSettingPath = os.path.join(paths[appIndex], paths[profileIndex])
IndexError: list index out of range

Python version: 3.7.0 (v3.7.0:1bf9cc5093, Jun 27 2018, 04:59:51) [MSC v.1914 64 bit (AMD64)]
QGIS version: 3.16.1-Hannover Hannover, b381a90dca

Python Path:
n:/QGIS3~1.16/apps/qgis/./python
r:/Users/timv/AppData/Roaming/QGIS/QGIS3\profiles\default/python
r:/Users/timv/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins
n:/QGIS3~1.16/apps/qgis/./python/plugins
n:\QGIS3~1.16\apps\Python37
n:\QGIS3~1.16\apps\Python37\Scripts
N:\QGIS 3.16\bin\python37.zip
n:\QGIS3~1.16\apps\Python37\DLLs
n:\QGIS3~1.16\apps\Python37\lib
N:\QGIS 3.16\bin
n:\QGIS3~1.16\apps\Python37\lib\site-packages
n:\QGIS3~1.16\apps\Python37\lib\site-packages\win32
n:\QGIS3~1.16\apps\Python37\lib\site-packages\win32\lib
n:\QGIS3~1.16\apps\Python37\lib\site-packages\Pythonwin
r:/Users/timv/AppData/Roaming/QGIS/QGIS3\profiles\default/python

What is the cause and/or remedy for the above?


texts:
Bit of an open question but....

Is it possible to design a form in qt designer that will be used as a form against a qgis layer but interacts with a different table in postgres/postgis as a lookup table.

For example, I have a postgis layer that i'm editing in qgis to capture the locations of shops. Some of the fields to capture are address, postcode and parish. Could I do a lookup to another postgres/postgis table to enable these fields to auto-populated my shops form field?

Si


texts:
Here is my CartoDB map:
https://tekim.cartodb.com/viz/1f42eaae-cad0-11e5-8e5c-0ecfd53eb7d3/public_map

I have enabled (Options -> Scroll Wheel Zoom).  Once I do this, I am able to use the scroll wheel to zoom as long as I am logged into CartoDB. However, the people with whom I am sharing the map cannot use the scroll wheel to zoom, and if I am logged out, I can't either. I am using Chrome and I have tested in IE 11 and both exhibit this issue. What might be causing this issue?


texts: Where can I get high resolution multispectral satellite imagery?

texts:
I'm new to GIS, so please be patient if my question is a bit naive. I would like to regularly monitor an area in the Czech Republic in multispectral image depth, with update at least once in a few days. I am willing to pay for it.

I checked the webpages intelligence-airbusds.com and sentinel-hub.com so far; both seem to be promising, but I was not able to get down to details and high-resolution samples.

Could you give me some hints where to start, how to see examples and evaluate the quality, before subscribing and paying?

More generally, which services and platforms would you recommend?


texts: Is there a way to import custom shapefiles for use with Gisgraphy?

texts:
I installed gisgraphy and looked at the data importer wizard -- it only has two options: Import Geonames and/or Import OpenStreetMap.

I want to make Gisgraphy work with shapefiles.  Since Gisgraphy uses PostGIS I was wondering; Is there a way to import shapefiles into PostGres/PostGIS directly and have it work with Gisgraphy?


texts: Getting Error while installing grass extension

texts:
I want to install grass extension r.stream.order in my grass 7.But,while installing I am getting below error:

g.extension -s extension=r.stream.distance url=
Fetching <r.stream.distance> from GRASS GIS Addons repository (be patient)...
Compiling...
Installing...
/usr/bin/install: cannot create regular file
'/usr/lib/grass72/docs/html/r.stream.distance.html':
Permission denied
make: *** [install] Error 1
WARNING: Installation failed, sorry. Please check above error messages.


texts: How to find lowest elevation change using raster data and w.r.t muliple polygons

texts:
I am trying to complete a spatial analysis which yields areas (parcels of land) which are comprised of polygons (in an attribute table) that correspond to the lowest land slope regions or smallest elevation change.

I've got the DEM, and I attempted using Spatial Join Tool to merge the attribute data together to the pink polygons, but I'm not sure how to spatially analyse it in order to yield the results I need.

Can someone please point me in the right direction?

I am using ArcMap 10.3.1.


texts: Converting 3-band raster of a topo map into Unique Values with Colormap or by Reclassifying using ArcMap

texts: I'm using a 1998 USGS topo map for a geologic map base in ArcMap 10.5.1. Right now, the topo raster is in 3-bands with an unsigned 8-bit pixel type. Strictly for cartographic reasons, I need to strip away the green from the topo, and keep the other colors. In order to do that, I converted the raster to an 8-bit signed pixel type by using Copy Raster. Then, I needed to convert the 3-band raster into a single-band raster. For that I used Image Analysis and made the raster grayscale. I can then build an attribute table for the new raster and can apply unique values to the pixels. However, I'm having trouble correlating the values in the new raster to the colors in the original raster. I think I need to build a colormap, but don't know how.


texts: Is it possible to recover ArcGIS Online tiles of a zone, then reproject and republish them in a Map Service?

texts: I'd like to use ArcGIS Online basemaps but not is the Web Mercator projection. Because of that I want to know if it's possible to reproject and republish them in a Map Service.


texts: How do I load 1 billion points into PostGIS?

texts:
I have a file geodatabase containing research data that I need to load into PostGIS. The GDB contains one feature class with one billion records with three columns (the FID, a FK, and the geometry (single points)).

I'm using this to load the GDB into PG:

ogr2ogr -f "PostgreSQL" PG:"host=localhost port=5433 dbname=XXX user=XXX password=XXX" C:\XXX\XXX.gdb -lco SPATIAL_INDEX=NO --config PG_USE_COPY YES

While loading, I also disable auto-vacuum on the table.

After loading the file, I create a spatial index.

The problem I'm encountering is the first 400 million records load very quickly but performance eventually degrades to crawl. It goes from millions of points in mere seconds to taking several minutes to load 100K points. How do I address this or tune my system?

I imported smaller GDBs ranging from 1 million to 500 million records using the above steps. At various points I had over 2 billion records in various tables within the same database. No performance issues and all GDBs loaded w/o issue.

Windows 10
16GB RAM
512GB SSD hard drive, 300GB free

I tried to load a subset of the points using some options I found here, except with ogr2ogr, but performance was horrible even when selecting "FID < 10". I quickly ruled out that option.

ogrinfo out.gdb -al -sql "select * from eea_1Kgrid where fid = 1" -dialect OGRSQL

Thoughts?

update 1: The points will not be rendered all at once. Small subsets may be rendered but this is more about getting all the data into the database for analysis and summary.


texts: Create a definition query for groups of layers: Is it possible?

texts: I have a large number of groups created, with each group having many layers.  I have looked for a way to set a definition query for all the layers in each group but haven't found anything.  Is it even possible, maybe with model builder or python?  I have ArcMap 9.3 and 10.


texts: Import shapefile as table in MS SQL Server (QGIS 3)?

texts: I have made a successful connection with MS SQL Server in QGIS 3.  I can see all of the tables. I know with QGIS 2 and PostgreSQL you had to use DB Manager to import a shapefile as a table.  Well, DB Manager doesn't have the option for MSSQL... Is there an alternative plugin?


texts: Reading and Exporting OSM Road Data as Network of Road Segments [Adjacency Matrices]

texts:
I wish to extract from an OSM file a road network, and export it as an array (adjacency matrix) where the nodes are the road intersections and the edges represent the road segments.

I don't want to just visualise it, but also export the resulting adjacency matrix as a file with a network data structure in at least one of the following formats: numpy file (.npy), pajek file (.net), gephi file (.gexf), csv file (.csv). Alternatively, a list of nodes and edges with geometric info would do.

I seek a solution using either Python or QGIS.

I was hoping for something like this Matlab package, especially the 3rd item in the description. Unfortunately, I don't have Matlab installed on my machine.


texts: Calculating percentage of two hold in two different shapefile

texts: I have two big datasets of different species. One shp is point which shows the location of species, and the other file is in polygon of same places and has got attribute table of public and free holds. I wanted to see how many percentage of the first one (point shp) is in free hold and how many is in public hold?


texts:
I'm trying to make the red text 'un-wrap' itself and use all the room available.

My code is shown below. (I've tried playing with the margin as well but no luck)


texts:
Does Bently have a way to import an image of a highway as point, lines, curves, etc. segments? I.e. digitise an image of a plan?

If something can do most of the hard work it would really help otherwise I have lots of tracing to do in Bentley or 12D model. If I can just get the lines that would really help.

Maybe there is a utility application I could use just to represent it as a dwg then I can bring it into Bentley like that.


texts: How to erase polygon from another in QGIS

texts: I want to erase the pink polygon from the blue polygon. I've tried clip, difference, and symmetrical difference. Clip produces an empty polygon. Difference produces just another blue polygon.


texts: Batch convertion with raw to GeoTIFF images

texts:
I have several 12-bit raw images acquired with Survey3 MAPIR spectral camera. There is a MAPIR Camera Control (MCC) software for make the images batch convertion to GeoTIFF format, but doesn't work very well.

Does anybody know how I could do this convertion in batch, with ANACONDA, QGIS or something else?


texts: Free data with good accuracy of Europe and Mediterranean

texts:
I'm searching for free data like Europe's coastline and Mediterranean island's names etc. but with good precision. I found some shp file from this site
http://www.naturalearthdata.com/downloads/
but the data it was too generalized and the scale small.
Free data usually is without any warranty about their accuracy. I need it for non commercial use.


texts: How do I measure shortest distance where it crosses +/- 180 (the International Date line)

texts:
I'm trying to measure the shortest distance between two sets of points on the globe in QGIS.  One set of points is close to -180 and the other close to +180,  QGIS Distance matrix is giving me a measurement of approx 31,000 KM so is obviously measuring the long way around.  I want the short way around measurement that crosses the international date line.

Any suggestions?  I've seen some people suggest re-projecting the map so that it splits at the prime meridian. However, I have a script that automates this calculation across several thousand data-sets, so I ideally want to solve this in a way where I don't have to re-project.  I currently use WGS 84 / EPSG:4326 and on the fly transformations enabled.


texts:
I want to calculate the cell distances between all cells within a feature class grid. Exactly the same as the game Minesweeper. I have included an example below. The original data contains 11000 cells with different labels but I have simplified in the example:

For each cell I want to find how many cells apart they are. Example here is for cell A13. I would then iterate this for every cell in the grid to produce the final csv.

I have tried cost distance & allocation in model builder but neither will work as the cells are rectangles and not squares which is returning decimal places in the final result.


texts: Find for each point-feature in layer the nearest point in other point-ayer over a network

texts:
I'm looking for a solution for following problem:

given a point layer A (10.000 points)
given a point layer B (15.000 points)
given a line layer C, as network (waternetwork)

I want to have for each point in layer A the nearest point from layer B. That nearest distance must be counted over a network, layer C, and put in a new attribute column in layer A.
I have found some ideas using QNEAT3, but it processes the distances from all points from A to all points from B, which is a huge overload.
Thanks for replies.


texts: Where is the reference for functions available in CartoCSS

texts: I'm playing with Tilemill and i see that for some properties, like marker-transform, i can use functions. I want to use rotate() which takes 3 arguments (i learned this by inputing just one so it reported an error). What i can't seem to find anywhere is the reference for these functions. Can someone point me to it?


texts:
Can I read a raster map which has been constructed in mapserver as WMS with rasterio?

I tried  :

from rasterio.io import MemoryFile

def print_metadata(dataset):
     print(dataset.profile)

url ='http://localhost:8080/?map=/maps/ivm3.map&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&LAYERS=tr_ivm_gray&STYLES=&SRS=epsg:3857&BBOX=2857613.741389,4274927.875099,4989229.633477,5176940.449967&WIDTH=1160&HEIGHT=540&FORMAT=image/tiff'

tif_bytes = open(url,'rb').read()

with MemoryFile(tif_bytes) as memfile:
     with memfile.open() as dataset:
         print_metadata(dataset)

and it gives :

FileNotFoundError: [Errno 2] No such file or directory: 'http://localhost:8080/?map=/maps/ivm3.map&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&LAYERS=tr_ivm_gray&STYLES=&SRS=epsg:3857&BBOX=2857613.741389,4274927.875099,4989229.633477,5176940.449967&WIDTH=1160&HEIGHT=540&FORMAT=image/tiff'

another code :

raster=rasterio.open(url)
print (type(raster))

print ( raster.meta )
print ( raster.count )

array=raster.read()

this gives:

<class 'rasterio.io.DatasetReader'>
{'driver': 'WMS', 'dtype': 'uint8', 'nodata': None, 'width': 1073741824, 'height': 454363579, 'count': 3, 'crs': CRS({'init': 'epsg:3857'}), 'transform': Affine(0.001985222000710666, 0.0, 2857613.741389,
       0.0, -0.0019852220040462356, 5176940.449967)}
3
Traceback (most recent call last):
  File "rio1.py", line 20, in <module>
    array=raster.read()
  File "rasterio/_io.pyx", line 322, in rasterio._io.DatasetReaderBase.read
MemoryError

raster.read() tries to read too much data.
What can I do here?


texts:
I have opened in QGIS a georeferenced raster, a Google Maps layer (from OpenLayers plugin), 1-2 vector layers, and if I save or toggle editing the vector files QGIS freeze, just only the QGIS, the rest applications of my PC work fine.

I saw here some similar questions without a clear answer, but, I have to refer it somewhere. Here is a picture of my system usage while I was trying to turn off the toggle edit.

  PC details: Debian 8, processor Intel i7, RAM 4GB


texts:
I have following elevation data of a region (the higher the altitude, the brigher the points):

How can I find out, which areas are filled by water, and which not?

Can I assume that if the elevation is below X meters, the place is filled with water? If yes, what is X equal to?


texts: Handing a Bounding Box (bbox) to an overpass query

texts:
I would like to incorporate the results of overpass-queries into a model, which in the end gives me a clipped union of these results (all roads within a specific area).

The challenge is to automate the query by handing it a bbox based on the area of interest.

Is there a way to do so? Could "Oriented minimum bounding box" help me in this?

Further specification:

I am using the graphical modeller of QGis to union several different types of roads (motorway, primary, etc) which I downloaded from OSM via the QuickOSM extension for QGis. I am doing several downloads due to possible crashes of QGis when forcing the processing/download of to much data via QuickOSM.

Right now the data has to be downloaded manually, but since I will have to do this several times for different regions I would like to automate the whole process further by including the download into the model.

My question: Is there a way to automatically hand a bbox based on my current area (available as shape) to an overpass API query?


texts: Given two lat/long/time points I need to extrapolate an entire path along the earth

texts:
Given two points and the times between them, I'd like to draw a line around the earth.
Everything I've tried led to some crazy lines, due to my inexperience with longitudinal movements.

Here is the code I'm using:

$orbx = $actualx - 1;
$orby = $actualy;
$steps = 360;
$fullorbit = "$actualx,$actualy,$elevation";
while ($steps-- > 0) {
  $orbx = getDegreesLat($orbx);
  $orby = getDegreesLon($orby);
  $fullorbit .= " $orby,$orbx,$elevation";
  $orbx -= 1;
  $orby += 0;
}
$fullorbit .= " $actualx,$actualy,$elevation";

Of course, my line looks perfect as long as I'm translating only latitude, how can I make this work for both latitude and longitude?


texts:
In ArcMap 10.3 I have a raster file with a polygon grid overlayed on it:

If the pixels in the colored raster have values such as blue = 5, yellow = 1, green = 3 etc, is there a way to output a table that will show the total number of blue, yellow, and green raster cells within the polygon cell (or something like 20 blue/5 cells so blue cells = 100 and so on) . Something similar to tabulate area would be great except I just need the total number/ total value.

This seems like an odd question to me so feel free to ask for clarification. It is difficult to think of all that I need to explain.


texts: Combining Microsoft Excel and QGIS in an elegant way

texts:
I've been working with an Microsoft Excel Spreadsheet where the rows are different investments and the columns are different variables/criteria (a column could for example be the year of the investment, a brief description of the investment or what budget the investment belongs to).
This makes the use of column-filters and pivot tables for reporting/summarising easy.
In addition, most of these investments are georeferenced (most are represented as a polyline and some as points) - however, not all of them are.
The goal is to link the information in the Excel sheet with QGIS so that I have the benefit of the filtering/sorting/etc. of the rows/columns while in Excel as well as exporting atlas maps displaying the investments including key information in QGIS.
So far I've solved this by assigning each investment an unique ID/key in the Excel sheet. I've then made a new sheet named "QGIS" that automatically gathers the most important information from the "main" sheet.
More precisely the columns are:

"ID" (string), unique ID/key for the investment
"2021" (int), budget cost for 2021. (there is column for each year: 2022, 2023, 2024, 2025 and so on)
"description" (string), a brief description of the investment
"budget" (string), what budget the investment belongs to.

When I'm done editing in Excel, I've saved the "QGIS"-sheet as an ".csv" file.
In QGIS I've made two layers: a point-layer and a line-layer. I've assigned each layer an ID column and then I've joined the Excel sheet to each layer on the ID column.
I've then made multiple duplicates of the "master" point-layer and line-layer based on budget categories, investment year, etc.
This works, but the routine is easily prone to errors, and it's very inefficient.
Is there a more elegant way to connect Excel with QGIS?
I'm very interested in how you would solve this problem.


texts:
I got some shapefiles out of a WFS service and the points and polygons are not where they supposed to be.

After taking a closer look in the data I came to the conclusion, that the X an Y coordinates are switched, which can be seen by the extents of the shapefiles.

Is there a way to switch the coordinates back?


texts: Creating subset of TIFF image and generating its NDVI in Python?

texts:
I have two Raster images, one from band4 with a B4 at the end and another with from band5 with B5 at the end. I want to subset the B5 raster to 800x600, then display it and save it as a GeoTiff. Then I want to compute the NDVI (I assume I'll need both the B4 and B5 to do this, but not sure). Then I want to display the NDVI subset of the B5 raster. Display it and save it as a GeoTiff.

How would I create something like a 800 x 600 pixel subset of a TIFF raster image? I want to also take that TIFF and generate an NDVI image for that subset.

NOTE: I am working with a Landsat image. The image has B5 at the end of the file title.

What I've done so far:

import rasterio
from rasterio.windows import Window
import matplotlib.pyplot plt # for later use

with rasterio.open('MyRasterImage.tif') as src:
w = src.read(1, window=Window(0, 0, 800, 600))

This runs fine without error after help from this site.

I want to display this using Spyder or Jupyter notebooks. So I thought to use matplotlib and did the flowing code:

# Plot
plt.imshow(w)
plt.show()

Doing this generates a 800x600 matplotlib window, but it's all purple, not sure why its producing this.

Now I want to be able to display this 800x600 image. Then after that I want preform an NDVI on that subset 800x600 image. Then display the subset 800x600 image with NDVI showing.

I know the formuala is: NDVI = (NIR - red) / (NIR + red)

But how do I extract out NIR and red here from this single Landsat image?

My attempt at that:

band1 = dataset.read(1)
band2 = dataset.read(2)
band3 = dataset.read(3)
print(band[2])

When I run that code for the bands I get the error:

rasterio indexerror: band index 2 out of range (not in (1,))

When I run this code:

print(w.count)

It returns '1'.

So this means that the Landsat image only has one band? But in order to do NDVI don't I need 3 bands?

I am thinking of writing some code like this to get the NDVI from that raster. But not sure how to go about extracting out the bands:

# We handle the connections with "with"
with rasterio.open(bands[0]) as src:
    b3 = src.read(1)

with rasterio.open(bands[1]) as src:
    b4 = src.read(1)

# Allow division by zero
numpy.seterr(divide='ignore', invalid='ignore')

# Calculate NDVI
ndvi = (b4.astype(float) - b3.astype(float)) / (b4 + b3)

This code doesn't work because bands isn't define as anything so I don't know how to define bands to get the NDVI.

After this I am not sure how to go about both displaying and saving the rendered image.


texts: How to write an SQL query to select recently edited features with Editor Tracking enabled?

texts:
Conditions: ArcGIS for Desktop 10.2 and ArcSDE 10.1

Is there a way to select recently edited features using the automatically populated fields that Editor Tracking uses?

The new Editor Tracking feature allows geodatabase managers to track changes to features. There are four new fields that track the creation user/date and the last edit user/date. I'd like to use SQL to identify recently edited ones for other tasks.

The psuedocode would be something like:

SELECT FROM [feature] WHERE [editor tracking field] LIKE [2013-11-04]


texts:
This is related to Utility Network Analyst tool. I would like to build a button to accomplish the following:

A Utility Network Analyst solve-function for downstream with flagging. Ofcourse, the user will set the flag and select the downstream options. When the “new” button is clicked, it will perform the “Solve” function in the Utility Network Analyst tool plus get the objects' information affected by the downstream.

I am very new to Esri. Can anyone help me?

I am using VS2010 C#.


texts: How to identify Unposted/Outstanding (A&D) edits of a version in Enterprise SDE Database using SQL interface?

texts:
I am working on ArcSDE10.2.2 Oracle11G(11.2.0.4) enterprise geodatabase with user base 600+ and 2000 versions( Created  from DEFAULT).
For effective GDB compression and optimized performance delta count need to be controlled and there is user dependency for posting the version edits to default. I want to enlist all versions having outstanding number of edits which has been posted yet. I want to check it using SQL query/python script etc.

Any suggestion inputs please share.


texts: Binding between QGIS custom form (Qt Designer) and date/time field in SpatiaLite database table?

texts:
For my attribute table I created a custom form in Qt Designer. The table has a date field and a time field.

In Qt Designer there are widgets for entering date and time but what data type is expected for the fields in my attribute table in QGIS?

The values always stays NULL in every case, when I choose date, datetime or text. Under Layer properties when I choose the edit widget "calendar" there is a calendar popup even when I use my custom form.


texts: How create String primary key with GeoTools?

texts:
I'm developing a gis app with geotools 19-SNAPSHOT. I want to get some gis data from another app. After getting data, save them in local db. Received data has a String primary key as follow:

@Id
private String code;

First, I want to get featureType,then create table and finally getting data and save them.

When I create a table with SimpleFeatureType, it create int fid primary key, But I want to it create String primary key named code. How Can I do it with GeoTools?

Update:

I create table as follow:

SimpleFeatureType getReceivedFeatureType()
{
     return getDataStore(cfgs).getSchema("myTableName");
}

DataStore getDataStore(cfgs)
{
     URI uri = new URI("db-url");
     Map<String, Object> params = new HashMap();

     params.put("dbtype", "postgis");
     params.put("host", url.getHost());
     //Other Db Connection Configs

     return DataStoreFinder.getDataStore(params);
}

store = getDataStore(cfgs);
SimpleFeatureType sft = getReceivedFeatureType();
store.createSchema(sft);


texts: Reading a RasterIO MemoryFile triggers a GTiff driver warning/error message

texts:
I have recently upgraded RasterIO from version 1.0.5 to 1.0.12. After this upgrade I now get CPLE_NotSupported in driver GTiff does not support open option MODE when reading a MemoryFile.

My code to write then read the memory file is

from rasterio.io import MemoryFile
with MemoryFile() as memfile:
    with memfile.open(**kwargs) as dst:
        dst.write_mask(clust_mask)
        dst.write(bands)

    with memfile.open() as src:
        for i, ea_band in enumerate(tmp_dst.read(masked=True)):
            print (ea_band.mean(), ea_band.std())

The kwargs used to create the file is

{'count': 2, 'crs': CRS({'init': u'epsg:28354'}), 'dtype': 'float32', 'driver': u'GTiff', 'transform': Affine(2.0, 0.0, 300146.0,
       0.0, -2.0, 6181978.0), 'height': 319, 'width': 358, 'nodata': -9999.0}

Is this and error or a warning and should I be worried about it? It seems to write/display the data correctly.
If it's just a warning how do I prevent it from being displayed?


texts: AttributeError: 'NoneType' object has no attribute 'split' from arcmap only

texts:
I'm running a simple python script from arcmap that sends an html post and processes the response. This script runs just fine from pythonwin but for some reason I receive the following error when running the script from an arcmap toolbox: AttributeError: 'NoneType' object has no attribute 'split'

The python code and stack trace is below. I'm stumped why it would run outside of arcmap but not within. I'm sure it's something simple I'm missing. Hopefully someone can point it out. :)

for row in cursor:
        arcpy.SetProgressorLabel("Loading {0}...".format(row))
        url = "https://myurl.com"

        payload = "streetAddress="+row[0]+"&city="+row[1]+"&state="+row[2]+"&zip="+row[3]+"&apikey="+apikey+"&format=csv&census=false&censusYear=2000|2010&notStore=false&version=4.01"

        headers = {
            'accept': "application/x-www-form-urlencoded",
            'content-type': "application/x-www-form-urlencoded",
            'access-control-allow-origin': "*"
            }

        response = requests.request("POST", url, data=payload, headers=headers)
        responseList = response.text.split(',')

StackTrace:

Traceback (most recent call last):
  File "C:\Development\PythonAddins\TAMUGeocoder.py", line 75, in <module>
    response = requests.request("POST", url, data=payload, headers=headers)
  File "C:\Python27\ArcGIS10.1\lib\site-packages\requests\api.py", line 53, in request
    return session.request(method=method, url=url, **kwargs)
  File "C:\Python27\ArcGIS10.1\lib\site-packages\requests\sessions.py", line 456, in request
    prep = self.prepare_request(req)
  File "C:\Python27\ArcGIS10.1\lib\site-packages\requests\sessions.py", line 389, in prepare_request
    hooks=merge_hooks(request.hooks, self.hooks),
  File "C:\Python27\ArcGIS10.1\lib\site-packages\requests\models.py", line 293, in prepare
    self.prepare_url(url, params)
  File "C:\Python27\ArcGIS10.1\lib\site-packages\requests\models.py", line 362, in prepare_url
    host = host.encode('idna').decode('utf-8')
  File "C:\Python27\ArcGIS10.1\Lib\encodings\idna.py", line 157, in encode
    labels = dots.split(input)
AttributeError: 'NoneType' object has no attribute 'split'


texts: QGIS autopopulate attribute table with spatial intersection

texts:
I am working with QGIS.

The aim of the task is to label (auto populate) the 'Road From' and 'Road To' columns in the attribute table based on which roads the street intersects with.

Looking at the example below. Central Ave is divided into segments based on the intersections as highlighted in red below. I want the attribute table to auto populate this row in the Road From and Road To columns to ‘Deakin St’ and ‘Barton St’.

Is it possible to write a code in the field calculator to tell the Road From column to fill according to the ‘Road name’ of the line it intersects with.

I guess the issues (besides can this be done) are that Central Ave intersects with itself as well as Deakin St. So the code will have to specify ‘does not equal itself’. Another one is can you tell it to fill the Road To colum … from the other road it did not put in the Road From Column.

To answer a question I got: No, I am not familiar with python (although when I was researching I did notice people using it) ... and Yes, I did look around at the plugins available. Do you think that Python is the potential answer ?


texts:
I am using geomintersect code in QGIS to only retrieve a target field from a target point layer if the target point layer intersects with the start point or endpoint of the line layer. or Vice versa (retrieving a target field from a target line layer if only the start point or endpoint intersects with the point layer).

I tried to use this function startpoint(geomintersects('targetLayer','targetField')) but I got no result.

what is the most recommended function to retrieve a field when only a point intersects a start or endpoint of a line?


texts: image segmentation of RGB image by K means clustering in python

texts:
I want to segment RGB images for land cover using k means clustering in such a fashion that the different regions of the image are marked by different colors and if possible  boundaries are created separating different regions. I want something like :

from this :

Is it possible to achieve this by K-means clustering?
I have been searching all over internet and many tutorials do it by k means clustering but only after converting the image to grey scale. I want to do it with an RGB image only. Is there any source that could help me begin with it?
Please suggest something.


texts: How to snap a road network to a hexagonal grid in QGIS?

texts:
I'm trying to use QGIS 2.14 to snap a road network to a hexagonal grid, but I'm getting strange artifacts.

I've created a hex grid with MMQGIS, cells are approx 20 x 23 m. I've buffered the road network by 1m and densified it so there is a node every few meters. You can see what I'm trying to achieve below. As you can see, I can get it to work in some cases:-

blue is the densified road (a buffered line)
red is the 'hexified' version - this is what I want to find
the grey is the hex grid

I then used the new Snap geometries feature to snap the nodes to the closest hexagon corner. The results are promising, but there seem to be some edge cases where the line expands out to fill the hexagon (or part of it):-

The reason for the buffer is that Snap geometries doesn't let you snap to a layer whose geometry is different. For example, you can't snap nodes on a LINE layer to points on a POINT layer). It seems to be happiest snapping POLYGON to POLYGON.

I suspect the roads expand out when one side of the buffered road line jumps to one side of the hex cell, and the other side jumps to the other side of the hex cell. In my example, the roads which cross west-east at an acute angle seem to be the worst.

Things I've tried, without success:-

buffering the road network by a tiny amount, so it remains a polygon but is very thin.
densifying the hex cells (so there are nodes along the edges, not just at the corners)
varying the maximum snapping distance (this has the largest effect, but I can't seem to find an ideal value)
using LINE layers, not POLYGONs

I find that if I change to using just LINE layers, It works for a while, then crashes. It seems to save its work as it goes - some lines have been partially processed.

  Does anyone know of any other way to snap points on a line to the nearest point on another
  line/polygon layer, ideally without needing to use postgres/postgis (although a solution with postgis would be welcome too)?

EDIT

For anyone who'd like to have a go, I've put a starter QGIS project here on Dropbox. This includes the Hex Grid and Densified lines layers. (The road network is from OSM, so can be downloaded using QuickOSM e.g. if you need to get the original to undensify the roads).

Note that it's in OSGB (epsg:27700) which a localised UTM for the UK, with units in meters.


texts: In PostGIS: how to split linestrings into their individual segments?

texts:
I have one table with many linestrings and I want to make a new table containing the individual segments of all the linestrings, no matter if they are grouped line per line or not.

This is my table, really simple:

table myline (
  the_geom geometry(LineString,4326),
  id integer,
)

I already tried with no results code like this:

SELECT ST_AsText( ST_MakeLine(sp,ep) )
FROM
-- extract the endpoints for every 2-point line segment for each linestring
(SELECT
  ST_PointN(geom, generate_series(1, ST_NPoints(geom)-1)) as sp,
  ST_PointN(geom, generate_series(2, ST_NPoints(geom)  )) as ep
FROM
   -- extract the individual linestrings
  (SELECT (ST_Dump(ST_Boundary(ll.the_geom))).geom
   FROM mylines ll
   -- eliminate 0 length linestring
    where st_x(st_startpoint(ll.the_geom))<> st_x(st_endpoint(ll.the_geom))
   ) AS linestrings
) AS segments;

and this:

select st_astext(st_makeline(st_pointn(a.line, a.idx), st_pointn(a.line, a.idx+1))) as txt
from (select the_lines as line, traj_id as idx from mylines) as a

I have to say that I created the lines from a set of GPS points I have from the dataset, so another solution could be create segment lines from the points directly.


texts:
I have a erosion map and land use map in DWG files, made of regions (boundary) and hatches.

How to convert regions or hatches to polygons?

This is the picture of the map in DWG file.


texts: How to change the map scale units for ArcGIS Online printable maps?

texts:
I have a web app created from the ArcGIS Online Basic View template. I added print functionality to the web app as seen in figure 1. You can see in figure 2 the scale units are set to miles/kilometers by default in the printable map. Is there a way to convert the units of measurement in the print map scale to feet and/or meters?

Figure 1

Figure 2


texts: No Popup for bottom layers, only top layer works

texts:
I have an OpenLayers Map. On this map, I have three overlay layers.
The first 2 layers come from one GeoServer Layer/data source, with two alternate styles.
The 3rd layer is a separate GeoServer Layer / data source.

For example, I have a GeoSever Layer called 'Roads'. For my 'Roads' Layer,
I have a 'Highway' style to show features with a highway attribute, and a 'Street'
style to show features with a street attribute, but both are from my 'Roads' layer.
The third layer is a called 'Lights', and is a separate set of features.

The issue is that in OpenLayers, If I have only the 'Roads/Highway' Layer visible, I can
successfully do a WMSGetFeatureInfo for my popup. If I only have the 'Roads/Street' Layer
visible, I can do a WMSGetFeatureInfo for my popup.
But, when I have both the 'Roads/Highway' and the 'Roads/Street' layers simultaneously
visible, I can only get a response from the top most drawn layer (painters model). The
other layer beneath the resonsive layer, doesn't give me featureInfo.

Note: when I have either the 'Roads/Highway' or the 'Roads/Street' layer visible with
'Lights' layer, then I can get a featureInfo response from both the 'Lights' layer and
whichever 'Roads' layer I have on.

One further note which I am not sure is relevant, is that I have my GeoServer styles
connected to a Workspace.


texts:
I am trying to publish a model from ModelBuilder so that I can use it on my GIS-Web application, so to start with I have created a model with parameters for buffer, and tested it on ArcMap (Desktop) and it works fine, but when I published it on server, and want to run it, an error appears :

  esriJobMessageTypeError: ERROR 000816: The tool is not valid.

I have attached some captures

Edit 1

I have followed these instructions : 1) I have created a folder that contains a toolbox and two folders, one named scratch and the other ToolData, 2) I have added the shapefiles into scratch folder. 3) I created a file geodatabsae called Schema withing the Tooldata folder 4) I have added the shapefiles into the Schema gdb, 5) I have created the model inside the toolbox. 5) I ran the model and it works correctly but when I publish it, it stop working with the error that I mentioned, I think the problem is with the inputs, I don't know why , but when publishing the service, a message appears to confirm uploading the shapefiles and I accept so the process can be done


texts: How to normalize choropleth maps of startups by state?

texts: I want to display on a choropleth map the startups of a particular domain that operate in that state. I have a count of startups by state. How do I normalize so that the map is not represented incorrectly?


texts: Linear referencing using virtual layers in QGIS

texts:
I am trying to replicate the behaviour explained on this page.
There are two initial layers :

One network layer (line geometry), called base_lines (fields : LB_idRiver)
One event layer (non-geometrical table), called linear_events and containing information about the beginning (eFrom) and the end of the event (eTo) (fields : fid, LE_from, LE_value)

From there, I need to, using virtual layers in QGIS :

Create two point layers showing the position of the beginning and the end of the events along the base_line layer
Create a new line layer extracting the parts of the base_line layer between the points.

So far, my query for the first task looks like this :
SELECT
    e.fid,
    e.'LE_FROM' AS Debut,
    e.'LE_idRiver' AS RouteID,
    e.'LE_Value' AS val,
    r.'LB_idRiver',
    ST_Line_Interpolate_Point(r.geometry, (e.'LE_FROM' - ST_M(ST_StartPoint(r.geometry))) / (ST_M(ST_EndPoint(r.geometry)) - ST_M(ST_StartPoint(r.geometry)))) as geom
FROM "Linear_events" as e
Join "base_lines" as r ON (e.'LE_idRiver' = r.'fid')

The table is created, but the geometry creation fails (NULL gets returned for each part of the geom alias).
What changes have to be made to the query to make it work as intended?
Once this query works, how is it possible to get the second one working?
(On a side-note, I'd be very interested in advice concerning virtual layer creation, as the interface is not really comfortable in QGIS)

I finally made it work based on the advice here :
For the first query :
SELECT
    e.fid AS Event_ID/*:int*/,
    e.'LE_FROM' AS Debut /*:real*/,
    e.'LE_idRiver' AS RouteID /*:int*/,
    e.'LE_Value' AS val /*:int*/,
    r.'LB_idRiver',
    ST_Line_Interpolate_Point(r.geometry, e.'LE_FROM' /st_length(r.geometry)) as geom /*:point:2056*/
FROM  "Linear_events"as e
Join    "base_lines" as r ON (e.'LE_idRiver' = r.'fid')

And the second query :
SELECT
    e.fid /*:int*/,
    e.'LE_idRiver' AS RouteID /*:int*/,
    r.'LB_idRiver' AS RivName /*:text*/,
    e.'LE_FROM' AS Debut /*:real*/,
    e.'LE_TO' AS Fin /*:real*/,
    e.'LE_Value' AS val /*:int*/,
    ST_Line_Substring(r.geometry, e.'LE_FROM' / st_length(r.geometry), e.'LE_TO' / st_length(r.geometry)) as geom /*:linestring:2056*/
FROM  "Linear_events"as e
Join    "base_lines" as r ON (e.'LE_idRiver' = r.'fid')


texts: Edit default properties in "identify features" tool QGIS 2.18.2

texts:
I have a few rasters and want to see their values. When I use the "identify features" tool, I first get the option of "identify all" (see attached image) and only when I click on that I can see the results, but in "Tree" view (that is the default)

Is it possible to change it so I won't have to click on "identify all" every time, so that the results will appear automatically? and that the "view" mode will be set to "Table"?


texts: Font size in HTML in QGIS' print composer only starts to change at 6

texts:
font size 10

font size 9

font size 8

font size 7

font size 6

font size 5

font size 4

font size 3

font size 2

font size 1

What's wrong with my code?


texts: Does osm2po support anchor tags with a colon?

texts:
I have seen some tags in OpenStreetMap that contain colons. Example:

cycleway:right = lane

I have tried to add those tags as anchors "as is" in osm2po.config. For example:

wtr.tag.cycleway:right.lane =    1,  82, 15,  bike

However, I get the following error:

FATAL  Exception at Sat Jul 08 22:19:12 GMT 2017
   Message: java.lang.NumberFormatException: For input string: "right.lane =    1"
   Class: java.lang.RuntimeException
   Stacktrace:
   de.cm.osm2po.converter.Tiler.open(Tiler.java:118)
   de.cm.osm2po.Main.runTiler(Main.java:227)
   de.cm.osm2po.Main.main(Main.java:158)

Is this supported in osm2po 5.1.0? If not, what is the best way to set a clazz for edges with this kind of tag?


texts: Difference between region and polygon in MapInfo Pro?

texts:
I am new to using MapInfo 12 and have recently noticed the convert to regions and convert to polyline function.

What is the difference between a region and a polygon?


texts:
I am trying to create a simple map in UTM32 (EPSG:25832) in R with leaflet and proj4leaflet, but I can't get it to work. I found a url for the basemap, but it could be a different one if you like - I don't really know how to find one.

library(leaflet)
leaflet(
  options =
    leafletOptions(
      worldCopyJump = FALSE,
      crs=leafletCRS(
    crsClass="L.Proj.CRS",
    code='EPSG:25832',
    proj4def='+proj=utm +zone=32 +ellps=GRS80 +units=m +no_defs',
    resolutions = 0.00028*c(
      17471320.7509, 8735660.37545, 4367830.18772, 2183915.09386, 1091957.54693, 545978.773466, 272989.386733, 136494.693366, 68247.3466832, 34123.6733416, 17061.8366708, 8530.9183354, 4265.4591677, 2132.72958385 ),
    origin =c(0.105946948013, 56.8478734515))
)) %>%
  addTiles(
urlTemplate = 'http://sg.geodatenzentrum.de/wmts_webatlasde/tile/1.0.0/webatlasde/default/DE_EPSG_25832_ADV/{z}/{y}/{x}.png',
options = tileOptions(minZoom=0,maxZoom=14,continuousWorld = TRUE)) %>%
  setView(lng = 12.526906, lat = 55.680339, zoom = 14)

I am not quite sure how to set the resolution and origin either.

The example here for UTM33 is working fine.


texts: Write spatial operations to PostGIS from QGIS or ArcGIS

texts:
I need to write my spatial operations to a (new) postgis layer (as the data is to large for a shp file). Can this be done in the QGIS browser or from ArcGIS without SDE?

When you go (for example) --> vector--geoprocessing--clip or geometry--multipart to singlepart you can only save to a shp format eventhough the input data is from postgis.

I have tried the postgis manager and you can edit tables, import/export data but I can't see any way of writing the data to postgis directly.

I have qgis lisboa and postgis 2.5 as well as arcgis advanced (10.1) but not ArcSDE. If necessary a python solution is fine (maybe using @dimitris response in Writing a plugin for Qgis using python to access PostGIS and write geometries

=== POSTGIS Attempt

I am new to Postgis. So as a test I have run the following from http://postgis.net/docs/ST_Dump.html

SELECT (ST_Dump("LowRez_HighHz".the_geom)).geom AS the_geom

FROM
  public."LowRez_HighHz";

It runs with the following output

  Total query runtime: 68839 ms. 47591 rows retrieved.

and then I get an output but has it really removed all the multiparts? I doesn't look like it's writing.

So when I add the following

-- Break a compound curve into its constituent linestrings and circularstrings
SELECT ST_AsEWKT(a.geom), ST_HasArc(a.geom)
  FROM ( SELECT (ST_Dump(p_geom)).geom AS geom
         FROM (SELECT ST_GeomFromEWKT('COMPOUNDCURVE(CIRCULARSTRING(0 0, 1 1, 1 0),(1 0, 0 1))') AS p_geom) AS b
        ) AS a;
          st_asewkt          | st_hasarc
-----------------------------+----------
 CIRCULARSTRING(0 0,1 1,1 0) | t
 LINESTRING(1 0,0 1)         | f

I get

  ERROR:  syntax error at or near "st_asewkt" LINE 11:
  st_asewkt          | st_hasarc
                         ^

  ********** Error **********

  ERROR: syntax error at or near "st_asewkt" SQL state: 42601

How do I correct this?

The table details are

FROM http://www.bostongis.com/postgis_dump.snippet I also tried the following

SELECT "LowRez_HighHz", ST_Multi(ST_Collect(f.the_geom)) as singlegeom FROM (SELECT "LowRez_HighHz", (ST_Dump(the_geom)).geom As the_geom FROM "LowRez_HighHz" ) As f GROUP BY "LowRez_HighHz"

and I get

  ERROR:  out of memory DETAIL:  Failed on request of size 15183219.

  ********** Error **********

  ERROR: out of memory SQL state: 53200 Detail: Failed on request of
  size 15183219.

I have 16GB RAM and 1TB free space...what's happening?

=== Testing Micha's input

I can't seem to run the code in one go. I get

  NOTICE:  CREATE TABLE will create implicit sequence
  "lowrez_highhz_single_gid_seq" for serial column
  "lowrez_highhz_single.gid" ERROR:  Table 'LowRez_HighHz_single' does
  not occur in the search_path CONTEXT:  SQL statement "SELECT
  AddGeometryColumn('','', $1 , $2 , $3 , $4 , $5 )" PL/pgSQL function
  "addgeometrycolumn" line 4 at SQL statement

But when I run

-- Create a new table
CREATE TABLE LowRez_HighHz_single (gid serial, Shape_Leng numeric, Shape_Area numeric);

It creates it as

  NOTICE:  CREATE TABLE will create implicit sequence
  "lowrez_highhz_single_gid_seq" for serial column
  "lowrez_highhz_single.gid" Query returned successfully with no result
  in 60 ms.

When I then comment out the first line and add

-- Make it a spatial table
SELECT AddGeometryColumn('LowRez_HighHz_single','geom',4283,'POLYGON',2);

I get

  ERROR:  Table 'LowRez_HighHz_single' does not occur in the search_path
  CONTEXT:  SQL statement "SELECT AddGeometryColumn('','', $1 , $2 , $3
  , $4 , $5 )" PL/pgSQL function "addgeometrycolumn" line 4 at SQL
  statement

What am I doing wrong? I also tried public."LowRez_HighHz_single"

=== Update on using all lower case.

I tried it but fails...might be because the table has mixed case currently.

the error code is

  ERROR: relation "lowrez_highhz" does not exist SQL state: 42P01
  Character: 327

I also tried it on a dataset that is already all lowercase.


texts: Where to extract Esri ASCII Raster Files?

texts:
In my spare time I am working on interpolating algorithms for Digital Elevation Models in C++. I always obtained my elevation data set as the ESRI ASCII Raster format because it is super simple.

http://resources.esri.com/help/9.3/arcgisengine/java/GP_ToolRef/spatial_analyst_tools/esri_ascii_raster_format.htm

In the past I always obtained my Data from http://gdex.cr.usgs.gov/ which is no longer supported. It was extremely easy for me to obtain this .asc-Files on their site and now it is down. And the new recommended extractor https://lpdaacsvc.cr.usgs.gov/appeears/ seems not to provide .asc-Files. Do you have a suggestion where to obtain this data without being forced to use a third party software to convert from another format to ESRI ASCII Raster format.


texts: What's the meaning of marker-fill-opacity: 0.1*[value]?

texts:
I'm using heatmaps in CartoDB and I would like to modify CartoCSS to reflect the value of a column. I have noticed this line in CartoCSS:

marker-fill-opacity: 0.1*[value];

I don't understand completely what [value] refers to but maybe there is a way to point to my column, something like:

marker-fill-opacity: 0.1*[my_column];

Is there something like this?


texts: Offset Google markers at the same address

texts:
What can I incorporate into my current code below to offset markers, as markers at the same address don't show.

Here's example http://www.mediwales.com/mapping/test

UPDATE:
I've changed the way the map works by running it from a database.


texts: How to align with QGIS two layers that have the same projection?

texts:
I am using QGIS 2.12.0 to create forest maps. I usually succeed to create maps from shapefiles with different projections but for some reason I cannot succeed to align the layer of one shapefile.

As described in the help file and solutions found over the Internet I:

Assigned a projection using the right-click button
Assigned a similar projection to the project
I turned on “on the fly” retroprojection and reloaded the shapefile that I cannot align
Replaced the projection of the other layer with the one of the shapefile that I cannot align
Started my map using the shapefile I cannot align
Erased the .prj and the .qpj files to force QGIS to use the projection I want
Created a new copy (using the "save as" option) of the shapefile that I cannot align using the projection I want

I obviously check if QGIS accepts the change in projection. In all cases, QGIS accepts the changes.

What else should I do to succeed?


texts: How to create Circle using X,Y Coordinate and radius by ArcPy?

texts:
I have been given three parameters X and Y coordinates and a distance or the radius of the circle.
I will have to draw a circle using ArcPy by these three parameters where the center is the point of X-Y coordinate and radius is distance.

e.g: Center (20.5937° N, 78.9629° E)
     Radius 5 unit


texts: How can I find internal overlaps in my polygon feature class ideally via modelbuilder tools?

texts:
I have a polygon layer containing "depth curves" for the making of a skin of earth kind of map. I want to check whether or not my feature class contain any overlaps internally between these polygons.
What kind of tools can I use for this? Ideally i would like to end up with a feature class containing the bits of overlaps, so i have the possibility to check graphically where my features overlaps. And if I could somehow implement this in a larger model, things would really look good.


texts: How to create a Typographic Map with Mapnik

texts:
So i am trying to make a typographic map using mapnik. I want the street names along a particular street not to overlap, but I want an intersection between two streets to over lap. Currently some streets don't have repeating labels. see attached image for whats happening. How do I accomplish this?

labelStyle = mapnik.Style()
text = mapnik.TextSymbolizer(mapnik.Expression("[FULL_NAME]"), "DejaVu Sans Book", 2,     mapnik.Color("#171717"))
text.label_placement = mapnik.label_placement.line
text.label_spacing = 2
text.force_odd_labels = True
text.min_distance = 0
#text.max_distance = len((mapnik.Expression("[FULL_NAME]"))
text.character_spacing = 1
text.minimum_path_length = 0
text.label_position_tolerance = 100
text.allow_overlap = 1
text.avoid_edges = 0
r.symbols.append(text)
#labelStyle.rules.append(r)
s.rules.append(r)
m.append_style('Text',s)