texts: Intersection of a raster and polygon in ArcGIS 10.3

texts:
I have two data sets one raster and one polygon. I am interested to see learning the area where the two overlap. I was thinking that I could use the intersect function, but since one of the layers are a raster that will not work.

Is there another way for me to compute the area of the overlying areas?


texts:
this is my first question here and I have done extensive searches but still can't find the (full?) answer I'm looking for. Hopefully someone can help!

I am loading a number of personal geodatabases (mdb) to PostGIS using OGR2OGR (and GDAL). I am able to get the whole database to load with no issues, but for my requirements I need to store the original filename alongside each record in each table (so that I can query/delete/filter at a later stage).

At the moment I am doing this using a simple OGR2OGR script that loads all tables in the MDB (over 250 of them)...this is about the limit of my ability right now so any solution needs to be fairly simple if possible?

My current simple script:
ogr2ogr -f "PostgreSQL" PG:"host=** active_schema=** user=** dbname=** password=** port=5432" *.mdb -skipfailures -append -a_srs EPSG:27700

Thanks.


texts: Switching career from web application to GIS developer?

texts:
I am a software developer with over 20 years building applications mostly in the web application space for various companies.   I am trying to get a foot hold into the GIS software development space and I am looking for advice on the best way to achieve that.

I've been contemplating obtaining a GIS Certificate and taking classes using the latest GIS tools out there.    I was hoping someone could make a suggestion as to an option or options I may take to do so.


texts: QGIS crashes when digitizing new fetaures

texts:
I can't currently digitize polygon features in QGIS 2.01 or 2.3.
I can draw a set of points but as soon as I right click to finish the system crashes. It still works fine in 1.8.

Windows 7 64


texts: If it is possible to create an animation in ArcGIS Pro with buildings rising up ( having the building footprints and DSM of the area as available data), then how do I do it?


texts:
I'm trying to modify the icon example from to display multiple icons, each with its own set of attributes that display when the icon is clicked.  The issue I'm having is that if icon1 is clicked, and then icon2 is clicked, before clicking anywhere else, icon2's popover contains the contents of icon1's popover.

I've tried various things, like running the popover('destroy'), but nothing seems to work.

I've created a jsfiddle to illustrate the problem, and included a hacky fix that I've commented out.


texts: Adding WMTSLayer from service to ArcMap using ArcObjects?

texts:
I want to add all the layers in my WMTS service to ArcMap using ArcObjects, but I am able to get only one layer.

Does anyone have a solution for this?


texts: Python OGR SQL - Finding invalid geometry

texts:
I previously asked this question@ regarding finding files that contain an invalid geometry using OGR.

The following SQL command was suggested to me, which seemed to work a treat: -sql "select * from filename WHERE not ST_IsValid(geometry)"

I am now trying to get this to work under python using the following code:

from osgeo import ogr

driver = ogr.GetDriverByName('MapInfo File')
dataSource = driver.Open(file, 0)
layer = dataSource.GetLayer()

result = dataSource.ExecuteSQL("select * from file WHERE ST_IsValid(geometry)")
featureCount = result.GetFeatureCount()

print (featureCount)

Upon running the code however, I get the following error message: ERROR 1: Undefined function 'ST_IsValid' used.

Does anyone know how I go about checking ST_IsValid under python?


texts: ArcHydro catchment area of specific stream segment

texts:
I am given the task to create the catchment areas of specific streams/lakes. To sum it all up, I am looking to find the catchment areas for lake 1, stream, lake 2 and downstream of lake 2 all individually.

When looking at the stream, I only want the stream's watershed minus that of lake 1. In other words, that which falls from roughly start to finish. There should be 4 catchments areas: one for each element previously mentioned and in the following picture.

Is this possible at all? Here are the steps I have gone through:

Initial data: DEM, Stream Polygon Layer, Lake Polygon Layer

Available resources: ArcGIS 10.5 basic, Spatial Analyst, ArcHydro plugin

1- DEM leveling: I've modified the raw DEM slightly so the lakes are a flat surface and not a slightly odd surface.

2- DEM reconditioning: modified the DEM with linear feature for, what i understand, is slightly change it's structure to allow water to pass through known streams with ease.

3- Fill sinks: Filled sinks in the DEM to avoid water getting "stuck" in processes.

4- Flow direction

5- Flow accumulation

6- Stream definition (this is where I think I might not be getting it): The stream definition takes as input the flow accumulation, which is a grid containing cells with a value in cells flowing into them.

There is a user-defined area for the area size I'm guessing, which I used 50km2. I've tried with many values all not giving results I might want so I continued with 50km2 to see if something would change and give me a result I would like.

7- Stream segmentation: This gave me a grid with cells having unique values per stream, fairly straight-forward.

8- Catchment grid delineation: Creates catchment areas with the grids created through the previous steps 6 and 4. The results vary depending on the area chosen in step 6. The smaller the ares, the more catchment areas I get.

9- Catchment Polygon processing: Creates a polygon feature class with the catchment areas created in the step 8.

10- Drainage line processing: Converts my stream link grid into drainage lines (doesn't give me many drainage lines. I think I might get more with a chosen area for step 6 lower than 50km2.

11- Adjoint Catchment Processing: This tool did not work for me. I get an error. I'm thinking it might be because I have too little catchment areas and they already are as "joined" as they can be for this tool to work properly. That is my guess, it could definitely be something else.

Here is the error i get when running this tool:

And here is the rest of the error code (it was too long for the screenshot)

Now, is everything I'm doing here in vain? Are the steps followed here ultimately going to lead me to my objective explained previously of getting the catchment areas of specific locations? I feel like I am building a tool but i can't quite grasp where it is going.


texts: Proportional Coloring within polygon for multiple variables

texts:
I'm new to QGIS and I am doing some demographic analysis of Ethiopia. I have a vector map with each of the country's administrative zones:

Within the attribute table, I have data for religious demographics in each zone.

What I would like to do is have color swatches appear within each polygon (five different colors each representing the different religion), based on the proportion of the population that follows that religion. So, if one zone was 50% Protestant and 50% Muslim, that polygon would have equal amounts of red and green swatches.

Similarly, if one zone was 30% orthodox, 50% muslim and 20% protestant, there would be three color swatches present, with the color representing muslim being predominant.

Here is also a screenshot of my attribute table, as requested:

The info I'm interested in displaying is to the right: the percentage of a zone's population that follows a certain religion. For many of the areas, it is overwhelmingly a single religion, but a good amount show a lot of variation over two to three different religions.


texts:
I am working with a DTM of 5x5m for a large area (already clipped it). I calculated the slopes and want to maintain only the slopes above 2 degrees. I already used the raster calculator to set all values above 2 to 1 and the rest to 0. However when I convert to a vector, I get a way too large vector file which results in only empty shapefiles once I try to do anything with it in the geoprocessing tools.

What I need is to be able to get an intersect (or clip) with another shapefile.

Is there a way to dissolve the different pixels with the same value before vectorizing the raster or some other way to be able to get an intersect between this way too large raster and my vector file?


texts: failed gdal script from QGIS Georeferencer using Thin Plate Spline (TPS)


texts: How to convert existing lines to polygons with Arcobjects, preferably in C#?

texts: Does anyone know how to convert existing line geometry to polygons with arcobjects as it works with the geoprocessing tool "feature to polygon"?


texts: Generating the skeleton of polygon using PyQGIS

texts:
I develop a plugin for my project which includes some useful tools. I reviewed some answers in GIS SE about generating skeleton of polygon. I couldn't find any solution which used PyQGIS. HCMGIS plugin has a tool "Skeleton / Medail Axes" but it uses only processing tools.
I want to generate the skeleton of a polygon using pure PyQGIS without using Processing tools. After doing some digging I found voronoiDiagram method defined in QgsGeometry. I don't figure it out how to use it correctly.
I used voronoiDiagram in the following script, but it throws an exception: AttributeError: 'list' object has no attribute 'disjoint'.
What am I doing wrong? How can I improve/change this code to get skeleton of polygon?
My code is:
# get polygon layer
polygonLyr = QgsProject.instance().mapLayersByName("polygons")[0]

# create new layer
centerlineLyr =  QgsVectorLayer("Linestring?crs=EPSG:3857", "Centerlines", "memory")
QgsProject.instance().addMapLayer(centerlineLyr)

# layer has one polygon
polygonGeom = list(polygonLyr.getFeatures())[0].geometry()

# create voronoi lines
lineCollection = polygonGeom.voronoiDiagram(edgesOnly=True).asMultiPolyline()

#select lines and add them to layer
centerlineLyr.startEditing()
for line in lineCollection:
    if not line.disjoint(polygonGeom):
        f = QgsFeature()
        f.setGeometry(line)
        centerlineLyr.addFeature(f)
centerlineLyr.commitChanges()


texts:
On QGIS 3.10, after I saved and close my file, and then reopen it.

The layers are unidentified, so I browse for the file location (shapefiles).

When I press on browse, the folders that appear is not the location of the shapefiles.

In other words, where does QGIS store the shapefiles?

I searched every folder of QGIS in program files, and I saw not of the Python shapefiles.


texts: OpenLayers drag event slow to respond (hard to move)

texts:
We have an application that loads data from the postgis and draws on the map in layers of vectors. When many elements are loaded, does the map get heavy to use the drag event, slow to respond to user interaction? Does it have to do with the amount of elements loaded on the map? I have seen applications with much more elements drawn but not so heavy for navigation with the drag event. The problem is only in the Drag event.

var DragZoom = new ol.interaction.DragZoom({
    condition: ol.events.condition.always
});

Load Feature on map (we have many:
$.ajax({
            url: '/Map/GetPolygonsByQuadrante/',
            beforeSend: Processando("Obtendo dados..."),
            data: JSON.stringify({
                MinLat: MinLat1,
                MinLong: MinLong1,
                MaxLat: MaxLat1,
                MaxLong: MaxLong1,
                NivelZoom: NivZoom
            }),
            cache: false,
            async: true,
            type: "POST",
            dataType: "json",
            contentType: "application/json; charset=utf-8",
            success: function (dataFromServer) {
                ClearLayerFromMap(VectorPolygon);
                MessLbl("lblMessage", "Ok!");
                try {
                    var estilo;
                    Polygons = new Array(dataFromServer.length);
                    // Loop over each object
                    for (var $i = 0; $i < dataFromServer.length; $i++) {
                        var ativo = GetVisibleByLayer(dataFromServer[$i].Tipo);
                        if (ativo === true) {
                            var object = dataFromServer[$i];
                            //verifica se o o layer existe
                            //adiciona a linha num array do mapara para posterior identificação
                            var pol = {
                                Id: object.Id,
                                Tipo: object.Tipo,
                                Pontos: object.Pontos
                            };
                            Polygons[$i] = pol;
                            var coordinates = null;
                            coordinates = [object.Pontos.length + 1]; //tem que voltar ao polygono noi ponto inicial
                            for (var $t = 0; $t < object.Pontos.length; $t++) {
                                coordinates[$t] = [object.Pontos[$t].Longitude, object.Pontos[$t].Latitude];
                            }
                            coordinates[length - 1] = [object.Pontos[0].Longitude, object.Pontos[0].Latitude];
                            estilo = GetStyleByLayer(object.Tipo);
                            var init = 'pol-';
                            DrawPolygon(init + object.Id, estilo, coordinates, SourcePolygon);
                        }
                    }
                    ProcessandoClose();
                    MessLbl("lblMessage", "Tot. polígonos:  " + dataFromServer.length.toLocaleString('pt-BR'));
                    ProcessandoClose();
                    AtualizandoMapa = false;
                    map.render();
                } catch (e) {
                    MessLbl("lblMessage", "...");
                    ProcessandoClose();
                    TrataErro(e);
                }
            },
            error: function (json) {
                if (json !== undefined > 0 && json.responseText === "secao_invalida") {
                    window.location.href = "/Home/";
                } else {
                    MessLbl("lblMessage", "...");
                    ProcessandoClose();
                }
            }
        });


texts: USGS Earthquake GPS Metadata: Buffer Zones & Python

texts: I'm making a python web crawler that gets USGS earthquake GPS data, and I'm wondering how I should process the latitude longitude pairs to see if it falls within a specified radius of a specified GPS location. What is the mathematical model---pseudocode really---that needs to be used for this problem?


texts: Trouble using ogr2ogr -sql commands from QGIS to rename attribute fields

texts:
I have a large number of files whose attribute fields are named consistently, but incorrectly. I would ordinarily use ogr2ogr.exe to fix this from the command line using something like:

ogr2ogr "C:\Temp1\OutDataSet.shp" "C:\Temp1\asdf\InDataSet.shp" -sql "SELECT End AS LineEnd, begin AS Begin FROM InDataSet"

This creates a new .shp file, renames attribute fields from End to LineEnd and from begin to Begin , also drops many unwanted attribute fields from the original shapefile.

I was hoping that I could achieve the same thing from QGIS 2.18 processing toolbox > [OGR] Miscellaneous > Execute SQL, as It is easy to set up a batch by using the QGIS forms. Has anyone had any success with this tool for this task?


texts:
I tried to process the following query with PgAdmin III, but it doesn't work and the hint provided doesn't make sense to me :

UPDATE grille_segments_finale
SET shipping_cost = 1000000
WHERE st_intersects(grille_segments_finale.geom::LineString, countries.geom::MultiPolygon) is true ;

And I get:
ERREUR:  le type « linestring » n'existe pas
LINE 6: WHERE st_intersects(grille_segments_finale.geom::LineString,...
                                                         ^

********** Erreur **********

ERREUR: le type « linestring » n'existe pas
État SQL :42704
Caractère : 191

Could you help me please? The geometry "LineString" is the right geometry type of my table ("grille_segments_finale") so I can't really find out where I am wrong...


texts: Keeping original edge IDs in Output Feature from sDNA analysis?

texts:
I’d like to keep the original edge IDs for each link in processing a network with sDNA using the Integral Analysis function. I’m using the free version of sDNA in QGIS, and from what I understand, new edge IDs are produced in the Output Feature from the Integral Analysis function.

I'm processing a network derived from OSM and would like to keep the original IDs for further analysis.


texts:
When I georeference satellite images with a CRS for my region (Hartebeeshoek94), I get a scale of 1:10 and larger. But when I use WGS84 the scale appears correct, e.g. 1:50 000.

I am using QGIS v2.13.3. For the Transformation settings, I use Thin Plate Spline for Transformation Type, Linear for Resampling Method and Hartebeeshoek94 for the Target SRS.

What I am doing wrong?


texts: C# Abstraction for Web Feature Service (WFS) Call

texts:
I'm currently working on a WFS client-side implementation in C# for a thick client application (Windows desktop application). The goal is to be able to choose some "layer" combinations to be shown overlaid on a map. The rendering will use existing function calls in the thick client. The C# / .NET requirement is pretty firm.

At this stage I'm looking at using standard .NET function calls to make HTTP GET calls to GetCapabilities, standard XmlReader to parse the capabilities, build an abstraction of that, then make more HTTP calls based on attributes from that abstraction. Its not that hard (probably because I haven't got to the GML parser part yet...) but I wondered if there might be a better way. In particular, I'm concerned that my client might work for some servers and not others (based on some experiences implementing WMTS in OWSlib).

I wondered if anyone knew of a reasonably "alive" open source product / project that can provide a C# abstraction for WFS on the client side.

There have been similar questions (in particular Connect to Web Feature Service using python/C#, which has a pretty rubbish answer with some nice follow-up banter; but its not very specific in its requirements) and some that looked superficially similar, but dealt with server side implementation or were actually for an admin aspect (e.g. GeoServer REST).

The nice people at Mass GIS have made available a C# connector, which is basically what I'm looking for, but it hasn't been updated in a long time (2007, by the looks). An ongoing project based on that would be interesting.

I'd prefer something that didn't "bundle" the rendering part (to make the final part a bit "lighter"), but if it wasn't that hard to extract, it might be interesting too.

In the ideal world, there would be OWSLib for C#.


texts: Downloading GIS Map for Pakistan into QGIS

texts:
I'm looking to add canals and bunds to my raster and have found some GIS maps at:

Canals

Bunds

Bunds Image

Network Plan Map Image - looks to contain both canals and bunds

I'm not clear on how to download or create map or shapefile of these maps. Is there a way to do this in QGIS? Or is there some other way to open the maps on QGIS?

A similar problem at World map with clear difference between land, sea and inland waterways, essentially I am looking to add inland waterways/canals and bunds and just need a good resource for maps containing this I could use in QGIS, if not a way to use the maps I have found.


texts: Keeping unchanged observed values after kriging interpolation

texts:
Is there a way to force kriging to keep unchanged measured values of observation points?

For example: if the observed value of a location is 1, after the interpolation it could be 0.9 or even 0.5, depending on the other locations' values. I'd like to have 1 for that observed location after the prediction too.

I know that, considering the kriging formula, this is a stretch but predicting values where it is actually known seems to me a stretch as well. I found that the only prediction method that preserves observation points' values is IDW but this is not the best method to predict values for the other prediction locations.

I'm using R and the package "intamap".

I'm using the "automatic" method provided by "intamap" package. I need to use it for two reasons: the first is for the very strong skewness of my data and the second reason is that kriging interpolation is the final step of method I'm developing that is devoted to be user friendly for inexpert persons.

"Intamap" package claims to meet these two needs. The reason why I need to keep original data unchanged is that observation points represent the probability values of a phenomenon to occur. Pay attention, I'm not talking about indicator kriging, because observed values are just probabilities. So I need to provide inexpert users the geographical distribution of different degree of probabilities. If I have two observation points with 1 and 0.5 probability values, I need these values to be unchanged and that kriging interpolates the values in the prediction locations between them.

I hope to have been clear in my explanation. How can I proceed now?


texts: qgis2 composer - non-rectangular inset map?

texts:
I'd like to have a near-full-page rectangular map object, then add a non-rectangular inset map object with 6 vertices.  Any suggestions on how to accomplish this, other than making two actual maps then brute-force post-processing the images?
Thanks!


texts: Adding consecutive unique ID for distinct value sets and sorting by date

texts:
I have a layer with features, which have a unique ID, called "Object_id". A second column is "Date", a third called "info" (see below).

What I need is to fill the fourth column "NEWID" with unique numbers for each set of features that have the same "info"-value. I already found solutions to do that in these posts: Post 1 and Post 2. But: I need to sort the "NEWID" ascending by the "Date"-column, just as shown in the table above (the lower "Date"-values need to get the lower "NEWID"-values).
The example above was filled freehand, just to clearify how it should look...


texts: How to geocode location using twittersearch for leaflet map of tweets

texts:
I'm using twittersearch to iterate over user input topics and location of interest that will then be visualized in a leaflet map. So using 'place' I can get a map, no problem. I'm trying to increase my odds of getting points by using 'location' in another for loop. Do I need to geocode the location data, if so what am I missing?

Geocode piece..

from geopy import geocoders

def geo(location):
    g = geocoders.Nominatim()
    loc = g.geocode(location)
    return loc.latitude, loc.longitude

original 'place' piece that produces results..

for tweet in ts.search_tweets_iterable(tso):
    tweetText = ( '@%s tweeted: %s' % ( tweet['user']['screen_name'], tweet['text'] ) ).encode("ascii", "ignore")
    print tweetText
    if tweet['place'] is not None:
        (lat, lng) = geo(tweet['place']['full_name'])
        tw = tweet['text']
        user= tweet['user']['screen_name']
        #map has markers placed at the location of the tweet, the pop up text is the contenet of the the tweet and the poster/user
        #this is in the form of username said tweet content
        folium.Marker([lat, lng], popup= user + " said " + tw,icon=folium.Icon(color='red',icon='info-sign')).add_to(mapLeafletPython)
        print( '@%s tweeted: %s' % (user, tw) + 'And they said it from (' + str(lat) +', ' +str(lng)+')').encode("ascii", "ignore")

this is the location piece that produces this error...AttributeError: 'NoneType' object has no attribute 'latitude'
which occurs on the (lat1, lng1) line.

for tweet in ts.search_tweets_iterable(tso):
    tweetText = ( '@%s tweeted: %s' % ( tweet['user']['screen_name'], tweet['text'] ) ).encode("ascii", "ignore")
    if tweet['user'] is not None:
        (lat1, lng1) = geo(tweet['user']['location'])
        tw = tweet['text']
        user= tweet['user']['screen_name']
        #map has markers placed at the location of the tweet, the pop up text is the contenet of the the tweet and the poster/user
        #this is in the form of username said tweet content
        folium.Marker([lat1, lng1], popup= user + " said " + tw,icon=folium.Icon(color='blue',icon='info-sign')).add_to(mapLeafletPython)
        print( '@%s tweeted: %s' % (user, tw) + 'And they said it from (' + str(lat) +', ' +str(lng)+')').encode("ascii", "ignore")

"tweetSearchUser" is the input from user as well as "where"

last piece that I think you need for testing besides variable ts which is the token..

tso = TwitterSearchOrder() # create a TwitterSearchOrder object
tso.set_keywords([tweetSearchUser]) # let's define all words we would like to have a look for
#tso.set_count(200) #limit return records to a number
tso.set_include_entities(False) # and don't give us all those entity information

Using this as a reference for twitter api and json format

http://socialmedia-class.org/twittertutorial.html

https://github.com/MichaelMartinez/blog/blob/master/content/iPython-Notebooks/Python-Twitter-Tools.md

I'm baffled here...

Updated with error...


texts: Creating a linestring from points using shapely

texts:
I am creating a linestring from points using shapely, but I could not generate the linestring using the points on my local system. This code is running fine on Google Colab. Is there any package dependency which I don't have installed on my system?
import shapely
from shapely import speedups
from shapely import geometry
from shapely.geometry import shape, Point, LineString, Polygon

a = LineString([(0, 0), (1, 1), (1,2), (2,2)])
a

Error massage
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-18-02826fee6746> in <module>
----> 1 a = LineString([(0, 0), (1, 1), (1,2), (2,2)])
      2 a

c:\python37\lib\site-packages\shapely\geometry\linestring.py in __init__(self, coordinates)
     46         BaseGeometry.__init__(self)
     47         if coordinates is not None:
---> 48             self._set_coords(coordinates)
     49
     50     @property

c:\python37\lib\site-packages\shapely\geometry\linestring.py in _set_coords(self, coordinates)
     95     def _set_coords(self, coordinates):
     96         self.empty()
---> 97         ret = geos_linestring_from_py(coordinates)
     98         if ret is not None:
     99             self._geom, self._ndim = ret

c:\python37\lib\site-packages\shapely\speedups\_speedups.pyx in shapely.speedups._speedups.geos_linestring_from_py()

ValueError: GEOSGeom_createLineString_r returned a NULL pointer


texts: Header row in Excel is not header row when imported into ArcMap?

texts:
I am importing an Excel table with X/Y coordinates to plot on a map (point data).

The original Excel table has a header row as its top row..The row is 'frozen' in Excel.

When I import the table and view it in ArcMap (10), the header row of the original Excel sheet is gone, and the first feature row is now the header row in ArcMap.

How can I avoid this?


texts: Splitting feature into layers which are overlapped in same layer in QGIS

texts:
My script needs to select each feature in layer. While selecting  features I found some of feature are overlapped. Feature can be selected but cannot  highlight in map canvas.
Is there any plugin or tools which can be used to separate overlapping feature into new layer so that one layer has all non overlapping feature and another layer has all overlapping features.


texts: How do I programatically edit the layer 'Display Expression'

texts:
I'm wondering if someone has or knows where a tool is to edit the 'Display Expression' of a Layer, or more specifically, a TableView, in ArcGIS Desktop. I know that it can't be accessed with Python, but it should be editable with ArcObjects.  But this isn't something I know how to do yet.

I'm using a Python Add-in to add a Layer into ArcMAP which has a relate with a couple TableViews that I add at the same time.  The trouble is that it always resets the field that is displayed in the identify window to something unhelpful.

Thanks for any suggestions.


texts: How can I convert a MULTILINESTRING geodatabase file into a coordinate system with lat / long for individual points in a CSV file?

texts: I have tried ogr2ogr, and QGIS but nothing seems to work for converting a MULTILINESTRING into individual lat/long coordinates.  I am doing this in order to get it into a flat file so I can rearrange the data for conversion into a data format that requires specific columns for latitude and longitude.  Ideally it would be in CSV format.


texts: Multiple python import from feature class into feature dataset in geodatabase

texts:
Is it possible to create a python script to import shapefiles/feature classes, from different sources into specific feature datasets in geodatabase, but in one script?

I managed to write the script to make multiple feature datasets in geodatabase. This is the code:

>>> import arcpy
    arcpy.CreateFileGDB_management("D:/GIS_Temp", "TEST.gdb")
    fdList = ["Datase_A", "Dataset_B", "Dataset_C"]
    for fd in fdList:
    arcpy.CreateFeatureDataset_management("D:/GIS_Temp/TEST.gdb", fd, "D:/GIS_Temp/Projection.prj")

And now i would like to import shapefiles/feature classes into specific datasets. For example, from "folder A", to "Dataset_A", "folder B" to "Dataset_B", and so on.

I could do it one script by one, like this:

>>> import arcpy
    arcpy.env.workspace = 'D:/GIS_Temp/Folder_A'
    arcpy.FeatureClassToGeodatabase_conversion(["File_X", "File_Y"],
    'D:/GIS_Temp/Test.gdb/Dataset_A')

But is it possible to do all it in one script. I am really new to python scripting, so I really don't know to combine those separate scripts and do the task.

If I use this code, i can get shapefiles into root of database, but not in feature dataset:

>>> import arcpy
    from arcpy import env
    import os
    env.workspace = "D:/GIS_Temp/Folder_A/"
    fcList = arcpy.ListFeatureClasses()
    for fc in fcList:
    arcpy.CopyFeatures_management(fc, "D:/GIS_Temp/Test.gdb/" + os.sep + fc.rstrip(".shp"))

And if I add name of dataset like this, it won't run:

arcpy.CopyFeatures_management(fc, "D:/GIS_Temp/Test.gdb/Dataset_A" + os.sep + fc.rstrip(".shp"))

This could work:

>>> import arcpy
    arcpy.env.workspace = 'D:/GIS_Temp/Folder_A'
    arcpy.FeatureClassToGeodatabase_conversion(["File_X", "File_Y"],
    'D:/GIS_Temp/Test.gdb/Dataset_A')

>>> import arcpy
    arcpy.env.workspace = 'D:/GIS_Temp/Folder_B'
    arcpy.FeatureClassToGeodatabase_conversion(["File_X", "File_Y"],
    'D:/GIS_Temp/Test.gdb/Dataset_B')

but only from one folder, to one dataset. I would like to write ONE code that would allow me to import defined shapefiles/classes, into specific feature datasets.


texts: Identifying X-Y coordinates from a raster

texts:
I need to find the high and low values of multiple raster datasets individually. Then determine the x,y location of these values.  I am using ArcGIS 10.1, have built the mosaic datasets I need, and have determined the locations of the high and low values (estimations are fine).  However I cannot determine the X,Y location for these values (points).

I know there is a "Find x-y" tool, but I need the opposite.  Is there a way to point click and return the x-y for that location?

With vector data I could populate a field with x-y values, however I do not have the time to cut out the raster datasets individually and turn them into vector data to do this.

There must be another way?


texts: ArcPy find points that are paired in time and space

texts:
I have a shapefile that that is time stamped and tracks two animals. I would like to create a field that denotes "together" based on proximity and location. And example of the data is:

time       animal        data
13:00:04      A            x
13:00:02      B            x
22:07:03      A            x
23:08:04      B            x

Together is defined as being within 20 meters of each other and points falling within 5 seconds of each other. So the first pair would be classed as 'together' but the second pair would not, assuming all these points were in the same location. I have looked for similar questions that use ArcPy but could not find any.


texts: How to store lines and polygons in JSON documents?

texts: Looking at the increasing NoSQL movement  and considering that databases like MongoDB offers a new perspective in flexible data storage for GIS. What is the best way to store lines and polygons in JSON documents to take advantage of 2d indexes and spatial functions?


texts: We sent client Esri files and now they want to see polygons

texts:
We are a voter list company. Every ten years new political boundaries are written for federal and state offices. All of the corresponding shapefiles are ESRI ready.

One of our clients is getting voter data from us and wanted to supplement it with shapefiles. We do not use ESRI ourselves nor do we understand the raw data and polygons as we are quite illiterate in this area---we simply are giving the raw ESRI data to an end user.

We sent the folder of shapefiles to a client who is also unschooled in ESRI and now they want to "see the Polygons". See sample data that we set the client at:

http://www.srrproject.org/resources/redistricting-shapefiles/#shapefiles-content.

My questions is the same as I posed to them: can they not "see the polygons" themselves from the shape data that we provided? Educate me on how these shapefiles are are to be employed.

NEW as of Oct 10:
Hi @Kirk, @RyanDalton, @Rayner, @PROBERT, @Aaron, I understand this is all free advice here in this wonderful exchange but I want to say that if what I am asking is above and beyond free, we can discuss offline a small consultative gig for you by contacting us at: anne.horrigan@eMerges.com. But here is my specific challenge: The client has 50,000 complete addresses across the US. They want to link each address with its corresponding congressional, state senate and state house district. I am looking for step by step---aka, For Dummies--- instructions on how the client may accomplish this task on a regular basis using the shapefiles mentioned earlier. For example, what software, what to load, how to import & export.


texts: Trying add mixed geometry types to same drawingTools GeometryLayer, get Error: Geometries parameter must be an array of ee.Geometry or GeoJson objects

texts:
var collection = ee.FeatureCollection(mycollection);

var myGeometry = collection.filter(predicates).geometry();

myGeometry.evaluate(function(geomList) {
        // Add the geometry as a layer of the drawing tools.
        drawingTools.addLayer([geomList], 'myGeometryLayer', 'red');
      });

  Error: Geometries parameter must be an array of ee.Geometry or GeoJson
  objects

Error only occurs if myGeometry is of more than one type (in my case, Points and MultiPolygons)

How can I add myGeometry to the drawingTools layer when it is of mixed type please?


texts:
I am working on an extension to OpenLayers to pre-cache tiles in a local Web SQL datastore to support mapping on mobile device when connectivity is dead.

I would like to pre-seed my cache with a set of tiles. Currently this has to be done by visiting the area in question and using the cacheWrite control.

This would be done by fetching the tiles in question from my MapServer WMS  and storing them locally in the Web SQL database. These would then be available for the cacheRead control to use on tileerror.

So, my question is: given a bounding box, is there an easy way for me to construct these calls using OpenLayers and hence populate the cache with the required tiles for a required zoom level?

A sample single tile call is:-

"http://mymapserver/pub/mapcache/mytheme?LAYERS=my_layer&FORMAT=image%2Fpng&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&STYLES=&SRS=EPSG%3A27700&BBOX=533371.66,169851.66,534011.66,170491.66&WIDTH=256&HEIGHT=256"

I can do this by panning the map about automatically, but the user experience is terrible ;-)


texts: Easiest way to create measured intervals between shape points

texts:
I need to import a scanned photo of the outcrop (geology project) and re-create the grid from which I took plug data- each plug point will be drawn as a shape point file.

The plugs were taken on a metre by metre grid layout. What would be the best way to ensure that the points are evenly spaced on the outcrop photos?- I guess I could mark the points on the scanned photo and then shapefile the marked points.


texts: Inputting x and y coordinates into table using geometry tokens via update cursor?

texts:
I am brand new to Python and am having trouble with this objective in my Intro to Python for GIS assignment:

Independent step… Create two new fields in the attribute table named X and Y    (or something handy).   Populate that field with the X,Y location of the oilspill   points themselves!

Note: The oilspills featureclass has many points already with set geometry but the coordinates are not given in the existing table.

I have successfully created the new columns called X and Y using the following code:

import arcpy

oilspills = 'S:\\Lab 2\\oilspills.shp'

arcpy.AddField_management(oilspills, "X", "DOUBLE", "", "", "13", "", "NON_NULLABLE", "NON_REQUIRED", "")
arcpy.AddField_management(oilspills, "Y", "DOUBLE", "", "", "13", "", "NON_NULLABLE", "NON_REQUIRED", "")

My issue now is that I can't figure out how to populate those new columns with the x and y coordinates for each point in the table. That said, I have managed to print out all the x and y coordinates for each point in the featureclass via the Python shell:

import arcpy

with arcpy.da.UpdateCursor('s:\\Lab 2\\oilspills.shp',['SHAPE@X','SHAPE@Y']) as cursor:
    for row in cursor:
        print row[0]
        print row[1]

But I'm not sure how to combine the first part of my code and the second so the new columns I've created get populated with those values. How can I do this?


texts: Showing one feature with two unique values as one symbol in ArcMap legend?

texts:
I am working in ArcMap 10.4. I have one feature that has two unique values, one is "simple fill" blue the other value is "10% simple hatch". I need to show the 2 values in the map but would like to only show the solid hatch in the legend. I know I could separate the feature and then just remove one from the legend.

Is there another way?


texts: Vertex tool is converting multi linestring into multi curve

texts:
After editing a PostGIS MutliLineString layer in QGIS with the vertex tool I can not commit the table edits because the MultiLineStrings have been converted to MultiCurves. The error message is:

Could not commit changes to layer road_cl_segmented
Errors: ERROR: 1 geometries not changed.
Provider errors:
PostGIS error while changing geometry values: ERROR:  Geometry type (MultiCurve) does not match column type (MultiLineString)

QGIS 3.16.3, PostGIS 3.1
I can not change the PostGIS geometry type. How can I prevent the vertex tool from changing the geometry type?


texts: Inserting table record without geometry into feature class using Python?

texts:
I am working with ArcGIS and Python and trying to append table rows into a feature class.

You cannot use append and you cannot just make a table view of the feature class and then try to append it in.

Any ideas?


texts: arcpy.da.SearchCursor storing data in a list with a query

texts:
I am currently working with a database that hold street names, x and y coordinates and many more fields. I would like, through a search cursor, to store the x,y coordinates in a separate list for every street name (unique string).

For example, I have 23 entries under Washington street (with 23 x and 23 y coordinates). I would like to find a way to create a list named washington and store the coordinates in that list, and do that for every single street in my database inside a for loop.

For now, all I got is 1 big list with all the data :

spatialjoin = arcpy.env.workspace + '\SpatialJoin_test.shp'

coord=[]

cursor=arcpy.da.SearchCursor(spatialjoin, ['Road_Name','X','Y'])

for row in cursor:

    print(row)

    coord.append(row)


texts:
I have never coded in Python before and I am trying to loop a simple map algebra equation across 365 rasters.  I am trying to calculate total water flow over cells but in the raster I have a negative value to indicate the water table is below ground: for my purposes those cells should read 0.  I have tried to write a for loop but it does not seem to be working.  Python does not throw any errors, but the output data set still has negative numbers. Any suggestions on why the loop isn't working to remove the negative values?

arcpy.env.workspace="C:/..../2013_q1_tiff_v2prov"
Fcs= arcpy.ListRasters()
From arcpy.sa import *
from arcpy.sa import *
print Fcs

for Raster in Fcs:
    Con(Raster <0, 0, Raster)

pleasework= CellStatistics(Fcs, "SUM")


texts: How to update a field in a CSV file/table based on where it is on the map

texts:
I have a shapefile which I downloaded that displays different region types. I also have a CSV file of geocoded addresses. I have loaded both into QGIS and can see the addresses scattered about on the map in different regions. I want to update my address file with the region code based on where it shows.
I had done something like this years ago in mapinfo using a query and object contains syntax, but not sure how to do this in QGIS.


texts: ArcGIS Server on Windows Server 2008 - .NET or Java Backend

texts:
I am finding myself having to "defend" why I have chosen to setup ArcGIS Server on .NET backend with IIS (vs Java and Tomcat).  The host server is running Windows Server 2008 R2.

I thought I remembered reading some type of ESRI slick stating that ArcGIS Server with .NET and IIS offers better performance vs Java and Tomcat, at least on Windows.  But after a couple hours of searching, I can't find any definitive evidence.

Can anyone state any hard, cold facts as to why ArcGIS Server on .NET and IIS is better on Windows?  Links, comments, personal experiences welcome.


texts:
I am creating a Leaflet application that uses the code from the Simple CSV template to set the bounds of the map to the selected markers.

The issue I have is that I want to set a start position that isn't the bounds of all of my markers. I would also like the map to return to this start position once the clear function is initiated as well.

Some of the key sections of code can be seen below but the repository is available by the link above.

var addCsvMarkers = function() {
hits = 0;
total = 0;
filterString = document.getElementById('filter-string').value;

if (filterString) {
    $("#clear").fadeIn();
} else {
    $("#clear").fadeOut();
}

map.removeLayer(markers);
points.clearLayers();

markers = new L.MarkerClusterGroup(clusterOptions);
points.addData(dataCsv);
markers.addLayer(points);

map.addLayer(markers);
try {
    var bounds = markers.getBounds();
    if (bounds) {
        map.fitBounds(bounds, {maxZoom: 15});
    }
} catch(err) {
    // pass
}
if (total > 0) {
    $('#search-results').html("Showing " + hits + " of " + total + " Projects");
}
return false;
};

&

$("#clear").click(function(evt){
    evt.preventDefault();
    $("#filter-string").val("").focus();
    addCsvMarkers();
});

EDIT: Is it possible for the setView that is called in when defining the map to override the fitBounds command?

var map = L.map('map', {
zoomControl: false,
minZoom: 3,
maxZoom: 18
}).setView([54.5, -4], 6);


texts: Google maps API find the nearest markers in my location

texts:
I have this code and I don't know how to include "find the nearest markers in my location".

I tried some codes but nothing didn't work perfectly.

I want if i click to button it shows markers in my location and remove all other markers. Is it possible?

function initMap() {
            var map = new google.maps.Map(document.getElementById('map'), {
              zoom: 3,
              center: {lat: 49.7437895, lng: 15.3360726}
            });

            var infoUser = [
                ["Italy"],
                ["Germany, Czech Republic"],
                ["Czech Republic"],
                ["Poland"],
                ["Prague, Czech Republic"],
                ["Italy"]
            ]
            var markers = [];
            for (var x = 0; x < infoUser.length; x++) {
              var adresa = infoUser[x][0];
              $.getJSON('http://maps.googleapis.com/maps/api/geocode/json?address='+adresa, null,
                function (adresa) {
                    return function(data){
                    var p = data.results[0].geometry.location;
                    var latlng = new google.maps.LatLng(p.lat, p.lng);
                    var marker = new google.maps.Marker({
                      position: latlng,
                      title: adresa,
                      map: map
                    });
                    markers.push(marker);
                  }
              }(adresa));
            }
        }


texts: Creating coded values with iteration from a field

texts: I'd like to find a way to add coded values to a Text domain using values from a table field, if such a thing is even possible. I'm using the Add Coded Value to Domain tool connected to two variables: Name and Description. While this hasn't worked yet, what I've tried was making the value of each variable the same as the two field names from my table with the name and description I'd like to have. My syntax for Name was %VendorName% and my syntax for Description was %VendorDescription%. All this does is create a string literal %VendorName% coded value. Is what I'm trying to do possible?


texts: What are the coordinates used in the Path parameter of the ArcGIS Polyline class?

texts:
I'm drawing a bunch of carrier routes from the USPS API to an ArcGIS map. I need to only show the routes that are within a certain distance from a certain point. I have the latitude and longitude for the center point.

I thought the paths used by Polyline were made up of latitude and longitude sets multiplied by 100000, and I tried using them as such but that didn't work. The paths are generated from point sets that look like this:  -9214713, 3256075. When divided by 100000 they look like longitude and latitude, but I guess they're not.

What are these coordinates used in the Path parameter of the Polyline class, and how can I convert them to lat and lng for the purpose of calculating distance?


texts:
I've been searching for the fastest way to load multiple CSV Files into QGIS 3.10.
My CSV-files consist of 5 columns:
Brand Name | Branch Name | Address | Latitude | Longitude

I am using Mac, and "files:///" might not work. Any suggestions?
I saw code related to this, and followed this. Here's what I came up with:
import glob, os

csv_files="Users/ajie/Documents/RestoBrand_Py/SAMPLE CSVs/"
os.chdir(csv_files)

for fname in glob.glob("*csv"):
   uri = "files:///" + csv_files + fname + "?delimiter=%s&crs=epsg:4326&xField=%s&yField=%s" % (",",",",",","x","y")
   name = fname.replace('csv', '')
   lyr = QgsVectorLayer(uri, name, 'delimitedtext')
   QgsProject.instance().addMapLayer(lyr)

Here's the error


texts: How to show selected feature top in map?

texts:
As shown in the piucture, there is a bunch of multilines in one layer. I selected one feature from the attribute table and wanted to look at its contour.

But the contour is overlapped by other features in this layer.

I'm looking for an option to bring the selected feature to the top in the upper picture. I'm using QGIS 3.10


texts:
I'm struggling with a rather odd problem with the python package OSMNx.
I have lots of GPS points in a pandas dataframe that were originally in EPSG 4326, but I've already projected them onto EPSG 32633 (UTM zone 33N). My aim is it to find the nearest OSM edge. For the sake of an example, let's look at the point (388860.7103015079, 5828739.389188443). One can look at this point here (might not work - maybe you have to enter the point manually).
Now I am using OSMNx and its get_nearest_edge() (and get_nearest_edges(), respectively) to find the nearest OSM edge. As the graph input, I use a rather small bounding box including the coordinates, as retrieved by, e.g.,
G = osmnx.graph.graph_from_bbox(52.607584, 52.573916, 13.357185999999, 13.424399999999999, network_type='drive', simplify = False),
but I guess it should also work, if you'd get the whole Berlin map.
G = osmnx.graph_from_place("Berlin, Germany", network_type="drive", simplify = False.
Normally, I have a whole pandas data frame of points and I run
osmnx.get_nearest_edges(G, df["lon_proj"], df["lat_proj"], method="kdtree", dist=0.0001)
For example's sake, let's run instead
osmnx.get_nearest_edge((388860.7103015079, 5828739.389188443))
Using the latter command, I end up with a tuple of OSM nodes specifying the edge I should be looking for: (1873509499, 732609150, 0). Now there are two problems. First, the OSM edge depends on the graph I am giving to that function (but note that both graphs G stated above contain the example point). Second, the nodes attained by this output are not really close to the original point (compare here).
Why won't the function give me the actual nearest edge, even though there are nodes much more close-by than the ones I am attaining? In the picture below you can see the node that is my desired output (orange) and the node attained by get_nearest_edge().
Edit: I figured that apparently osmnx.get_nearest_edge() takes a tuple of (lat, lon) coordinates (so latitude goes first), while osmnx.get_nearest_edges() does it the other way round. Now using the former works, but using the latter still gives super random results. Perhaps you guys have an idea about this.


texts: Creating TIN in standalone ArcPy script for ArcGIS Pro?

texts:
I'm trying to create a TIN from contour data. My program will actually create an empty TIN, but it will still give me an error. Below is my code:

import arcpy
arcpy.env.workspace = r"C:\Users\Aliza\Desktop\GIS Programming Fundamentals\Final_Project\Contours"
arcpy.CheckOutExtension("3D")
out_tin = "SampleTinssas"
height_field = "Elevation"
in_features = r"C:\Users\Aliza\Desktop\GIS Programming Fundamentals\Final_Project\Contours\Site_37_Contours.shp"
sr = arcpy.Describe(in_features).spatialReference
arcpy.CreateTin_3d(out_tin, sr, "Site_37_Contours.shp", "DELAUNAY")

arcpy.CheckInExtension("3D")

I know that the in_features portion of the tool is supposed to be in the form of a value table, but I have no clue how to format that. I know that my height needs to be included, but again, I have no idea where to start with that.

Below is my error message:

"C:\Program Files\ArcGIS\Pro\bin\Python\envs\arcgispro-py3\python.exe"
C:/Users/Aliza/PycharmProjects/finaltest.py
Traceback (most recent call last):
 File "C:/Users/Aliza/PycharmProjects/finaltest.py", line 8, in <module>
    arcpy.CreateTin_3d(out_tin, sr, "Site_37_Contours.shp", "DELAUNAY")
  File "C:\Program Files\ArcGIS\Pro\Resources\ArcPy\arcpy\ddd.py", line
 2922, in CreateTin
    raise e
  File "C:\Program Files\ArcGIS\Pro\Resources\ArcPy\arcpy\ddd.py", line
2919, in CreateTin
    retval =
convertArcObjectToPythonObject(gp.CreateTin_3d(*gp_fixargs((out_tin,
spatial_reference, in_features, constrained_delaunay), True)))
  File "C:\Program
Files\ArcGIS\Pro\Resources\ArcPy\arcpy\geoprocessing\_base.py", line 506, in
<lambda>
   return lambda *args: val(*gp_fixargs(args, True))
arcgisscripting.ExecuteError: ERROR 999999: Error executing function.
Failed to execute (CreateTin).

Process finished with exit code 1


texts: "Scale range" vs "Scale dependent visibility" when labeling features using QGIS

texts:
What is the difference between "Scale range" vs "Scale dependent visibility" settings when labeling features using QGIS?
From my perspective they are identical.


texts: How do I change a raster's geotransform with GDAL?

texts:
I'm trying to use GDAL to set the geotransform of a new raster that I've created by following the steps outlined in the GDAL API Tutorial.

# create the new dataset
driver = gdal.GetDriverByName('GTiff')
dataset = driver.Create('test_gt.tif', 60, 60, 1, gdal.GDT_Float32)

# check the default geotransform
print dataset.GetGeoTransform()  # prints (0, 1, 0, 0, 0, 1)

# try to alter the geotransform and ensure that it has been set
dataset.SetGeoTransform([0,1,0,0,0,-1])
print dataset.GetGeoTransform()  # prints (0, 1, 0, 0, 0, -1)

dataset = None  # closes the dataset

# Try reopening the dataset now and see if the geotransform has been set.
ds = gdal.Open('test_gt.tif')
print ds.GetGeoTransform()  #prints (0, 1, 0, 0, 0, 1)

The dataset.SetGeoTransform() documentation says that this should set the affine transformation coefficients (which, according to the dataset.GetGeoTransform() documentation is set to [0, 1, 0, 0, 0, 1] by default), but as you can see from my above code, the changes don't seem to actually take effect when I try to change them.

I've even tried flushing the new raster's cache to disk by dataset.FlushCache(), but this doesn't seem to save the changes either.

How can I have GDAL actually save the altered geotransform to disk?

I'm using GDAL version 1.6.3 installed from PyPI (via the command-line tool pip).


texts: Service Hosted for Arc Server- Editor Tracking Only Allowed for UTC Time?

texts:
I am hosting a service to my server then sharing it to multiple apps/AGOL
I want to track edits so I enabled editor tracking with database_time since all my editors are in my timezone and local.
However when I go to publish the service I get the error:

Error: "SEVERITY STATUS CODE DESCRIPTION NAME TYPE DATA FRAME
High [Re-analyze] 00129 Data frame has layers that record Editor Tracking dates in database time zone, but the service does not have time zone information Layers Data Frame Layers"

Is UTC only allowed for services? if not how do I use database time since it is already enabled?


texts: Two layers with same projection do not overlap

texts:
I have a boundary layer for an area, the layer properties are:

top: 4431527.680136 m
bottom: 4431199.984995 m
left: 569103.516495 m
right: 569438.079758 m

Projected Coordinate System:    NAD_1983_UTM_Zone_15N
Projection: Transverse_Mercator
False_Easting:  500000.00000000
False_Northing: 0.00000000
Central_Meridian:   -93.00000000
Scale_Factor:   0.99960000
Latitude_Of_Origin: 0.00000000
Linear Unit:    Meter

Geographic Coordinate System:   GCS_North_American_1983
Datum:  D_North_American_1983
Prime Meridian:     Greenwich
Angular Unit:   Degree

Now, I have a new line shapefile (that I received from someone else) that is supposed to overlay on top of my layer but it did not. The new shapefile did not have any coordinate system. I defined it with the GCS North American 1983. I tried projecting it to NAD_1983_UTM_Zone_15N but it does not work. Here are the layer properties for the new layer. I am confused about the extent in dd.

top: 5194.424316 dd
bottom: 4802.068848 dd
left: 4754.654785 dd
right: 5162.077148 dd

Geographic Coordinate System:   GCS_North_American_1983
Datum:  D_North_American_1983
Prime Meridian:     Greenwich
Angular Unit:   Degree

Does anyone have any suggestion?


texts: Removing house numbers from OpenStreetMap using OpenLayers?

texts: I want to remove all the house numbers of the base layer from OpenStreetMap using OpenLayers 2 and if it is posible the points of interest too.


texts: ArcMap Spatial Analyst - Kernel Density Estimation Produces Blank Rasters

texts:
My problem is with the Kernel Density operation in Spatial Analyst. ArcMap produces blank rasters for some reason when given a normal point feature class as input. I had no problem creating the heatmap I wanted in QGIS and then adding the raster back into Arc, but would like to avoid this inefficient process in the future. My inputs are attached. The same happens with the Point Density operation. Any ideas?

I have ensured all layers are the same projection as the data frame.


texts: How many active QGIS users are there in total?

texts:
I'm the CEO of MangoMap which is a online GIS data management and web map publication platform, we're big QGIS fans and make extensive use of it. We would really love to integrate MangoMap more tightly with QGIS by making the transfer of data and edits completely seamless between the desktop projects in QGIS and MangoMap data/maps in the cloud. We already allow QGIS to get layers from MangoMap using WMS but want to make it much easier to send data and edits from QGIS to MangoMap.

Before we get started though I would really like to know if there are any up to date numbers out there on how many active QGIS users there are at this time. I've seen an article on Spatial Galaxy from 2011 that estimated 100k, but it's my feeling that QGIS has matured and gained a lot of traction since then so the figure is unlikely still valid.


texts: Does GeoPandas offer some dedicated Equal Distance projection like it is available for example in FME during buffering (https://www.safe.com/transformers/geographic-bufferer/) or in ArcGIS when applying metric buffer?


texts: How to display variable number of images associated with one point on QGIS

texts:
I have a set of point locations. I have a set of images with reference numbers. I'd like a way to easily see the images associated with the location, using QGIS.

This is for my own use rather than necessarily for communicating to others, so it should stay on QGIS if possible and the solution should be 'live' (or at least should be very easily updated) rather than involving a one-off import/export of any kind.

At the moment these are vector images, but could easily be raster. On this occasion it's images of sign designs, but in the future might be something else - perhaps photos.

A key issue is that one location has anything from 0 - 5 images associated with it. I'm expecting to produce a text file associating locations with image references (i.e. with two columns, location reference and image reference). It's not beyond the bounds of possibility that I could something more imaginative (database etc). I'm beginning to explore the capacity of QGIS to work with one-to-many relationships. It looks like for ArcGIS the answer might lie with python code based labelling - is there an equivalent?

The issue is that I don't really have a clue what the options might be so can't investigate them, nor even how to ask this question in a really clear and directed way... but I'm sure that others have the same question.


texts: Calculating average NDVI for GPS coordinates using QGIS?

texts:
I have calculated the NDVI image and overlain GPS coordinates for a pack of wild animals over the top.

I would like to calculate the average NDVI of all the GPS points, and then calculate the average for all of the random points that I have also put over the top. I would like to do this for the entire map, and Zonal Statistics does not seem to be helping.


texts: gdal_translate JPEG to GeoTIFF using Python

texts: I was able to convert my GeoTIFF images (from the SpaceNet dataset) to JPEGs using gdal_translate but I would like to know how to convert the JPEGs back to GeoTIFF. I've noticed that gdal_translate also creates an auxillary file with extension jpg.aux.xml when I do the conversion. Is this to be used for the conversion from JPEG to GeoTIFF? If so, how can this be done?


texts: Group attribute tables values for use with conditionals - Python

texts:
I would like to group values in a field in an attribute table. The field contains duplicate values. These duplicate values should be grouped based on duplicate values, the field would look something like the arbitrary 'Field' below:

Field   Date    Flagged_Field
X      1/1/2014
X      1/1/2013
X      1/1/2015    Y
Y      1/1/2000
Y      1/1/2008    Y
Z      1/1/1999
Z      1/1/1998
Z      1/1/1992
Z      1/1/2003    Y

Values X, Y and Z which are duplicate records would then be grouped. Ideally I would like to proceed to use if, elif and else to compare other fields within the group containing the value X for Field and once finished move on to the group contain the value Y and lastly Z.

Here's what I have:

cursor = arcpy.SearchCursor(fc)
row = cursor.next()
for row in cursor:
    rowField = (row.getValue("Field"))
    listC = ([list(j) for i, j in groupby(rowField)])
    for groups in listC:
        print str(groups)

What happens is I get a list of unicode for each letter in the field. I believe if I can somehow compare the field values against each other as opposed to the characters in just one field I can get close to what I need with this.

I'm trying to gather all like values based on attributes, specifically one field. I would then like to, using these like values compare other attributes of the data. For example I have a date field which, based on like values in the Field shown in my above example I would like to flag the record containing the most recent date, see revisions above.


texts: How to calculate and compare RMSE between two DEMs?

texts:
I'm doing a project based on using Lidar point cloud data and I have access to a DEM that cover all of Sweden. But my goal is to produce several DEMs covering a small area using different interpolation methods. For this I'm working with the LAS dataset toolbox. When that's done I want to compare my own produced DEMs with the DEM that cover whole of Sweden and check for RMSE error.

I'm quite clueless how to solve this since I'm not a mathematician and I'm new to working with Lidar.

I'm using ArcGIS for this.


texts: Field Types in ArcMap, QGIS and their respective modelers

texts:
I am aiming to replicate an ArcMap Model in QGIS Graphical Modeler.

ArcMap model functions give a much wider range of options (see fig 1) for the type of field to be created. Furthermore, the field calculator function when accessed through the GUI offers different options than those offered in the Model Builder function. I am creating a model in which these field types are of particular importance. I need to ensure the field I am creating using the QGIS Modeler field calculator function is Double.

Primary Question: How can I ensure I replicate a "Double" field type in the QGIS Model Builder Field Calculator function?

ArcMap Model Field Type

QGIS Model Field Type for Field Calculator Function

QGIS Field Calculator Types


texts: Finding area of three overlaying rasters in ArcMap?

texts:
I have to find the area of three overlaying rasters which are DEM,LULC and Total Evapotranspiration at the same time so that I can find the landuse landcover area and evapotranspiration of land at a particular elevation.
I have used Tabulate area when just 2 of these overlaying rasters come but don't know which tool to use when three rasters are coming at the same time.


texts:
I want to anticipate that I know nothing about programming and I am a new university student in the world of QGIS.

Today I tried to turn on QGIS on my mac(Finder version: 10.13.6) to continue to work for a lab at university, but during startup it crashed saying: 'Unexpected quit of QGis'.

I tried to uninstalled it completely (Version: 3.4.13) and reinstalled it but happened anyway, also doing the same process with the last version of QGis downloadable (Version: 3.10 I think), but the result was always the same.

Can someone help me pls? I'm desperate.

The complete error text is quite long.
The following link is a Word document with the error text, on my google drive:

https://drive.google.com/open?id=1l4lZKFLAvkF0vyiEJRp0MM4ocn6hrEeL

I attach some screenshots of it, if it can help somehow.


texts: Split line at newly created point feature

texts: I have one Point layer and Multiple line layers. I always wants to Maintain Topology between them. For which, If user adds any new feature in Point Layer, exactly snapped with line feature, the line feature should split at that point. I checked the manual way like 'split feature' tool available, where I need to split line first and place my point feature to that location. I checked Class: QgsSnappingUtils, to identify if the point is snapped to line before splitting it but not getting any work around. Please suggest how to proceed. Thank you.


texts: Starting Tile Sever Rendering Error: OSMBright.xml

texts:
I was following the entire process found here: https://switch2osm.org/serving-tiles/manually-building-a-tile-server-14-04/

I have an error though on Testing your tileserver part where i run the following commands

sudo -u username renderd -f -c /usr/local/etc/renderd.conf

Encountered this error: renderd[6166]: An error occurred while loading the map layer 'default': Could not create datasource. Required parameter 'type' is missing  encountered during parsing of layer 'land-low' in Layer at line 48 of '/usr/local/share/maps/style/OSMBright/OSMBright.xml'

What have I missed?

Edit: Here are the lines in OSMBright.xml, the first line here is line 46 so line 48 is the line with srs

<Layer name="land-low"
status="on"
srs="+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0.0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs +over">
    <StyleName>land-low</StyleName>
    <Datasource>
        <Parameter name="dbname"><![CDATA[osm]]></Parameter>
        <Parameter name="extent"><![CDATA[-20037508.34 -20037508.34 20037508.34 20037508.34]]></Parameter>
        <Parameter name="file"><![CDATA[/usr/local/share/maps/style/osm-bright-master/shp/simplified-land-polygons-complete-3857/simplified_land_polygons.shp]]></Parameter>
        <Parameter name="srs"><![CDATA[+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0.0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs +over]]></Parameter>
    </Datasource>
</Layer>


texts: Edit attribute table for shapefiles in folders and subfolders with arcpy?

texts:
Continuously to Merge shapefiles in folders and subfolders with arcpy I try to edit attribute table for 20 shapefiles that spread in big a folder and are divided to a lot of sub folders. All shapefiles are called "migrashim" and the name of the filed that i want to change is "structuretype". I try this code:

import arcpy,os,sys,string,fnmatch,arcpy.mapping
from arcpy import env

rootPath = r"C:\Project\layers"
pattern = 'migrashim.shp'
features = arcpy.UpdateCursor(r"C:\Project\layers\*.shp")
counter = 0
for root, dirs, files in os.walk(rootPath):
    for filename in files:
        if filename == "migrashim.shp":
            for feature in features:
                if feature.structuretype == 3:
                feature.structuretype = 4
                features.updateRow(feature)
        counter = counter + 1
print 'edit'
print counter
del feature,features

but i got an error:

IOError: "C:\Project\layers\migrashim.shp" does not exist

What wrong with my code?


texts: Running PyQGIS script directly (without using interface)?

texts:
I wrote the following script, which runs fine from the QGIS python console (it loads two shapefiles, intersects them, and then finds the areas of the new regions):

import qgis.utils
import qgis.core
from qgis.analysis import *
from PyQt4.QtCore import QVariant
import os
print(os.getcwd())

layer1 = iface.addVectorLayer("/Users/casta2k/ShapefilePractice/precincts/USA_precincts.shp", "precincts", "ogr")
layer2 = iface.addVectorLayer("/Users/casta2k/ShapefilePractice/data_EPSG_4326/Bnd_2015_q1_region.shp", "zipcodes", "ogr")

overlayAnalyzer = QgsOverlayAnalyzer()
overlayAnalyzer.intersection(layer1, layer2, "/Users/casta2k/ShapefilePractice/intersect/intersectOutput.shp")

layer3 = iface.addVectorLayer("/Users/casta2k/ShapefilePractice/intersect/intersectOutput.shp", "intersect", "ogr")

layer = iface.activeLayer()
provider = layer.dataProvider()

areas = [ feat.geometry().area()
      for feat in layer.getFeatures() ]

field = QgsField("area", QVariant.Double)
provider.addAttributes([field])
layer.updateFields()

idx = layer.fieldNameIndex('area')

for area in areas:
    new_values = {idx : float(area)}
    provider.changeAttributeValues({areas.index(area):new_values})

I'd like to be able to run the code directly from a python script, without having to open the GUI. I've tried searching the internet, but I've found nothing helpful. Do I have to start by running QgsApplication.initQgis()? I thought I had to, but I get a segmentation fault. Also, I tried changing iface.addVectorLayer() to QgsVectorLayer(), which seems to work fine in the Python console, but it seems like the files aren't loaded when I run the python script directly.

I tried to include what @Paulo suggested, but I'm still getting errors. The code now starts as follows:

import sys
import os

sys.path.append("/Applications/QGIS.app/Contents/Resources/python")

from qgis.core import *
import qgis.utils
from qgis.analysis import *

# supply path to where is your qgis installed
QgsApplication.setPrefixPath("/Applications/QGIS.app/Contents/MacOS", True)

# load providers
QgsApplication.initQgis()

But on the last line (QgsApplication.initQgis()) I'm getting the following error:

QCoreApplication::applicationDirPath: Please instantiate the QApplication object first
QCoreApplication::applicationDirPath: Please instantiate the QApplication object first
Segmentation fault: 11


texts: QGIS 3.4 many-to-many relations in forms

texts:
On QGIS Attribute Forms, I want to display related Features which are connected via a many-to-many relationship table:

Tables / QGIS Layers:

Feature a: a_id, geometry, some more attributes

Feature b: b_id, geometry, some more attributes

Join table: join_a_to_b: a_id, b_id

Now I define one-to-many-Relations via Project Properties:

Relation 1: Parent: Feature a, Child: Join Table referencing a_id

Relation 2: Parent: Feature b, Child: Join Table, referencing b_id

This way I can show all connected join table entries in the attribute form, but I do not manage to show connected Features of type "b" directly in an "a" attribute Form.

From the docs, this should be possible this way - any hints what I'm doing wrong? Swapping Parent and Child wont work either.

Using QGIS 3.4.3 on Windows 7 (same behaviour with Linux) on top of PostgreSQL


texts:
I want to perform an analysis on urban areas in the UK. As, indeed, urban areas can be defined in multiple ways, I would like to rely on some well-accepted definitions and polygons.

I found this EU project that has nice data about functional urban areas in the Union, but the source data is rather old (2000).

Then I found UK city regions, one of the many admin units in the UK, but I can't find shapefiles for them (and they look very large).

Moreover, I would like to find a point dataset with cities and towns (a gazetteer) in the UK, with population information in order to use them as labels dynamically (e.g., show cities > 1M, etc.). The OS Open Names product is way too detailed, and does not seem to contain a handy categorisation between small and large centres.

Any tips?


texts: Why does GeoServer (UI) not refresh the layer after a succesfull REST PUT?

texts:
I have a layer configuration in Geoserver (space1:my_layer) which I want to change.
I change the layer via REST PUT on http://localhost/rest/workspaces/space1/datastores/store1/featuretypes/my_layer
Afterwards I do a GET on the same URL and the configuration is changed. Perfect.
Then I launch the UI / web of GeoServer and then still my old configuration is shown there. No matter if I (hard) refresh the page, it doesn't update.

Is there a commit needed? or a refresh? is the REST response more right than the web response? Why should I do?


texts:
So I have a .net core app which currently throws an error in console saying:

Access to XMLHttpRequest at 'https://app.test.mk/arcgis/rest/services/karpos/Karpos_Vector_Izvodi/MapServer?f=json&dpi=96&transparent=true&format=png8' from origin 'https://localhost:44348' has been blocked by CORS policy: The 'Access-Control-Allow-Origin' header contains multiple values 'https://localhost:44348, *', but only one is allowed.

and this is what my Startup file looks like:
public void ConfigureServices(IServiceCollection services)
        {
            services.AddCors(o => o.AddDefaultPolicy(builder =>
            {
                builder
                    .WithOrigins("https://localhost:44348/",
                                 "Access-Control-Allow-Origin", "Access-Control-Allow-Credentials")
                    .AllowAnyOrigin()
                    .AllowAnyMethod()
                    .AllowAnyHeader();
            }));
            services.AddControllersWithViews();
            services.AddRazorPages().AddNewtonsoftJson();
        }}

public void Configure(IApplicationBuilder app, IWebHostEnvironment env)
        {
            if (env.IsDevelopment())
            {
                app.UseDeveloperExceptionPage();
            }
            else
            {
                app.UseExceptionHandler("/Home/Error");
                // The default HSTS value is 30 days. You may want to change this for production scenarios, see https://aka.ms/aspnetcore-hsts.
                app.UseHsts();
            }
            app.UseCors();

            app.UseHttpsRedirection();

            app.UseStaticFiles();

            app.UseRouting();

            app.UseAuthorization();

            app.UseEndpoints(endpoints =>
            {
                endpoints.MapControllerRoute(
                    name: "default",
                    pattern: "{controller=Home}/{action=Index}/{id?}");
            });

        }

and this is at the top of my View, in a script tag.
require(["esri/config"], function (esriConfig) {

        esriConfig.defaults.io.corsEnabledServers.push("http://app.test.mk");
        esriConfig.defaults.io.corsEnabledServers.push("https://localhost:44348/");
       // esriConfig.defaults.io.corsEnabledServers.push("http://app.test.mk/arcgis/rest/services/karpos/Karpos_Vector_Izvodi/MapServer");
       // esriConfig.defaults.io.corsEnabledServers.push("https://app.test.mk/arcgis/rest/services/karpos/Karpos_Vector_Izvodi/MapServer?f=json&dpi=96&transparent=true&format=png8");

    });

Anyone has any idea what the problem could be? I've been trying things for two days and the app keeps throwing the same error.


texts: Determining intersection between GeoJSON feature and given coordinates?

texts:
I'm overlaying a Leaflet map with features of GeoJSON, fetched with D3 and added as SVG elements to the DOM.

My GeoJSON contains too many features to render at once, so my goal is to add to the map only features, that intersect current cursor location, so I can have a smooth hover effect.

I can get current map coordinates of the cursor from a mousemove event, but I have to determine if current coordinates intersect a feature in GeoJSON and exactly what feature is intersected.

Thinking further, maybe I have to make a small buffer, so it's easier for the user to hit the closest feature.

Is there a fast way in Leaflet or D3 (or even pure JavaScript) to get the intersected feature from GeoJSON feature collection?


texts: Serving OSM tiles - mod_tile/renderd/postgresql problem: I sometimes get 404 Not Found

texts:
When I request tiles I get sometimes 404 Not Found. Usually, if I request the same tile again it appears alright (with 200 OK).

This is what I found in renderd logs:

renderd[734]: An error occurred while loading the map layer 'default': Postgis Plugin: ERROR: could not open file "base/16386/240582": Interrupted system call
in executeQuery Full sql was: 'SELECT * FROM ( SELECT way, way_area AS area, COALESCE(landuse, leisure, "natural", highway, amenity, tourism) AS type
FROM planet_osm_polygon
WHERE way_area > 100000
ORDER BY way_area DESC
) AS data LIMIT 0'

This is my 000-default.conf:

<VirtualHost *:80>
        # The ServerName directive sets the request scheme, hostname and port that
        # the server uses to identify itself. This is used when creating
        # redirection URLs. In the context of virtual hosts, the ServerName
        # specifies what hostname must appear in the request's Host: header to
        # match this virtual host. For the default virtual host (this file) this
        # value is not decisive as it is used as a last resort host regardless.
        # However, you must set it for any further virtual host explicitly.
        #ServerName www.example.com

        LoadTileConfigFile /usr/local/etc/renderd.conf
        ModTileRenderdSocketName /var/run/renderd/renderd.sock
        # Timeout before giving up for a tile to be rendered
        ModTileRequestTimeout 0
        # Timeout before giving up for a tile to be rendered that is otherwise missing
        ModTileMissingRequestTimeout 600

        ModTileCacheDurationMax 31536000
        ModTileCacheDurationMinimum 31536000
        ModTileMaxLoadMissing 999999999
        ModTileEnableTileThrottling Off
        ModTileBulkMode Off
        TimeOut 600

        ServerAdmin webmaster@localhost
        DocumentRoot /var/www/html

        # Available loglevels: trace8, ..., trace1, debug, info, notice, warn,
        # error, crit, alert, emerg.
        # It is also possible to configure the loglevel for particular
        # modules, e.g.
        #LogLevel info ssl:warn

        ErrorLog ${APACHE_LOG_DIR}/error.log
        CustomLog ${APACHE_LOG_DIR}/access.log combined

        # For most configuration files from conf-available/, which are
        # enabled or disabled at a global level, it is possible to
        # include a line for only one particular virtual host. For example the
        # following line enables the CGI configuration for this host only
        # after it has been globally disabled with "a2disconf".
        #Include conf-available/serve-cgi-bin.conf
        Alias /log/ "/var/log/"
    <Directory "/var/log/">
           Options Indexes MultiViews FollowSymLinks
           AllowOverride None
           Order deny,allow
           Deny from all
           Allow from all
            Require all granted
    </Directory>
</VirtualHost>

My postgresql.conf:
https://pastebin.com/veYxsui5

I'm running the setup in Docker container on c5.xlarge instance on AWS.

How to eliminate 404s?


texts: Extract raster values around point features with kernel

texts:
I have a point layer of known landslide locations, and a raster layer (30m) of vegetation cover (categorical). My goal is to extract the prevailing raster value around the point locations.

Almost all landslides are located within 30m of a road which means that if I use the "extract value to point' tool, my out put would only correspond to the raster class "Developed-Roads". This is not helpful since I need to look at cells further away from the road in order to obtain "true" vegetation values like "Tree Cover >= 70 and < 80%", for instance. I can't create buffers around the point locations and run zonal statistics, because there is too much overlap.

When I researched this problem, I also came across the "focal statistics" tool, but I don't see how that is helpful. Even though the tool allows you to define various neighborhoods that would certainly be an advantage, it doesn't allow me to input a point feature and a raster layer at the same time.

Can I use focal statistics, if so how? If not, what other options are there?

Check out the picture to get a better idea of the problem:


texts: Randomly altering raster map of habitat types?

texts:
I have a raster of habitat types for a specific area in Scotland. I need to create future habitat scenarios with changes in habitat to assess the population viability of a bird species.

For example, in the future there might be 10% more forestry in the area. I would like to alter the current map by randomly adding the forestry in blocks of a certain size. I am, so far, thinking along the lines of selecting random points from a raster which identifies areas where the forestry could occur and growing the correct sized blocks using some sort of cellular automata.

Does this seem like the best way of going about this? Is there a better method?

If this is the best way available, how could I do this in, preferably, R?
(I am currently looking at the rpoints function in "spatstat" along with the CellularAutomata package)

I also have access to GRASS, QGis and ArcMap 10 if there are simpler ways in any of them.


texts: Get field value: JSON response from feature service

texts:
I have JavaScript in work order management software:

load("nashorn:mozilla_compat.js");
importPackage(java.util)
importPackage(Packages.psdi.server)
var ctx = new HashMap();
ctx.put("url","https://something.com/arcgis/rest/services/something/zones/MapServer/8/query?geometry=123456%2C1234567&geometryType=esriGeometryPoint&spatialRel=esriSpatialRelIntersects&outFields=*&returnGeometry=false&f=pjson");
service.invokeScript("LIB_HTTPCLIENT",ctx);
var jsonResp = ctx.get("response");
service.log(jsonResp);

The script sends a request (with coordinates) to a feature service:

  https://something.com/arcgis/rest/services/something/zones/MapServer/8/query?geometry=123456%2C1234567&geometryType=esriGeometryPoint&spatialRel=esriSpatialRelIntersects&outFields=*&returnGeometry=false&f=pjson

The feature service (REST) returns information about the feature that the coordinates/point intersected with.

{
 "displayFieldName": "DAY",
 "fieldAliases": {
  "OBJECTID": "OBJECTID",
  "DAY": "DAY",
  "ZONE": "Area",
  "SHAPE.AREA": "SHAPE.AREA",
  "SHAPE.LEN": "SHAPE.LEN"
 },
 "fields": [
  {
   "name": "OBJECTID",
   "type": "esriFieldTypeOID",
   "alias": "OBJECTID"
  },
  {
   "name": "DAY",
   "type": "esriFieldTypeString",
   "alias": "DAY",
   "length": 10
  },
  {
   "name": "ZONE",
   "type": "esriFieldTypeString",
   "alias": "Area",
   "length": 6
  },
  {
   "name": "SHAPE.AREA",
   "type": "esriFieldTypeDouble",
   "alias": "SHAPE.AREA"
  },
  {
   "name": "SHAPE.LEN",
   "type": "esriFieldTypeDouble",
   "alias": "SHAPE.LEN"
  }
 ],
 "features": [
  {
   "attributes": {
    "OBJECTID": 327,
    "DAY": "Tuesday",
    "ZONE": "1",                     <-------zone is here
    "SHAPE.AREA": 5154109.08323004,
    "SHAPE.LEN": 10563.5744286387
   }
  }
 ]
}

I want to isolate the zone number from the JSON text so that I can insert it into a record/field.

How can I parse/isolate the zone number using JavaScript?


texts: How can I refer to two tables in ArcPy's CalculateField_management?

texts:
I tried to refer to the different table's fields just like I saw in How to Calculate Field in joined table using ArcPy? but it doesn't work.

There are two tables in my example. I have created table "c":

I generated a summarized table (called ttt) with right click on the field "helysz" and summarize. Then I added a field "hely", and copied the newly created "helysz" field's content into that, just for sure, for different field names. Then I added the ossz field to the second table "ttt". After this I linked the c's helysz field to the secondly created table's hely field. The table c's fields from the second OBJECTID field are from the linked table. They have appeared only after the linking. (ArcGIS shows the first one in this view).

This is the second table "ttt":

I wrote this ArcPy to go through the hely field with CalculateField_management and find those certain values from the original helysz field, then work with the values what are in one record with the found helysz.

ArcPy Code

The error message I get:

  raise e
  arcgisscripting.ExecuteError: ERROR 000539: Invalid field D:/JAN17ALG/joo/Felsotarkany_energetika_kiserleti.gdb/c.helysz
  Failed to execute (CalculateField).

hely and helysz means the topographical number. And one topographical number can have more buildings.
the data of the buildings would be the "c" table.
the data of the topographical numbers (or fields what have multiple buildings) are in the "ttt" table.

Why this doesn't work?


texts: Breaking dynamic text line in Data Driven Pages

texts:
I am using data driven pages and some of my dynamic text attributes are too long for my layout.

Is there a way to break the line into 2 or more lines with out editing the data?

Example:

FROM

Billy Joes Dry Cleaners LLC

TO

Billy Joes
Dry Cleaners LLC


texts: How to implement a "Snap to Road"-functionality?

texts:
I am having a problem with a little application I'm trying to develop. What I do is take a phone with me and drive around the city on a bike.

I collect all the Latitude/Longitude-pairs and at home I plot them on a map.

What I get, however, is zig-zag lines, because the GPS location is not perfectly accurate. Here's an example plot:

Obviously I was driving in a more or less straight line and not through the hotel. ;)

This is perhaps basic for GIS-experts, but here's my question: How do I get better results, with the red line (and thus the recorded points) being on the roads?

These are the things I've considered:

Rolling-average over measurements. E.g. take 5 coordinates and calculate average coordinate. Tried it and it doesn't work. It didn't work with the data shown in the picture. (But the line was smoother..)
Only accepting measurements where the device (iPhone) reports a high accurracy. Tried this and it helps a little bit, but still not good enough.
Correlating of GPS data with acceleration and direction sensors of the phone. I haven't tried this yet and it sounds non-trivial.
Using some sort of "Snap to Road"-algorithm.

I think it should be possible to get data on roads for a given location. (OSM?) Then perhaps the nearest road and the nearest point on that road from the current position could be calculated and the position could be adjusted to be ON the road. Sort of like the "snapping to grid" function of graphical interfaces.

If this is the right approach, what are some pointers for software and methods I can use?

I assume I will have to deal with the following problems:

Huge dataset for any given country or city (I'd have to slice tiny portions out, but also store the whole dataset on my own server?)
Transformation of formats..
High CPU/memory-usage(?)
A complex algorithm(?)

Looking at car navigation software it seems that this problem has long been solved. Are there any implementations that can be re-used, preferably of course for the iOS?

Thank you

UPDATE:

I have implemented a simple map matching algorithm:

I downloaded several hundred MB of road data from OpenStreetMap (OSM) and inserted them into my DB.
For each measured GPS-point I fetched all OSM-nodes in a certain radius.
For each fetched node I fetched all roads (series of nodes) from DB and treated them like line segments in 2d-vectorspace.
I calculated the closest point (distance from gps-measurement) for each segment and selected the one with the shortest distance.
I've changed the lat-lon coordinates of the measurement to those of the matched point.

Here is the actual result for the example plot above:

It's just slightly better.

I think the next step is to calculate multiple points in step 4 and use each of them in a probability / plausibility step 4a, where each one is checked for things like wether the bike really took a turn onto another street.

It's not real-time so it should be possible to take the future into account (e.g. to look at the next measurement that followed, or the next few.), and to calculate a probability score for each candidate point.

Anyway, just wanted to post a quick update. Perhaps this will help someone.


texts: PostGIS raster calculate average slope of 3x3 neighborhood around a geometry point

texts:
I have two tables

corners_final (vector point data)
dem_slope (1ft slope DEM, calculated in degrees)

With this simple query I am able to extract the slope value that intersects a given point

SELECT ST_Value(rast,geom)
FROM dem_slope
JOIN corners_final ON ST_Intersects(geom,rast)
WHERE cornerid=1160060

But I want to find the average slope around that point by 3x3 area.

I cannot figure out how to use ST_Neighborhood... Should I be using ST_MapAlgebra?

I do not see a lot of examples on these functions.


texts: How to build a map that highlights only one state?

texts:
I'm a journalist working on maps of several Alaska databases. I'm learning as I go and have a question that's probably very basic, but very important for my projects:

Using CartoDB, how can I create a map that displays only Alaska? (ie not the rest of the US or other countries)

A simple step-by-step answer to this would be immensely helpful!


texts: How to add a buffer on a point in OpenLayers 3

texts:
I am trying to add a buffer around some points shown on a map. I have tried to do something similar back in the days of OL2 and I remember it was not straight forward.

Is there a simple way to do this in OL3? I have found this:

ol.extent.buffer(extent, value, opt_extent)

but I am not able to find an example.


texts: How to Edit Domain and Have Edit Show in Attributes Dropdown

texts:
I have created a previous domain with the same user login in ArcGIS 10.3. I wanted to add a coded value to this domain and, have the new addition show up on the drop-down options in the attributes when I add a new feature. I am able to add it to the domain. The issue is the addition is not showing up in the attributes when the new feature is created.

I have tried just plain open the properties in the gdb and add the addition to the domain. I have tried to use the tool "Add Coded Value To Domain". The addition is not showing up in the attrubutes.  I have tried to convert the domain to a table and than table to domain, new domain name. This worked but, when I changed the field domain to this new domain, no drop-down is shown when editing this field for a new feature.


texts: Explaining ${z}/${x}/${y}.png for offline map?

texts:
Im trying to build map with OMS under Linux (Ubuntu) to made an offline apps. Im working with Ionic and building just for Android for write now.
I have the base code:

 var map;
        function init(){
            map = new OpenLayers.Map('map',
                    { maxExtent: new OpenLayers.Bounds(-20037508.34,-20037508.34,20037508.34,20037508.34),
                      numZoomLevels: 19,
                      maxResolution: 156543.0399,
                      units: 'm',
                      projection: new OpenLayers.Projection("EPSG:900913"),
                      displayProjection: new OpenLayers.Projection("EPSG:4326")
                    });
            var newLayer = new OpenLayers.Layer.OSM("New Layer", "/tiles/${z}/${x}/${y}.png", {numZoomLevels: 19});
            map.addLayer(newLayer);
            map.addControl(new OpenLayers.Control.LayerSwitcher());
            var lonLat = new OpenLayers.LonLat(1.145, 51.865).transform(map.displayProjection,  map.projection);
            if (!map.getCenter()) map.setCenter (lonLat, 11);
            }

I got my tiles made with Maperative here:

/home/in/Android/offline/tiles/12

Inside 12 folder there are three more folders, where the png files are they are called 2223, 2224, 2225.

How does the argument "/tiles/${z}/${x}/${y}.png" really work? I have search a lot but I havens see a clear explanation.
Shall I need to put the files in a special folder?
If I use a complete path (Windows example):

  var newLayer = new OpenLayers.Layer.OSM("New Layer", "file://localhost/C:/Users/Ed/Downloads/Maperitive/Maperitive-1141/Tiles/${z}/${x}/${y}.png"

when I turn build to an Android apps will break the link.

Is a relative path a better option?


texts: Stop Qgis layer from refreshing every time I pan

texts:
I've got a layer in a QGIS project which source is a PostGIS view. It's a quite complex view and needs like 20 seconds to load (it makes several spatial analysis).

My problem is that the QGIS project becomes unmanageable because every time I pan it loads the layer again and gets stuck another 20 seconds..

Is there any option to disable for this QGIS project, or even better, for this layer to stop requerying data every time I pan the map?


texts: Where to find Unified Soil Classification System .Style?

texts:
I'm looking for the Unified Soil Classification .style patterns to add to my style manager. I haven't had much luck with a variety of web searches for either ArcMap and AutoCAD.
Does anyone know where I might be able to find them?


texts: Map tiles not loading until a zoom is executed

texts:
The initial map is missing tile. No errors in the console. If I zoom though, the tiles refresh and all of them load. See attached screenshots.

I have used //js.arcgis.com/3.18 API.

Please let me know how to proceed on this.


texts: "Out of Bounds" CAD-originated shapefiles unable to georeference

texts:
I'm trying to import several polyline files from an AutoCad drawing into ArcMap. I keep running into an "out of bounds" error when trying to georeference the CAD layers to where they should be.

I don't believe the CAD files have a coordinate system. The shapefiles are located in the middle of the wilderness some 300k meters away from the expected location. When trying to georeference them, selecting Adjust does nothing and the Link Table shows 0 for the X and Y regardless of where I place the links.

I've tried using the Copy Features to Feature Class tool with environments set to use Union of Inputs for the extent and disable M and Z values. The resulting line files are also unable to be georeferenced.

When I try to edit, select all, and then Move using a set X,Y point, I'm told the file can't be moved because it is out of bounds.

I've tried creating new GDB's with a set coordinate system..but that didn't work either.

What could I try next?

I'm trying to place the CAD files into NAD_1983 State Plane MO East in Feet. I am not sure of the XY coordinates of the CAD file.

Edit: The original CAD file has no coordinate system, per the surveying company. I have a folder of exported shapefiles from the CAD file that I cannot georeference because any resulting feature classes made from those shapefiles are "out of bounds."


texts: Is there anything like http://jsonlint.com or http://geojsonlint.com for WKT? How can I validate a complex WKT string, say e.g. a Multipolygon?


texts: How to extract value from a raster based on lat and lon?

texts:
I have several netcdf files that I can read using the netCDF4 library in python. The data is gridded in 4320 rows and 2160 columns.

I also have a list of lat/lon (in decimal degrees) and I am trying to find the value of the gridcell in which the coordinates lie.

Currently, I am converting the coordinates into gridcell boundaries by a simple linear transformation and querying the netcdf file using these new points. Is this approach valid?

# prelim
from netCDF4 import Dataset
import numpy as np

# coordinate data
lons = [40.822651, 40.69685, 40.750038, 40.798931, 40.74612]
lats = [19.83832, 20.007561, 19.97426, 19.86334, 19.843889]

# find boundaries of gridcells corresponding to lat/lon
lons_grid = np.ceil(12*(lons + 180))
lats_grid = np.ceil(12*(lats + 90))

# find value in gridcell for given lat/lon
fnc = Dataset(filename.nc, 'r')
lat = fnc.variables['latitude'][:]
lon = fnc.variables['longitude'][:]
variable = fnc.variables['variable'][:]
print variable[lons_grid, lats_grid]

If not, is there a way to intersect point data with a raster in python without using arcpy? (My netcdf files are very large and cannot be read with arcpy.)


texts: Panning & Zooming do not work when a vector layer feature fills the entire map

texts:
I can't seem to find a solution to allow panning or zooming when a vector layer polygon fills the entire map area. My current solution is to use the zoom out button to zoom out far enough to find an open area of the map to click on to pan and zoom. Is there a way to propagate the pan/zoom capabilities through the vector layer features?

Even in the OL code (https://github.com/openlayers/openlayers/blob/v6.5.0/src/ol/control/Zoom.js), I don't see any reference to vector layers preventing zooming:
zoomByDelta_(delta) {
    const map = this.getMap();
    const view = map.getView();
    if (!view) {
      // the map does not have a view, so we can't act
      // upon it
      return;
    }
    const currentZoom = view.getZoom();
    if (currentZoom !== undefined) {
      const newZoom = view.getConstrainedZoom(currentZoom + delta);
      if (this.duration_ > 0) {
        if (view.getAnimating()) {
          view.cancelAnimations();
        }
        view.animate({
          zoom: newZoom,
          duration: this.duration_,
          easing: easeOut,
        });
      } else {
        view.setZoom(newZoom);
      }
    }
  }
}


texts: Add multiple GDB files in subdirectories to mxd using python

texts:
I have GDB files located in many sub-Folders and sub-sub-Folders. All Sub Folders are located in one large directory
i asked it for layer file in
Add multiple Layer files in subdirectories to mxd using python
I use this code:

import arcpy
import os

workspace = r"C:\Project\layers"
mxd = arcpy.mapping.MapDocument(r"C:\Project\Project.mxd")
df = arcpy.mapping.ListDataFrames(mxd, "*")[0]
walk = arcpy.da.Walk(workspace, datatype="Layer")
for dirpath, dirnames, filenames in walk:
    for filename in filenames:
        arcpy.mapping.AddLayer(df, arcpy.mapping.Layer(os.path.join(dirpath, filename)))
        print filename
    mxd.save()
del mxd

How can i do it for GDB files ?


texts: Search only all the nodes that have tags

texts:
I'm trying to execute queries in Overpass Turbo looking for only nodes that contain a tag definition,

for example i'm running a query like this:

<osm-script output="json">
 <query type="node">
  <bbox-query {{bbox}}/>
 </query>
  <print mode="body"/>
  <recurse type="down"/>
  <print mode="skeleton"/>
</osm-script>

this return all nodes with and without tags,

how to query only nodes with tags?


texts:
ST_Overlaps, ST_Contains, ST_Touches, ST_Within works for TIN with point data. But, when I try with 2 TIN it gives error. Could you please help me understand the reason for failure and how to fix it. I need comparisons only between TIN objects.

SELECT ST_Overlaps
('TIN Z (((0 0 0,0 0 0.5,0 0.5 0.5,0 0 0)),((0 0.5 0,0 0 0,0 0.5 0.5,0 0.5 0)))'::geometry,
 'POINT(0.5 0 0.5)'::geometry);

SELECT ST_Overlaps
('TIN Z (((0 0 0,0 0 0.5,0 0.5 0.5,0 0 0)),((0 0.5 0,0 0 0,0 0.5 0.5,0 0.5 0)))'::geometry,
 'TIN Z (((0 0 0,0 0 0.5,0 0.5 0.5,0 0 0)),((0 0.5 0,1 1 1,1 0.5 0.5,0 0.5 0)))'::geometry);

  ERROR:  Unknown geometry type: 15 - Tin
  CONTEXT:  SQL function
  "st_overlaps" statement 1

-Note: ST_3dIntersects works with 2 TIN objects