texts: WFS-T request works in localserver but fails in remote server

texts:
I have a local server with a geoserver 2.7. I use python requests in order to add geometry in a postgis store. It works great. This is how it looks like:

headers = {'Content-Type': 'application/xml'}
url = 'http://localhost:8080/geoserver/gn/wfs'
xmlstr = """<wfs:Transaction service="WFS" version="1.1.0"
  xmlns:wfs="http://www.opengis.net/wfs"
  xmlns:gml="http://www.opengis.net/gml">
  <wfs:Insert>
   <dkar>
    <the_geom>
     <Point srsName="http://www.opengis.net/gml/srs/epsg.xml#4326">
      <pos>16.259765625 38.16911413556086</pos>
     </Point>
    </the_geom>
   </dkar>
 </wfs:Insert>

"""

status_code = requests.post(url, data=xmlstr, headers=headers, auth=('user', 'pass')).status_code

This worked fine but when I try to execute the same code on the server (with geoserver 2.5) I get an error:

No such feature type http://www.opengis.net/wfs:dkar

I did a DescribeFeature and I got:

<?xml version="1.0" encoding="UTF-8"?><xsd:schema    xmlns:xsd="http://www.w3.org/2001/XMLSchema" xmlns:geonode="http://geonode.org/" xmlns:gml="http://www.opengis.net/gml" elementFormDefault="qualified" targetNamespace="http://geonode.org/">
 <xsd:import namespace="http://www.opengis.net/gml" schemaLocation="http://staging.geonode.wfp.org:80/geoserver/schemas/gml/2.1.2/feature.xsd"/>
 <xsd:complexType name="dkarType">
 <xsd:complexContent>
  <xsd:extension base="gml:AbstractFeatureType">
    <xsd:sequence>
      <xsd:element maxOccurs="1" minOccurs="0" name="the_geom" nillable="true" type="gml:PointPropertyType"/>
      <xsd:element maxOccurs="1" minOccurs="0" name="countries" nillable="true" type="xsd:string"/>
      <xsd:element maxOccurs="1" minOccurs="0" name="numbers" nillable="true" type="xsd:int"/>
    </xsd:sequence>
  </xsd:extension>
 </xsd:complexContent>
</xsd:complexType>
<xsd:element name="dkar" substitutionGroup="gml:_Feature" type="geonode:dkarType"/>
</xsd:schema>

Any idea what can cause this error?


texts:
[Using QGIS 2.12 Lyon on Mac with OS X 10.10.5]

I am struggling to georeference some old hydrographic survey maps that are in the form of TIFF files (.tif). If I open the georeferencer tool and set ground control points (GCPs) from the map canvas, then run the georeferencer, it offers to save the GCPs in a file. I have done so every time.

Now, if I change the project CRS, is the file of GCPs not going to work? That is, does the saved file of GCPs only work in the original project CRS? I ask because I have changed the project CRS and then loaded the old GCP file into the georeferencer tool, and in the image it appears that everything is fine. But georeferencing still doesn't work.


texts: Using UpdateCursor and Where Clause to manipulate attributes?

texts:
I am working on a practice problem, and i have hit a wall. Here is what the question is asking.

  Write a script that updates the attribute of SIZE with the values
  ‘Large’ or ‘Small’, depending on the value of the TOT_ENP (the number
  of passengers) field. If the number of passengers is larger than
  100,000, then it’s a large airport (SIZE = ‘Large’); otherwise, it’s a
  small airport (SIZE = ‘Small’).

Heres what I have thus far:

import arcpy

arcpy.env.workspace = "E:/lecture08data2"

fc = "airports.shp"
cursor = arcpy.da.UpdateCursor(fc, ["SIZE"], '"SIZE" = \'Large\' and TOT_ENP > 100000 or "SIZE" = \'Small\' and TOT_ENP < 100000')
for row in cursor:
    cursor.updateRow()

del row
del cursor


texts: How to use OSM Bright with Mapnik Renderd

texts:
I have set up an OSM tile server using instructions at this blog -
http://seshagiriprabhu.wordpress.com/2013/07/21/building-an-openstreetmap-tile-server-on-ubuntu-12-04-lts/

The slippymap renders fine and work nice.
Now, I want to change the style of the map to OSM Bright. I do not want to (rather, do not have technical expertise to) make any changes to the template.
I wish to learn how to render my maps usng OSM Bright or any other Carto CSS templates.
I have searched around the net for quite a few days, but cannot find any simple step by step instructions on how to make my mapnik Renderd to display maps with OSM Bright template.

My instance runs on Ubuntu 12.04. Please give me a step by step guide on how to make this happen.


texts: How to force QGIS Centroid Fill inside the Polygon?

texts:
I would like to use a centroid fill to enable me to place a letter inside a polygon. The advantage of using a centroid fill is that it also shows up in a legend, whereas labelling the layer would mean I have to put a text label on the legend in the composer.

The problem I have with the centroid fill is that the centroid can fall outside of the polygon, and as such my letter will not be within the polygon I need it to be.

I have included an example that shows where the centroid fills have worked, and one where it hasn't (the label shown in red):

Is there a way to force the centroid fill to be inside of the polygon?

Or is there another way to achieve this?


texts: Create output folder for merge result of Model Builder on QGIS 3.10

texts:
I have export my model builder to python's script and now i would like to save the result to an created folder. But the folder name is based on a layer entity. Example: I have 3 layers (city, store, population) as input which gives me 2 results. I would like to save these 2 shapes in 1 folder which will have the name of a city. If I take Washington, Boston and Detroit, and launch my algorithm, it will take me out (thanks to the iteration), 6 layers in all. So I would like my 2 layers created for Washington to be automatically placed in the Washington folder and the same for the other 2 cities.

I tried that but it doesn't work : Create a folder adding a record
I have an error of the type: iface.activeLayer(city): too many argument.
I don't understand where is the problem.

EDIT : Now I can get the algorithm to work. This one creates the files of the cities concerned (so 3 files). But I can't put the files in the folders that concern them ...

My code

from qgis.utils import iface
import os

france_substation_source_rte = iface.activeLayer()
selected_features = france_substation_source_rte.selectedFeatures()
idx= france_substation_source_rte.fields().indexFromName('code_poste')
path = r"C:\Users\utilisateur\Documents\Maps_GIS\dataset\data_analysis\automation_results\_"
for feat in selected_features:
    attrs = feat.attributes()
    if not os.path.exists(path + attrs[idx]):
        os.makedirs(path + attrs[idx])


texts: Forcing gdal_retile to produce full tiles when original file's size is incompatible?

texts:
I have a set of large geotiff (DEM + orthoimagery) that I need to generate image pyramids from. So far gdal_retile from osgeo4w is working very well for my use case with one exception:  I need all my tiles to be of uniform pixel size (eg: 1025x1025). gdal_retile will, as its default behavior, produce smaller tiles at the edges if the original file's pixel dimensions are not a multiple of the tile size.

Alternatively, is there a gdal command that will let me pad the DEM and orthoimagery with 'empty' / 'transparent' data out to an even multiple of my tile size)?


texts: ArcMap Field Calculator Python string concatenation malfunction

texts:
I am trying to concatenate three fields into one using Python with the rjust padding function.

My function in the code box is the following:

def AutoSchluessel(gmk, zae, nen):
       return "{}0{}{}".format(gmk, str(zae).rjust(7,"0"), str(nen).rjust(5,"0"))

And the Expression is:

flcode = AutoSchluessel( !gmkgcode!, !zaehler!, !nenner!)

What i am getting is for the fields is the following result:

  gmkgcode = 5902;
  zaehler = 32;
  nenner = 17

  flcode = 5902.0000032.0017.0

It should in my opinion though result in something without points (59020000003200017).
What am I missing?


texts: Change the colors of polygons on a choropleth map with leaflet

texts:
I have data from different countries (=polygons) from 1990 to 2012 and depending on the data value I want to map the changes with a timeline enabling the user to visualize the changes (with colors) from 1990 all the way to 2012.

I am using leaflet for the first time, so a little more detailed answer would be great.

I got the shapefile to show on the map, but when I try to the style function it doesn't work (no colors). the column "greenhousegaz_1990" contains the data I want to show.

I tried with the ID column just to test and it worked fine, colors did show correctly and in range.

function getColor(d) {
    return
           d > 900000  ? '#762A83' :
           d > 750000  ? '#af8dc3' :
           d > 600000  ? '#e7d4e8' :
           d > 450000  ? '#f7f7f7' :
           d > 300000  ? '#d9f0d3' :
           d > 150000  ? '#7fbf7b' :
                         '#1B7837';
}

function style(feature) {
    return {
        fillColor: getColor(feature.properties.greenhousegaz_1990),
        weight: 2,
        opacity: 1,
        color: 'white',
        dashArray: '3',
        fillOpacity: 0.7
    };
}


texts: Creating line pattern fill using geometry generator (PyQGIS)

texts:
I'm trying to make a line pattern fill in QGIS using the geometry generator and a Python function.

The code is inspired from this answer : Save line geometry as points geometry in QGIS (pyqgis)

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

@qgsfunction(args='auto', group='Custom')
def fill_line_pattern(x, feature, parent):

#create a base line that is 45 degrees
base_xmin = 2285865.0
base_xmax = 2669093.0
base_ymin = 1017594.0
base_ymax = 1399851.0

#calculating interval between parrallels
countX = int((2800000.0-2285865.0) / xInterval)
#empty list
pt=[]
pts = []
intersected_lines = []
#raw lines creation
for xOff in range(countX+1):
    ptXmin = base_xmin + xOff*(xInterval)
    ptXmax = base_xmax + xOff*(xInterval)
    ptmin = QgsPoint(ptXmin,base_ymin)
    ptmax = QgsPoint(ptXmax,base_ymax)
    pt = ([ptmin, ptmax])
    pts.append(pt)
#intersect lines with polygons
#for f in feat :
for p in pts :
    rawlines = QgsGeometry.fromPolyline(p)
    if feature.geometry().intersects(rawlines):
        geom = feature.geometry().intersection(rawlines)
        intersected_lines.append(geom.asPolyline())
#create lines
for line_ok in intersected_lines:
    return QgsGeometry.fromPolylineXY(line_ok)

The issue is that the first line only is crossing the polygon...

How can I have all the lines crossing the polygons?

In comparison to the inspirational code above, there's no loop for the polygon features.

In the function, the feature called is a QgsFeature which can't use getQgsFeatures().

If a for loop is used :

for f in feature :
    for p in pts :
    rawlines = QgsGeometry.fromPolyline(p)
        if f.geometry().intersects(rawlines):
            geom = f.geometry().intersection(rawlines)
            intersected_lines.append(geom.asPolyline())

This error message appears:

  Evaluation error: 'int' object has no attribute 'geometry'

EDIT

The code below does almost what I want, it seems a little bit buggy because not all the intersections appear... The other disadvantage is the long time QGIS takes to render it on the canvas and in the composer.

On the other hand, the pdf output size is so light and is competitive with other software.

Maps are made to be shared: QGIS biggest issue is PDF output file size

The inspiration of the source : https://www.lutraconsulting.co.uk/blog/2015/06/05/qgis-function-editor/

Lutra Consulting uses a cache but I didn't find the way to do it...

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

def _layer_intersection(layer, xInterval):
#make 45 degrees lines from layer extent
    ext = layer.extent()
    xmax = ext.xMaximum()
    ymax = ext.yMaximum()
    xmin = ext.xMinimum()
    ymin = ext.yMinimum()

    pts_hash_box = []
    line_in_polygon = []

    x_base_min = xmin - (ymax-ymin)
    countX = int(((xmax-xmin)+(ymax-ymin)) / xInterval)

    for xOff in range(countX+1):
        ptXmin = x_base_min + xOff*(xInterval)
        ptXmax = xmin + xOff*(xInterval)

        ptmin_hash_box = QgsPointXY(ptXmin,ymin)
        ptmax_hash_box = QgsPointXY(ptXmax,ymax)

        pt_hash_box = ([ptmin_hash_box, ptmax_hash_box])
        pts_hash_box.append(pt_hash_box)

    multipolylines = QgsGeometry.fromMultiPolylineXY(pts_hash_box)

    #intersection between 45 degrees virtual lines and the polygon features
    for layer_feature in layer.getFeatures() :
        if layer_feature.geometry().intersects(multipolylines):
            pts = layer_feature.geometry().intersection(multipolylines).asMultiPolyline()
            for pt in pts :
                line_in_polygon.append(y)

return QgsGeometry.fromMultiPolylineXY(line_in_polygon)

@qgsfunction(args='auto', group='Custom')
def line_pattern_fill(layer_name, xInterval, feature, parent):
for layer in QgsProject.instance().mapLayers().values():
    if layer.name() == layer_name :
        return _layer_intersection(layer, xInterval)


texts: Setting domains for feature class using ArcGIS Online

texts:
I am looking to setup a basic app that allows field staff to fill out a form and attach a photo in order to track certain field features.

The issue I am having is that ideally I would like to setup domains in order to have drop down menus for filling out the forms. I can certainly set this up through ArcGIS Desktop, however my licensing level does not allow for setting up photo attachments.

If I create a feature class on-line through ArcGIS Online, I can setup photo attachments, however it doesn't appear that I can setup domains via ArcGIS Online. Am I wrong in assuming I can't set domains through ArcGIS Online?

My only options seem to be either an app with no photo attachments but easy to fill in form or an app with photo attachments but a harder to fill in form.

Is it possible to set domains via AGOL?


texts: Calculate file path with field calculator

texts:
Doeos anyone have script out there to calculate the file path name with the field calculator....I am trying to add a field with the file path name so I can keep track of the original source directory of .shp's...

Thanks!


texts:
I am trying to use dynamic IDs so I can create a bounding box over the edges.
I do have a starting road, which I set up before the query. The closest road_id from a city's coordinates. (i.e. 923342 in this case).
From this road I need to calculate the distance to other cities.
The next query works, but it is way too slow (about 15s), since it selects all the 1.2M roads.
SELECT city.id, city.road_id, city.name,
(
    SELECT SUM(roads.km)
    FROM pgr_dijkstra(
        'SELECT id, source, target, cost, reverse_cost FROM roads',
        (SELECT source FROM roads WHERE id = 923342),
        (SELECT target FROM roads WHERE id = city.road_id)
    ) AS dijkstra
    INNER JOIN roads ON roads.id = dijkstra.edge
) AS distance
FROM cities AS city
WHERE city.id IN (2136, 1984, 1850)
ORDER BY distance ASC

Now, in order to speed things up I created a bounding box over the edges:
SELECT city.id, city.road_id, city.name,
(
    SELECT SUM(roads.km)
    FROM pgr_dijkstra(
        'SELECT id, source, target, cost, reverse_cost FROM roads
        (
            SELECT ST_Expand(ST_Extent(geom_way), 0.1) as box
            FROM roads
            WHERE source = (SELECT source FROM roads WHERE id = 923342)
            OR target = (SELECT target FROM roads WHERE id = 12345678)
        ) AS bbox
        WHERE roads.geom_way && bbox.box
        ',
        (SELECT source FROM roads WHERE id = 923342),
        (SELECT target FROM roads WHERE id = city.road_id)
    ) AS dijkstra
    INNER JOIN roads ON roads.id = dijkstra.edge
) AS distance
FROM cities AS city
WHERE city.id IN (2136, 1984, 1850)
ORDER BY distance ASC

(12345678 is a test value retrieved manually).
Now the execution time is around ms. That's great!
The final request is to dynamically change that 12345678 value and replace it with city.road_id. I have no idea how to do that. I've tried using WITH clause, so I can select the roads within the bounding box previously, but I was unable use that virtual table inside the SQL text from pgr_dijkstra.
Looking for something along the lines of:
...
    SELECT SUM(roads.km)
    FROM pgr_dijkstra(
        'SELECT id, source, target, cost, reverse_cost FROM roads
        (
            SELECT ST_Expand(ST_Extent(geom_way), 0.1) as box
            FROM roads
            WHERE source = (SELECT source FROM roads WHERE id = 923342)
            OR target = (SELECT target FROM roads WHERE id = city.road_id) -- [using dynamic ID]
        ) AS bbox
        WHERE roads.geom_way && bbox.box
        ',
        (SELECT source FROM roads WHERE id = 923342),
        (SELECT target FROM roads WHERE id = city.road_id)
    ) AS dijkstra
    INNER JOIN roads ON roads.id = dijkstra.edge
...


texts: GeoServer in Windows x64 - war file installation versus Windows installer?

texts: After running a GeoServer war file installation in Tomcat for some time in Windows Server 2008 x64, I've become increasingly unhappy with Tomcat's chronic memory issues and leaks. Now that the GeoServer Windows installer appears to run in a 64-bit environment as a Windows service when used in tandem with a 32-bit JVM, I'm looking into abandoning the war file in favor of the Windows installer. What are the pros and cons of the Windows installer versus the war file in a 64-bit Windows environment?


texts: How can I create a new layer from the current selection using PyQGIS?

texts:
Is it possible to create a duplicate layer, only displaying the selection from the layer you are duplicating in QuantumGIS?
Similar to ArcGIS -> Create Layer from Selected Features.

I looked at "How to duplicate a layer"

iface = qgis.utils.iface
vl = iface.activeLayer()
iface.addVectorLayer( vl.source(), vl.name(), vl.providerType() )

but can't see how to duplicate only the selected features?


texts: Focal Mean of neighbors using ST_DWithin on points in PostGIS

texts:
I'm trying to calculate focal means on a grid of points in PostGIS.

I know that a raster calculation would be more efficient, but I only need to run the algorithm once to get the data, so I don't mind letting it run and waiting.

If you think there's a much faster way to do this, including conversion to corresponding format and back, please let me know. I'm primarily interested in understanding how to do this in PostGIS (because I'm tyring to learn PostGIS right now), but I'm open to hearing other ideas if you think this is an atrocious thing to do.

If I understand this correctly, the query below should, for each point, find the neighboring points that fall within 360 projection units and give me a count and average value drawn from the neighboring points.

select
    a.grid_id,
    a.pt,
    count(b.pt) as count,
    avg(b.co2perday) as co2perday
from avgs_pts as a
join avgs_pts as b
    on ST_DWithin(a.pt, b.pt, 360)
group by a.grid_id, a.pt;

Does this do what I think it does?
I've already created a spatial index using CREATE INDEX pt_index ON avgs_pts USING gist (pt) and VACUUM (analyze) avgs_pts. Is there anything else I can do to make this query more efficient in PostGIS?


texts: QGIS crashes when running script from Processing Script Editor

texts:
Encountering this weird bug. I'm using QGIS 3.4.7. I am trying to write a script but it crashes when I run it from the editor. But it does not crash when run from the Processing Toolbox. Right now the script does not do anything yet. It just takes an input and outputs it.

Sometimes, instead of crashing, I get this Python error

Traceback (most recent call last):
File "C:/OSGEO4~1/apps/qgis-ltr/./python/plugins\processing\gui\AlgorithmDialog.py", line 247, in on_complete
self.feedback.pushInfo(self.tr('Execution completed in {0:0.2f} seconds').format(time.time() - start_time))
NameError: free variable 'self' referenced before assignment in enclosing scope

Edited to include code.
I didn't include it before because the code is just a copy of the template and it just outputs the same input layer.

def processAlgorithm(self, parameters, context, feedback):
    lpi_lyr = self.parameterAsSource(
        parameters,
        self.INPUT,
        context
    )
    if lpi_lyr is None:
        raise QgsProcessingException(self.invalidSourceError(parameters, self.INPUT))

    feedback.pushInfo('CRS is {}'.format(lpi_lyr.sourceCrs().authid()))

    count = lpi_lyr.featureCount()
    total = 100.0 / count if count else 0

    feats = [feat for feat in lpi_lyr.getFeatures()]
    feedback.pushInfo('Total number of features: {}'.format(len(feats)))

    #PROCESSING
    mem_lyr = QgsVectorLayer("Point?crs="+lpi_lyr.sourceCrs().authid(), 'point', 'memory')
    prov = mem_lyr.dataProvider()
    prov.addAttributes([field for field in lpi_lyr.fields()])
    mem_lyr.updateFields()
    prov.addFeatures(feats)

    ###OUTPUT
    output_lyr = lpi_lyr
    (sink, dest_id) = self.parameterAsSink(
        parameters,
        self.OUTPUT,
        context,
        output_lyr.fields(),
        output_lyr.wkbType(),
        output_lyr.sourceCrs()
    )

    if sink is None:
        raise QgsProcessingException(self.invalidSinkError(parameters, self.OUTPUT))

    features = output_lyr.getFeatures()

    for current, feature in enumerate(features):
        # Stop the algorithm if cancel button has been clicked
        if feedback.isCanceled():
            break
        # Add a feature in the sink
        sink.addFeature(feature, QgsFeatureSink.FastInsert)
        # Update the progress bar
        feedback.setProgress(int(current * total))

    feedback.pushInfo('Finished writing sink')

    return {self.OUTPUT: dest_id}


texts: Create a table with open() in a loop using the SearchCursor and InsertCursor

texts:
I like to create a table with open() in a loop using the SearchCursor and InsertCursor, but I don't know how.
Before I used a few geoprocessing tools to get a field for percent. The Attributtable of the raster looks like this

So I need to write the field percent with the value 8 in a new table to get a result which looks like this:

This is where I'm working on with ArcMap10.2.1:

import arcpy
import os
from arcpy import env

env.overwriteOutput = True
env.workspace = r"D:\Julia_T\projekt\testdataA"

#folder to list
ascFileList = arcpy.ListFiles("*.asc")

print ascFileList
#works
for ascFile in ascFileList:
# geoprocessing steps
ascFileName = os.path.splitext (ascFile)[0]
print ascFileName
#works
rastFile = ascFileName + "_output.img"
print rastFile
#works
#run the tool
arcpy.ASCIIToRaster_conversion (env.workspace + "\\" + ascFile, env.workspace + "\\" + rastFile, "INTEGER")

# Process: Projektion definieren
arcpy.DefineProjection_management (env.workspace + "\\" + rastFile, "GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]")
# Process: Feld hinzufügen
arcpy.AddField_management (env.workspace + "\\" + rastFile,"Prozent","DOUBLE")
# Process: Feld berechnen
arcpy.CalculateField_management(env.workspace + "\\" + rastFile,"Prozent","( [Count] /53959) *100", "VB")
#tabelle bauen
with arcpy.da.SearchCursor(env.workspace + "\\" + rastFile,["Value","Prozent"],""""Value" = 8""") as sCur:
    outputfile = r"D:\Julia_T\output\output.txt"
    with arcpy.da.InsertCursor(outputfile,["Value","Prozent"]) as iCur:
               for row in sCur:
                print row [0], row [1]
                f = open (outputfile, 'w')
                f.write ("Value\nProzent")
                f.close()

How can I get the percent from the value 8 in a new table?


texts: Create buffer-lines in QGIS alongside Vector lines

texts:
I am currently doing a map of a farm, to see which roads needs to be repaired.

Those roads are either one sided fenced, on both sides fenced or not fenced at all. I want to make buffer lines around the already mapped roads (vector lines), which have fences. Something like this:

I hope this is possible to understand. Of course I can draw additional lines alonside those paths... but this is a hellofa work.

I work with QGIS :)

Thanks for your help in advance!

Lukas


texts:
How can I add Dark Matter from carto.com in my QGIS project?

As WMS by example or with a plugin.


texts:
Does anybody know if there is a way to statistical calculate the similarity of two polygons using PostGIS?

I want to know how similar is the polygon (red outline) with the polygon (black outline).


texts: Converting a SpatialPointsDataFrame into km units from m, in R

texts:
I read in a UTM shapefile to R using readOGR. I want to use the resulting SpatialPointsDataFrame to build a mesh in the INLA package, but I've been told that to help with numerical stability, I should convert the geometries to km from m.
So far, I've been using INLA:
grid2 <- readOGR(here::here("PNG_grid.shp"))
mycrs <- fm_sp_get_crs(grid2)
mycrs2 <- fm_crs_set_lengthunit(mycrs, "km")
grid3 <- inla.spTransform(grid2, crs0 = mycrs, crs2 = mycrs2)

...but get an error "Invalid target CRS for SpatialPointsDataFrame". Looking at mycrs2 it's pretty weird. The value in mycrs is "+proj=utm +zone=54 +south +ellps=WGS72 +units=m +no_defs", so I thought this might work:
grid3 <- inla.spTransform(grid2, crs0 = "+proj=utm +zone=54 +south +ellps=WGS72 +units=m +no_defs",
         crs1 = "+proj=utm +zone=54 +south +ellps=WGS72 +units=km +no_defs")

...but I get the same error. Is there anything in sp or rgdal that would help me?


texts:
I have a road network where each link is part of a reference number.
I am trying to sum up the length of roads per reference number and status.
My nested query fails and just sum the total length of roads for each status category and are not split into reference number group.

SELECT number,
(sum(ST_Length(geom)) FROM schema.roads WHERE roads.status IN ('G','V')) as not_started
(sum(ST_Length(geom)) FROM schema.roads WHERE roads.status IN ('T','V')) as started
(sum(ST_Length(geom)) FROM schema.roads WHERE roads.status IN ('D')) as finish
GROUP BY number

I believe I have to build a left join query for this, but using functions and group by within a left join query is a bit too complex with my current experience.
Could anyone please have give me some suggestions?


texts:
I would like to create a distance matrix and have as output the distances expressed as kilometers (km).

Here my current output, which I believe is in degrees and NOT in km:

My coordinate system (CRS) is EPSG:27700, OSGB 1936 / British National Grid, I am working with an Apple computer on QGIS 2.14


texts: Using writeOGR on a list of SpatialPolygons

texts:
This question has been asked before, but I find the answer confusing since there are no reproducible examples that link the question and the answer, and the answer doesn't actually use the writeOGR function.
Use lapply and writeOGR to write multiple shapefiles within a list

I have a list of spatial objects. I want to export all of the using writeOGR. How can I apply this function to the list of spatial objects?

The solution I've come up with applies writeOGR to the outer list element, and I suppose I need to apply it to the inner list element, but I don't know how to do this.

library(GISTools)
library(rgdal)
data(tornados)

torn_list <- list(torn, torn2)
lapply(torn_list, function(x) writeOGR(obj = x, dsn = "tempdir", layer = x, driver = "ESRI Shapefile"))

Error in as.character.default(layer) :
  no method for coercing this S4 class to a vector

Any suggestions?

To clarify, I do want to do this to a list of shape files that have different coordinate reference systems. I would like to create a different shapefile for every element of the list. I realize now I would need to have some way of making the "layer =" argument  for each element of this list as well.


texts: How to Change CartoDB Torque Cat Map Slider Range

texts:
Map slider range includes 83 month/years (2017, 2018, and 2019 are not there) but actual data only has 57 month/years see: https://sildeag.cartodb.com/viz/1e7ee1c6-0cbc-11e6-9822-0e787de82d45/map

Is there any way for a non-programmer to alter this behaviour so that slider range is for actual data values only using the "free" version?


texts:
I installed the latest version of QGIS on Mac OS 10.7.3 with 2 additional frameworks (GDAL Complete and GSL) but there is no main menu bar (not at the top of screen, not at the window of QGIS).

Where is the menu?


texts: Waiting for Python script tool to complete in ModelBuilder?

texts:
I have an Arcgis (10) model which has 2 process that need to be run sequentially.

The 1st process is a python script that creates multiple rasters, but that output is not known to ModelBuilder, therefore I can't use it as a precondition. Or at least not as outline in Example #1 in that page. Maybe there's another route I haven't discovered.

The 2nd process needs those rasters in order to operate. How do I let #2 know when #1 is complete so it may start execution?

To reproduce build a model with:

Process_1: a script which runs PolygonToRaster several times, converting each field in the attribute table to a raster (here, but field list will need to be edited).

Process_2: Run LessThanFrequency Spatial Analyst tool on the results (but really the tool doesn't matter).


texts: Pie Chart Slice Sizing in ArcGIS-Desktop

texts:
I have two fields I would like to show in pie charts on a map.

I'm having trouble understanding how field 2 is represented as a part of field 1. If a polygon in field 2 is over half of field 1, why is it being shown as almost 1/4 of the pie chart? - confused!

For example - this pie chart is showing the light blue field as 80,000 and the dark blue field as 70,000

I read in the help files that it converts your data to map units, so it is showing field 2 based on its proportion of the size it has set as the scale. I think.

But anyway, it is misrepresenting my data...I need to show field 2 as half the pie chart if the figure is half of field 1. How to get around this? I've tried setting is as a ratio but it still doesn't look right.

nm I think I get it after reading excel help. What I need to  do is convert fields 1 and 2 to a proportion of 100, so 80,000 would become 5 and 60,000 would become 4. Easy way to do this using field calculator?


texts: Extracting Grids containing 50% polygons (from land use layer) using QGIS?

texts:
I have a vector layer of Parcels for a certain area. I need to make grids of 100m * 100m over the parcel layer and then extract those grids which contain 50% or parcels by area.

I am currently using QGIS. The have tried taking an intersection of the Grids and the parcels and then run a selection query to select the ones where the area of the intersection (of parcel and grid) should be greater than 50m.

Is there any direct tool or plugin that serves the same purpose?


texts: Select lines in within a polygon (buffer) with PostGIS

texts:
I have one polygon layer "polygons" and one line layer "lines". I would like to create a layer with all line parts that are inside the polygon layer.
First I cut the lines in the intersection of the polygon layer. I can see that all lines are cut in the intersection of the polygons.
then I would like to select the line parts inside the polygon that has a length longer than approximate the double width of the polygon (buffer).
First image shows the cutted lines, I've highlighted one line for you to see. the black area is the polygons (buffer).

I use the following sql to select the "inside"-lines:
CREATE TABLE public."insidelines" AS SELECT l.*
    FROM "lines" l
        LEFT JOIN "polygons" p ON ST_Contains(p.the_geom, l.the_geom)
    WHERE p.the_geom IS NOT NULL AND ST_Length(l.the_geom) > 2000;

 ALTER TABLE public."insidelines" ALTER COLUMN the_geom TYPE geometry(LINESTRING, 32633) USING ST_SetSRID(the_geom, 32633);

the result is:

It is only the yellow line that was found! Any idea why?
Both layers are in UTM 33 (32633).
I'm using PostGIS ver 2.5


texts:
I just got my very first NMF file and I am stumped on how to get it into my database. I have been able to handle anything thrown at me so far using either QGIS and ogr2ogr to get my data into PostGIS.

It is a polygon layer, so just opening it in Excel is not enough and I do not have ArcGIS to run some extensions I heard may get the job done.

After Export to Excel, I can write some code to merge the related rows, make some WKT and construct some insert statements, but I would like to the extra work if possible.


texts: What is the easiest way to add up hourly NetCDF files to obtain monthy files?

texts:
I have about 1 terrabyte worth of hourly .nc files from 1951 to 2018. I've already reorganized the hourly files into a monthly directory. I now need to add up all the  hourly .nc files in each monthly directory to make monthly .nc file.
What is the quickest way to do this?

If there's no way to do it without coding, is it possible on python?
right now my files look like this:

C:\Users\Ismail\Documents\Data\1951\1\1951-1-1-0000.nc
C:\Users\Ismail\Documents\Data\1951\1\1951-1-1-0100.nc
...
...
...
C:\Users\Ismail\Documents\Data\2018\12\2018-12-31-2200.nc
C:\Users\Ismail\Documents\Data\2018\12\2018-12-31-2300.nc

I want it to look like this:

C:\Users\Ismail\Documents\Data\1951\1951-1.nc
C:\Users\Ismail\Documents\Data\1951\1951-2.nc
...
...
...
C:\Users\Ismail\Documents\Data\2018\2018-11.nc
C:\Users\Ismail\Documents\Data\2018\2018-12.nc

Thank you for your help!


texts:
I'm working with ArcGIS Desktop 10.7 ModelBuilder and was wondering if there is a way to automatically export the TIFFs my model creates into a rendered TIFF.

Below is a photo of what I'm trying to recreate.

Is there a code in arcpy that I can use to recreate this export function?


texts: QGIS crashing when trying to change style attribute

texts:
QGIS started crashing w Ver 2.18.3-1

If I open a project from a previous version it crashes w error QGIS bin has stopped working.

If opening only a new CSV file or postgresql DB layer it will load the new layer and crash when I attempt to open the properties of the layer.

I did re install ver 2.14  stable Corp release and can get project to open, but those projects will crash QGIS when opening the style attribute.

Has anyone one else experienced this or have a fix?

"qgis-ltr-bin.exe has stopped working"

Windows 10 Enterprise Error log:

Faulting application name: qgis-ltr-bin.exe, version: 0.0.0.0, time stamp: 0x58799528
Faulting module name: ntdll.dll, version: 10.0.14393.447, time stamp: 0x5819bc32
Exception code: 0xc0000005
Fault offset: 0x000000000004839c
Faulting process id: 0x23dc
Faulting application start time: 0x01d28853a2c42cdc
Faulting application path: C:\PROGRA~1\QGIS2~1.14\bin\qgis-ltr-bin.exe
Faulting module path: C:\WINDOWS\SYSTEM32\ntdll.dll
Report Id: acab3985-0a94-4541-8c50-66053a0995b8
Faulting package full name:
Faulting package-relative application ID:


texts: OpenLayers 4.6.5 GeoJSON WGS84 not mapped with a different projection

texts:
I have a map projected on EPSG:3857. I must add ocean areas around the globe. Areas come from different GeoJSON files that are added via a select. Most of geojson files have a projection WGS84 (EPSG:4326). Some of them have a orthographic projection (those in the Arctic and Antarctic). Via another select I should be able to change view projection (following the example).

I want that layers change their projection according to the view projection chosen. But I failed to do it, as when I change view projection with the code below, the layers remain visible with their own native projection.

My code is

//projection for north pole
proj4.defs("ORTHO","+proj=ortho +lat_0=85.5333333333 +lon_0=-72.53333333339999 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs ");

//area style
var defaultStyle =
 new ol.style.Style({
  fill: new ol.style.Fill({
  color: [251,183,82,0.8]
 }),
 stroke: new ol.style.Stroke({
 color: [251,118,74,1],
 width: 1
 })
});

//layers
/*NON POLAR*/
var source_ALL_NON_POLAR_GEOJSON_IMMA_WGS84 =
  new ol.source.Vector({
  url: 'polar_examples/new/ALL_NON_POLAR_GEOJSON_IMMA_WGS84.geojson',
  format: new ol.format.GeoJSON({
  }),
dataProjection: 'EPSG:4326'
});

/*POLAR*/
var source_ortho_North =
  new ol.source.Vector({
  format: new ol.format.GeoJSON({
  }),
url: 'polar_examples/new/ARCTIC_GEOJSON_FAKE_ORTHOGRAPHICnorth.geojson',
});

var layer_ortho_North =
  new ol.layer.Vector({
  source: source_ortho_North,
  style: defaultStyle,
  name: 'ortho_North',
  visible: true,
});

var ocean_map =
  new ol.layer.Tile({
    source: new ol.source.XYZ({
    url: 'https://services.arcgisonline.com/arcgis/rest/services/' +
    'Ocean/World_Ocean_Base/MapServer/tile/{z}/{y}/{x}'
    }),
  visible: true,
  });
var map = new ol.Map({
    layers: [],
    target: 'map',
    view: new ol.View({
      //projection:'EPSG:3857', //view projection mercator
      projection: "ORTHO", //view projection ortho
      center: ol.proj.fromLonLat([0, 90], 'EPSG:3408'),//ortho center
      //center: ol.proj.transform([-97.6114, 38.8403], 'EPSG:3857', 'EPSG:4326'),//mercator center
      zoom: 2,
      minZoom: 1
    })

//adding layers
  map.addLayer(ocean_map);
  map.addLayer(layer_ortho_North);
  // map.addLayer(layer_ALL_NON_POLAR_GEOJSON_IMMA_WGS84);

To test it I used to comment each projection and layer to be added, in order to see what happens (the select switcher are not implemented yet).

Reading about the problem I found that, when getting the source, it is possible to set a new projection type to the GeoJSON file features matching the one used for the map.

new ol.source.Vector({
 features: new ol.format.GeoJSON().readFeatures(source, {
 dataProjection: 'EPSG:4326', //geojson projection
 featureProjection:'EPSG:3408', //view projection
 })
});

However I don't have a clue about what to use as "source" for the readFeatures function. I tried to add the GeoJSON file name but it doesn't work.


texts: Not correct numeric value error using Multilevel B-Spline interpolation with SAGA and QGIS?

texts:
This is the generated log:

Processing algorithm…
Algorithm 'Multilevel b-spline interpolation (from raster)' starting…
Input parameters:
{ 'DATATYPE' : 0, 'EPSILON' : 0.0001, 'GRID' : '/Users/PassifloraTsui/Desktop/study_area/ssTemp.tif', 'LEVEL_MAX' : 11, 'METHOD' : 1, 'TARGET_OUT_GRID' : '/Users/PassifloraTsui/Desktop/study_area/ssTemp_BSpline.sdat', 'TARGET_USER_FITS' : 0, 'TARGET_USER_SIZE' : 30, 'TARGET_USER_XMIN TARGET_USER_XMAX TARGET_USER_YMIN TARGET_USER_YMAX' : '132.20166999999998,145.51116999999996,-21.831979999999955,-9.220450000000007 [EPSG:4326]', 'UPDATE' : False }

grid_spline "Multilevel B-Spline Interpolation (from Grid)" -TARGET_DEFINITION 0 -GRID "/var/folders/nb/kxpdv6xn68xbmxc5l_8rgbcw0000gn/T/processing_abf242d6ade948aaa3a64539df235165/dbc876f6b3d641d7ac65bc37ae745e79/ssTemp.sgrd" -METHOD 1 -EPSILON 0.0001 -LEVEL_MAX 11.0 -UPDATE false -DATATYPE 0 -TARGET_USER_XMIN 132.20166999999998 -TARGET_USER_XMAX 145.51116999999996 -TARGET_USER_YMIN -21.831979999999955 -TARGET_USER_YMAX -9.220450000000007 -TARGET_USER_SIZE 30.0 -TARGET_USER_FITS 0 -TARGET_OUT_GRID "/Users/PassifloraTsui/Desktop/study_area/ssTemp_BSpline.sdat"
'11.0' is not a correct numeric value for option 'LEVEL_MAX'.
____________________________

##### ## ##### ##
### ### ## ###
### # ## ## #### # ##
### ##### ## # #####
##### # ## ##### # ##
____________________________

SAGA Version: 2.3.1

____________________________
library path: /Applications/QGIS3.app/Contents/MacOS/lib/saga/
library name: libgrid_spline
library : Spline Interpolation
tool : Multilevel B-Spline Interpolation (from Grid)
author : O.Conrad (c) 2006
____________________________

Usage: saga_cmd grid_spline 5 [-GRID <str>] [-TARGET_DEFINITION <str>] [-TARGET_USER_SIZE <double>] [-TARGET_USER_XMIN <double>] [-TARGET_USER_XMAX <double>] [-TARGET_USER_YMIN <double>] [-TARGET_USER_YMAX <double>] [-TARGET_USER_FITS <str>] [-TARGET_TEMPLATE <str>] [-TARGET_OUT_GRID <str>] [-METHOD <str>] [-EPSILON <double>] [-LEVEL_MAX <num>] [-UPDATE <str>] [-DATATYPE <str>]
-GRID:<str> Grid
Grid (input)
-TARGET_DEFINITION:<str> Target Grid System
Choice
Available Choices:
[0] user defined
[1] grid or grid system
Default: 0
-TARGET_USER_SIZE:<double> Cellsize
Floating point
Minimum: 0.000000
Default: 1.000000
-TARGET_USER_XMIN:<double> Left
Floating point
Default: 0.000000
-TARGET_USER_XMAX:<double> Right
Floating point
Default: 100.000000
-TARGET_USER_YMIN:<double> Bottom
Floating point
Default: 0.000000
-TARGET_USER_YMAX:<double> Top
Floating point
Default: 100.000000
-TARGET_USER_FITS:<str> Fit
Choice
Available Choices:
[0] nodes
[1] cells
Default: 0
-TARGET_TEMPLATE:<str> Target System
Grid (optional input)
-TARGET_OUT_GRID:<str> Target Grid
Grid (output)
-METHOD:<str> Method
Choice
Available Choices:
[0] without B-spline refinement
[1] with B-spline refinement
Default: 1
-EPSILON:<double> Threshold Error
Floating point
Minimum: 0.000000
Default: 0.000100
-LEVEL_MAX:<num> Maximum Level
Integer
Minimum: 1
Maximum: 14
Default: 11
-UPDATE:<str> Update View
Boolean
Default: 0
-DATATYPE:<str> Data Type
Choice
Available Choices:
[0] same as input grid
[1] floating point
Default: 1
Execution completed in 2.38 seconds
Results:
{'TARGET_OUT_GRID': '/Users/PassifloraTsui/Desktop/study_area/ssTemp_BSpline.sdat'}

Loading resulting layers
The following layers were not correctly generated.

/Users/PassifloraTsui/Desktop/study_area/ssTemp_BSpline.sdat

You can check the 'Log Messages Panel' in QGIS main window to find more information about the execution of the algorithm.

and I've checked the documentation of this algorithm in SAGA website, which said the LEVEL_MAX should be an integer, not a float number, however, QGIS automatically converts it to a float number, which leads to the failure.


texts: Determine the Degree of Data Lossiness when Converting Data

texts:
I'm working through some of the ESRI tutorials on Data Interoperability and it briefly mentions the issue of data conversion being a lossy process.

Some questions that come to mind are:

What is being lost?
Does the degree of lossiness depend on data type and format?
Which conversion combinations are most affected?
Is there situations where conversion should be avoided at all costs?
What techniques are used to determine the degree of lossiness?

Feel free to add to the list or express any other concerns.


texts:
I want to make a custom application which includes an option just like the "Select Features by Polygon" feature in the QGIS. I am able to draw the polygon using following code:

self.rb = QgsRubberBand(self.canvas, QGis.Polygon)
vertex = QgsPoint(point.x(), point.y())
self.rb.addPoint(QgsPoint)

I want to add the feature in which the polygon sort of a sticks to the mouse pointer until the polygon is completed just like the "Select Features by Polygon" in QGIS which is more user friendly. Any ideas on how to do this?


texts: Clipping/ Zooming in on shapefile using R

texts:
I'm working on a shapefile of all the villages of Portugal from this link.
I've uploaded it to R:
#read in shapefile of portugal
portugal = readOGR(dsn = "./communities", layer = "Cont_AAD_CAOP2017", verbose=FALSE)

And I've got the coordinate reference system information
# get coordinate reference system information
st_crs(portugal)

Coordinate Reference System:
  User input: ETRS89 / Portugal TM06
  wkt:
PROJCRS["ETRS89 / Portugal TM06",
    BASEGEOGCRS["ETRS89",
        DATUM["European Terrestrial Reference System 1989",
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4258]],
    CONVERSION["Portugual TM06",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",39.6682583333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",-8.13310833333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["easting (X)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["northing (Y)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Topographic mapping (medium scale)."],
        AREA["Portugal - mainland - onshore."],
        BBOX[36.95,-9.56,42.16,-6.19]],
    ID["EPSG",3763]]

I now want to create a second shapefile using only a region of Portugal, as delimited by some coordinates, namely
W 8°46'00"-- W 7°40'00" / N 41°08'00" -- N 40°39'00"

How can I do this?
Once again, I'm using R


texts: Is it possible to add a map which is published to QGIS cloud, on another web page? As a script code, like customised Google maps for example


texts: Which commercial off-the-shelf GIS server products work on Linux?

texts:
I am looking for suitable COTS GIS server product for Linux platform.  The product should allow 500 Clients concurrently and should have Java / C++ development APIs.
So far I found ESRI ArcGIS server which supports for Linux.  I need few more products which have good performance like ESRI.  Kindly share the list of  available COTS (Commercial off the shelf) GIS server products.


texts:
Following the pgRouting workshop here, I have the routing query below:

$sql = "SELECT gid, start_id AS source, end_id AS target, the_geom,
                  distance(the_geom, GeometryFromText(
                       'POINT(".$lonlat[0]." ".$lonlat[1].")', 4326)) AS dist
                 FROM ".TABLE."
                 WHERE the_geom && setsrid(
                       'BOX3D(".($lonlat[0]-0.1)."
                              ".($lonlat[1]-0.1).",
                              ".($lonlat[0]+0.1)."
                              ".($lonlat[1]+0.1).")'::box3d, 4326)
                 ORDER BY dist LIMIT 1";

This query fails with the error: Query failed: ERROR:  Operation on two geometries with different SRIDs. Now, I've been looking around and came across this. I've tried the suggestions on this post, but the error persists. Any idea on how to fix this will be appreciated.

UPDATE:
Running SELECT st_srid(the_geom) FROM table_name against my database table returns -1. In my query and map, the assumed external projection is EPSG:4326. Is changing this value likely to help? Is there any harm in changing the srid value directly in the database? How can I go about changing the srid value in the database if it doesn't present any risk?

UPDATE:
I had to change the srid of the column from the UNKNOWN indicator of -1 to 4326. The answer at this StackOverflow post was very helpful in this. Also, the question in this SE GIS post was instrumental in helping me in the right direction.


texts:
I'm trying to register a listener event on "loadend" for my google maps layer on OpenLayers. But for some reason it's not firing. There must be something small i'm overlooking. What is wrong with my code?

Code:

<!DOCTYPE html>

    <meta http-equiv="Content-Type" content="text/html; charset=utf-8">
    <meta name="viewport" content="width=device-width; initial-scale=1.0; maximum-scale=1.0; user-scalable=0;">
    <meta name="apple-mobile-web-app-capable" content="yes">
        <title>OpenLayers Google (v3) Layer Example</title>
    <link rel="stylesheet" href="../js/OpenLayers/theme/default/style.css" type="text/css" />
    <link rel="stylesheet" href="style.css" type="text/css" />


texts: Finding min and max raster values using polygons in ArcMap?

texts:
I have a raster dataset that encompasses an area which has previously been classified into vegetation types. I have the classification polygons and want to lay them on top of a vegetation index image and find the range of vegetation index values within each polygon.

Is this doable in ArcMap 10.4?

The closest I've gotten is to create points for each pixel using Raster to Point, then perform a spatial join, but all that gives me is the mean values for each polygon.


texts: Can QGIS use a non-spatial Oracle Database?

texts:
I'm building my first plugin in qgis, that needs to read some coordinates from an oracle database (not oracle spatial database) and create a shapefile with this points. I'm searching but only find to in oracle spatial!
Is there a way to do this in a simple oracle database?


texts:
I am looking for a source of maps that identifies lane-level streets (i.e. maps that encode all lanes of a street). Mapbox has recently released Mapbox Drive, but I could not figure out how to access these maps. I hope if anyone could assist on finding how to get access to these maps.
Also, are there other resources that I can utilize to get lane-level details? I am aware that OpenStreetMaps encodes the number of lanes of a street, but does not have polylines for these lines.


texts:
All my layers in the qgs-project-file are working fine as WMS-layers, but I only can't provide the virtual layer.
The virtual layer is an ODBC-Connection of an Access-Database which contains one column with polygons as WKT and works like a charm.
My QGIS-server version is 3.4.4-2

Isn't it possible to provide a virtual layers via the QGIS-server as a WMS?

Or is there a solution?


texts: Getting inaccurate results converting from New York State projection to NAD83 with Python's pyproj

texts:
I've been trying to do a standard coordinate translation between New York State Long Island (ESRI:102718) and latitude/longitude (EPSG:4269)...but I keep getting the wrong lat/lng results.

I don't think this has anything to do with the pyproj library in particular, but maybe I'm missing something in the transformation process.

Sample data comes from NYC.gov: http://www1.nyc.gov/site/planning/data-maps/open-data/dwn-selfac.page

Its documentation says:

  The X coordinate of the XY coordinate pair which depicts the approximate location of the facility, program, or parkland property. (Also see Y COORDINATE). Facilities' X and Y coordinates must be used together for mapping selection. The XY coordinates are expressed in the New York Long Island State Plane coordinate system.

Here's a sample record:

```

Sample record:

"3861130100",1,"MN",104,"MANHATTAN COMMUNITY DISTRICT 4",731,1,"YAI - YOUNG ADULT INSTITUTE",
"460 W 34 St","10001",3807,7,72,,98,,,,91,57," 103","0103",3,2,10,5500,1007310001.000,
984533.000,213998.000,"L"
```

So, the python code:

import pyproj
NYSP1983 = pyproj.Proj(init="ESRI:102718")
x, y = (984533, 213998)
NYSP1983(x, y, inverse=True)

That inverse projection should result in the lat/lng coordinates. However, the result is nowhere near New York:

       (-65.75200831429146, 41.80318458470038)

According to Google geocoder, these are the proper coordinates:

    location: {
        lat: 40.7541045,
        lng: -73.9989178
    }

Here's another way to do transformations in pyproj, going from latlng to the specified X/Y projection, but still no dice:

import pyproj
NYSP1983 = pyproj.Proj(init="ESRI:102718")
lng, lat = (-73.9989178, 40.7541045)
NYSP1983(lng, lat)

The result:

            (300091.3883012832, 65232.600609979076)

Doing an explicit transform between systems doesn't help either:

import pyproj
NAD83 = pyproj.Proj(init="EPSG:4269")
NYSP1983 = pyproj.Proj(init="ESRI:102718")
x, y = (984533, 213998)
pyproj.transform(NYSP1983, NAD83, x, y)

         (-65.75200831429146, 41.80318458470038)

What am I missing here?

Edit

Seems that I'm missing some parameters that need to be set up, such as False_Easting...but I don't know how in pyproj (or in any system) how to specify that. Namely, how do I use what's in this PROJCS string (example here: https://gist.github.com/veltman/6be84d216376814cdc67) to initialize the desired projection?


texts:
Is there a way with Leaflet in R to have layer controls unchecked by default? I've got an addCircleMarkers layer that I want unchecked when the map loads. Here's the code snippet (I know it's not reproducible, but I'm hoping there's a simple way to do this that doesn't require running all the code).

  mymap%>%addCircleMarkers(data=mydata, weight = 1,color="red",radius = ~sqrt(DataPoint),group="NewLayer")%>%
  addLayersControl(overlayGroups = c("NewLayer"),options = layersControlOptions(collapsed = FALSE))

I found this posting that describes it using leaflet directly in js, but I'm not an expert coder in js and am relying on the R interface:
Check/uncheck overlays in groupedlayercontrol.js


texts: Voronoi polygons tool not working in QGIS

texts:
So I've been trying to build up some Voronoi diagram out of my point dataset but I can't seem to get it to work.

When I try it on QGIS 3.4 I get the following error:

  Traceback (most recent call last): File "/Applications/QGIS3.app/Contents/MacOS/../Resources/python/plugins/processing/algs/qgis/VoronoiPolygons.py",
  line 166, in processAlgorithm self.tr('There were no polygons
  created.'))
  _core.QgsProcessingException: There were no polygons created

If I try it on QGIS Mac if it makes a difference. Could use PostgreSQL if needed.


texts:
The tmap_icons() function within the tmap package allows users to call a png file into tmap to be used on a map. However, I can't seem to find out how to add custom symbols to the tmap legend. See the reprex below showing the issue:
library(tmap)
#> Warning: package 'tmap' was built under R version 3.6.3
library(sf)
#> Warning: package 'sf' was built under R version 3.6.3
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
library(tidyverse)
#> Warning: package 'tidyverse' was built under R version 3.6.2
#> Warning: package 'tidyr' was built under R version 3.6.2

tdir=tempdir()
Stateshpurl= "https://www2.census.gov/geo/tiger/GENZ2018/shp/cb_2018_us_state_5m.zip"
if(file.exists(paste(tdir,"/cb_2018_us_county_5m/cb_2018_us_county_5m.shp",sep=""))==FALSE){
 download.file(Stateshpurl, destfile = file.path(tdir,"State_Boundaries.zip"))
 unzip(file.path(tdir,"State_Boundaries.zip"),exdir = tdir)}
US_states=read_sf(paste(tdir,"/cb_2018_us_state_5m.shp",sep=""))

US_states = st_transform(US_states, 5070)

L48 = US_states %>%
 filter(STATEFP != "02") %>%
 filter(STATEFP != "15") %>%
 filter(STATEFP != "60") %>%
 filter(STATEFP != "64") %>%
 filter(STATEFP != "66") %>%
 filter(STATEFP != "68") %>%
 filter(STATEFP != "69") %>%
 filter(STATEFP != "70") %>%
 filter(STATEFP != "72") %>%
 filter(STATEFP != "74") %>%
 filter(STATEFP != "78")

icon_url="http://29.media.tumblr.com/tumblr_m0q2g8mhGK1qk6uvyo1_500.png"
download.file(icon_url, destfile = file.path(tdir,"Image.png"),mode = 'wb')

Icon=tmap_icons(file=paste(tdir,"/Image.png",sep=""))

tm_shape(L48)+
 tm_borders()+
tm_symbols(shape=Icon)

tm_shape(L48)+
  tm_borders()+
  tm_symbols(shape=Icon)+
  tm_add_legend(type="symbol",shape = Icon)
#> Error in shapes[i] - 999: non-numeric argument to binary operator

Perhaps I am going about this wrong, but it seems like there should be a way to add one of the custom symbols created with tm_icons() into the legend.


texts: Measure distance in km of polyline in QGIS

texts:
As much as I have looked numerous times for how to measure vector polyline in QGIS I can not find the way. This seems like an easy thing to do but the answer eludes me. Dah!
After I have traced out a route as a vector polyline on a map layer I sure would like to know how long in kilometres it is. The actual total line length with many points winding along a path.
And not by taking a Ruler tool to measure crudely the distance as you see on Google Maps.


texts: How to change OpenLayers Map styles with bright.json

texts:
I got problem with using Bright.json I caant load that styles to map.
So far I got

let vLayer = new ol.layer.Vector({
      source: new ol.source.Vector({
      url: "/www/bright.json",
      format: new ol.format.GeoJSON()
    })
  })

  let map = new ol.Map({
    target: "map",
    layers: [
      new ol.layer.Tile({
        source: new ol.source.OSM()
      }), vLayer
    ],
    view: new ol.View({
      center: ol.proj.fromLonLat([14.437800, 50.075539]),
      zoom: 12
    })
  });


texts:
I know that it is possible to make a line or polygon with dashed stroke. It can be done with the following way, for example:

    var myStyle = new ol.style.Style ({
      stroke: new ol.style.Stroke({
      color: 'rgba(221, 99, 0, 0.5)',
      width: 3,
      lineDash: [40,10] // <- [dashLength, dashSpace]
    });

I would like to know it is possible to put a pattern or a text (as a pattern) in the stroke drawn, in OpenLayers 3. I painted a picture which represents what I would like:

I tried to look for some options here but there is no any.

If it is can't be done with the standard OL3 libary, are there any external library for this matter?

I tried to look for some solution over this site but I could not find any. I am of the opinion that this "task" can only be done by making a "hook" in the canvas renderer of OpenLayers 3.


texts:
Just pulled QGIS from git and am trying to compile it, here's the relevant error:

[ 89%] Building CXX object python/CMakeFiles/python_module_qgis_core.dir/core/sipcorepart3.cpp.o
In file included from /home/noah/qgis/qgis-git/qgis/Quantum-GIS/python/core/qgsmaplayerregistry.sip:10:0:
/home/noah/qgis/qgis-git/qgis/Quantum-GIS/python/../src/core/qgsmaplayerregistry.h: In member function ‘virtual void sipQgsMapLayerRegistry::connectNotify(const char*)’:
/home/noah/qgis/qgis-git/qgis/Quantum-GIS/python/../src/core/qgsmaplayerregistry.h:252:10: error: ‘virtual void QgsMapLayerRegistry::connectNotify(const char*)’ is private
     void connectNotify( const char * signal );
          ^
/home/noah/qgis/qgis-git/qgis/Quantum-GIS/build-master/python/core/sipcorepart3.cpp:21714:46: error: within this context
         return;
                                              ^
make[2]: * [python/CMakeFiles/python_module_qgis_core.dir/core/sipcorepart3.cpp.o] Error 1
make[1]: * [python/CMakeFiles/python_module_qgis_core.dir/all] Error 2
make: * [all] Error 2

I'm currently downgrading python as suggested on the AUR page-- I don't understand the problem, I hope this is the right thing to do.


texts: Using inline variable substitution output and input name problem to join two Python script tools together in ModelBuilder?

texts:
I have a script I have exported from ModelBuilder and am trying to incorporate at the beginning a script kindly written by @BERA.  The process is to first convert a polygon shapefile to a raster and do a cost distance analysis, set null values in the resultant accumulated cost distance and then summarise the accumulated cost distance inside a polygon of interest. At the beginning, I need to iterate through the fields in an attribute table to provide the values for the cells in the cost surface raster that goes into the cost distance analysis, and the whole process to repeat for each raster. I am having difficulty with providing inline variables to give names for the input and output names for the cost analysis and the set null values and summarize zones. @BERA provided the means for the raster outputs to be named after the field names – which is perfect, but being completely new to python I can’t get the two scripts to flow together because of the input and output names.  I have provided the two scripts without my attempts!(sorry about the long paths).

This is written by @BERA:`

Import arcpy,os
polygons=r'C:\folder\polygons.shp'
sources=r'C:\folder\sources.shp'

fieldlist=['fodyexp1','fodyexp2','fodyexp3','fodyexp4'] #extend/edit to include all your fields

for item in fieldlist: #For each field specified in fieldlist perform all below:
FeatureToRaster_conversion(in_features=polygons, field=item, out_raster=r'in_memory\tempraster')
#Dont know what you want to do next, cost distance?:
outCostDist = CostDistance(sources, r'in_memory\tempraster')
outCostDist.save(os.path.join(r'C:\outputfolder','Raster_'+item)) #The outputs will be named after field, for example Raster_fodyexp1
#and so on...`

This is my cost distance script:

    # Import arcpy module
import arcpy

# Check out any necessary licenses
arcpy.CheckOutExtension("spatial")

# Set Geoprocessing environments
arcpy.env.scratchWorkspace = "C:\\rod.gdb"
arcpy.env.snapRaster = "scores_fody_m"
arcpy.env.extent = "534630.586094391 7813776.96082858 552599.503136767 7825252.24043838"
arcpy.env.workspace = "C:\\rod.gdb"

# Local variables:
habitat = "habitat"
fodyrange_current = "fodyrange_current"
AQ = "AQ"
hab_sn_1 =
"D:\\gis_data\\rodrigues\\Rod_analysis\\iterate_model\\hab_ras\\hab_sn1"
dist_sn_1 =
"D:\\gis_data\\rodrigues\\Rod_analysis\\iterate_model\\dis_raster\\dist_sn_1"
backlink_fmaq = ""

  nullcost_sn_1 =

"D:\\gis_data\\rodrigues\\Rod_analysis\\iterate_model\\setnull_cost_r\\nullcost_sn_1"

zone_sn_1 = "D:\\gis_data\\rodrigues\\Rod_analysis\\iterate_model\\zonal_tables\\zone_sn_1"

# Process: Polygon to Raster
arcpy.PolygonToRaster_conversion(habitat, "warb_mean", hab_sn_1, "CELL_CENTER", "NONE", "20")

# Process: Cost Distance
arcpy.gp.CostDistance_sa(fodyrange_current, hab_sn_1, dist_sn_1, "",
backlink_fmaq)

# Process: Set Null
arcpy.gp.SetNull_sa(dist_sn_1, dist_sn_1, nullcost_sn_1, "\"value\" = 0")

# Process: Zonal Statistics as Table
arcpy.gp.ZonalStatisticsAsTable_sa(AQ, "id", nullcost_sn_1, zone_sn_1, "DATA", "ALL")


texts: Leaflet / Google Map baselayer / markers not visible

texts:
New to Leaflet, currently trying to migrate a project from native Google Maps to Leaflet; having been drawn to Leaflet by the layer management and options of using better 'roadmap' style basemaps.

Have downloaded the plugin to use Google Map tiles as a base map for Leaflet,and its working great as a base layer.

However I have found a problem, and I'm not sufficiently experienced with Leaflet to understand why.  I have markers which are visible with any other base layer - but not with Google. Have tried forcing Z order high (1000 etc.) but it makes no difference. Opacity of the markers is full 1.0

Markers are visible on all other layers (e.g. OSM), when I switch base to Google (satellite or map) they're no longer shown, but switch back to a non-Google base and they're there again.

The reason I'm pursuing the Google base layer is that Google is (as far as I'm aware) the only provider of detailed consistent aerial photographic tiles for my region of interest. Wish to use other maps but have option to switch to a Google 'satellite' view via the layer manager.

Currently using Leaflet v0.4.5 ( but have also tried v0.4 ), same symptoms occur using different browsers (IE, FF etc.).

Would be grateful for any assistance.


texts: Extruding buildings from OS MasterMap data

texts:
I've downloaded OS MasterMap with the intention of extruding buildings to their respective heights. I wrongly assumed that I would receive shapefiles of the buildings which I'd then extrude using the height data in the attribute table but instead comes as a shapefile called 'TopographicArea'.

I've imported this shapefile and linked the included building heights csv file but how do I go about extracting just the buildings and extruding them to their heights?


texts: How to calculate the overlap of two polygons in PostGIS

texts:
I am trying to calculate the percentage overlap of two vector polygon layers and saving the result as a separate column. The code below is working:

CREATE TABLE parcels_all_shapefile_new AS
SELECT DISTINCT a.*, COALESCE((ST_Area(ST_INTERSECTION(a.wkb_geometry, b.wkb_geometry))/ST_Area(a.wkb_geometry))*100, 0) AS aw
FROM parcels_all_shapefile_flood2 a LEFT OUTER JOIN filter_ancientwoodlands b ON ST_Overlaps(a.wkb_geometry, b.wkb_geometry) OR ST_Within(b.wkb_geometry, a.wkb_geometry);

However I am getting an extra ~50,000 rows than the original table. The DISTINCT function has worked for me in the past to deal with this issue, when multiple polygons are intersecting with each other.

Can anyone see what I am missing? Also does anyone know how to deal with a situation when a polygon overlaps multiple polygons and you want to add the percentages together as a single result?

UPDATE: The below answer worked except I had to change INNER JOIN to LEFT OUTER JOIN, here is the code:

CREATE TABLE parcels_all_shapefile_aw AS
SELECT
    a.ogc_fid, a.poly_id, a.title_no,
    SUM(COALESCE((ST_Area(ST_INTERSECTION(a.wkb_geometry, b.wkb_geometry))/ST_Area(a.wkb_geometry))*100, 0)) AS aw
FROM parcels_all_shapefile_flood2 a
    LEFT OUTER JOIN filter_ancientwoodlands b ON
        ST_Overlaps(a.wkb_geometry, b.wkb_geometry) OR ST_Within(b.wkb_geometry, a.wkb_geometry)
GROUP BY a.ogc_fid, a.poly_id, a.title_no


texts: QgsTask.fromFunction not running on_finished method unless an exception is raised in the task

texts:
I'm trying to create a task to process some geometry intersections in a QGIS plugin I've created. I'm following 15.2.2 Task from function from the PyQgis Developer cook book exactly as its shown. I can get the tasks to run as expected but the completed function is never called unless an exception is raised from the doSomething function(see code from the cook book below).
When I run the cook book code from the Python Console in Qgis it works properly, showing the Task Name and iterations in the Log Messages after the task has finished. I've tried adding self to the start of each task and including them in the class that I'm calling the tasks from as suggested in this post.
I have used my own task classe before in the same plugin which have worked well but I'm at the point were I've invested a fair bit of time trying to get this way to work that I want to see it through. Is there any way I can get the specified function to run once a task has successfully completed within my plugin?
Edit: I actually looked through the code for my custom task class and I'm using signals to pass information back to the main thread as the finished method in the class was not being called upon task completion. The workaround is ok though it would be good to understand why it is not working as intended.
import random
from time import sleep

MESSAGE_CATEGORY = 'TaskFromFunction'

def doSomething(task, wait_time):
    """
    Raises an exception to abort the task.
    Returns a result if success.
    The result will be passed, together with the exception (None in
    the case of success), to the on_finished method.
    If there is an exception, there will be no result.
    """
    QgsMessageLog.logMessage('Started task {}'.format(task.description()),
                             MESSAGE_CATEGORY, Qgis.Info)
    wait_time = wait_time / 100
    total = 0
    iterations = 0
    for i in range(100):
        sleep(wait_time)
        # use task.setProgress to report progress
        task.setProgress(i)
        arandominteger = random.randint(0, 500)
        total += arandominteger
        iterations += 1
        # check task.isCanceled() to handle cancellation
        if task.isCanceled():
            stopped(task)
            return None
        # raise an exception to abort the task
        if arandominteger == 42:
            raise Exception('bad value!')
    return {'total': total, 'iterations': iterations,
            'task': task.description()}

def stopped(task):
    QgsMessageLog.logMessage(
        'Task "{name}" was canceled'.format(
            name=task.description()),
        MESSAGE_CATEGORY, Qgis.Info)

def completed(exception, result=None):
    """This is called when doSomething is finished.
    Exception is not None if doSomething raises an exception.
    result is the return value of doSomething."""
    if exception is None:
        if result is None:
            QgsMessageLog.logMessage(
                'Completed with no exception and no result '\
                '(probably manually canceled by the user)',
                MESSAGE_CATEGORY, Qgis.Warning)
        else:
            QgsMessageLog.logMessage(
                'Task {name} completed\n'
                'Total: {total} ( with {iterations} '
                'iterations)'.format(
                    name=result['task'],
                    total=result['total'],
                    iterations=result['iterations']),
                MESSAGE_CATEGORY, Qgis.Info)
    else:
        QgsMessageLog.logMessage("Exception: {}".format(exception),
                                 MESSAGE_CATEGORY, Qgis.Critical)
        raise exception

# Create a few tasks
task1 = QgsTask.fromFunction('Waste cpu 1', doSomething,
                             on_finished=completed, wait_time=4)
task2 = QgsTask.fromFunction('Waste cpu 2', doSomething,
                             on_finished=completed, wait_time=3)
QgsApplication.taskManager().addTask(task1)
QgsApplication.taskManager().addTask(task2)


texts: Which one is better ?SQL server 2012 or Oracle 11g

texts:
I want to learn a database server, SQL server 2012 or Oracle 11g . I surf in the internet to compare SQL server and Oracle for spatial data and found good websites such as  gisboston.it was very interesting for me but i want to compare oracle and sql server about job opportunity in GIS. which one is better to find a good job?
 which one is better to manage spatial data ?


texts: Showing label if greater than percentage using Arcade with ArcGIS Pro

texts:
I have a text field called PerVax1 that has percentages.

I want to create an expression that labels the features when the field has a percentage greater than 30%. Here's some things I've tried so far:
var percentage = $feature.PerVax1;
var prcnt = "30%";
IIf(percentage >= "30%"){
return percentage}

//resulted in error:
Invalid expression.
Error on line 4.
Semicolon or new line expected.

var precentage;
percentage = IIf($feature.PerVax1 >= "30%")
return percentage;

//resulted in error:
Invalid expression.
Error on line 2.
Wrong number of arguments


texts: How can I map polygon value distribution?

texts:
I'm surveying an area and have a z-value attached to each polygon. In comparison to the area I'm studying the polygons are quite small. What I'm trying to do is show the distribution of this value. I've been researching IDW, heatmapping and Kriging but none of these seem to quite fit the purpose. The main problem I've come upon so far is by using the centroids, which aren't evenly distributed, I am analysing more that spatial distribution than the value distribution. I'm hoping to produce a raster output, purely on a visual basis as they show the spread clearly.
If anyone has any suggestions of tools, techniques or even advice on where to start looking, I would be more than happy for your help.
I'm using QGIS mainly, but also have access to ArcGIS and some other pieces of software.


texts: I would like to create a feature class using an attribute table of another feature class. However, it does not use XY coordinates, so I can't create it that way. I thought I could export the data and make it a feature class but that doesn't seem to work.


texts: Removing white background from raster as well as cleaning up edges?

texts:
I am looking to remove white background from a tiff as well as cleaning up the edges. I have used display background values from the symbology tab as well as applying gamma stretch for edges.. I can not seem to do both at the same time..

I am using ArcMap.

I have a raster that has a white background, which I can remove, but leaves jagged white/black edges.

I am trying to remove for a better visual.

It is a multiband tiff. It was derived from drone photogrammetry software.


texts: How to delete a feature in Openlayers 3.8.2 with `del` Key

texts:
I'm using OpenLayers-3.8.2 for developing a GIS application. Users draw some features and can modify all existing features on the map. They can delete features. I can delete feature as follow:

user click some features
User click on a button and then: layer.getSource().removeFeature(feature)

Now I want the user to be able to remove features with del key.

What do I implement it?


texts: Count points in polygon tool keeps failing in QGIS

texts:
I keep getting an error and when i check the logs this is what I'm seeing. All layers are using the same CRS, I've tried v.clean just in case. I'm stumped.

  Uncaught error while executing algorithm          Traceback (most recent
  call last):
              File "C:/PROGRA~1/QGIS2~1.18/apps/qgis/./python/plugins\processing\core\GeoAlgorithm.py",
  line 203, in execute
                self.processAlgorithm(progress)
              File "C:/PROGRA~1/QGIS2~1.18/apps/qgis/./python/plugins\processing\algs\qgis\PointsInPolygon.py",
  line 80, in processAlgorithm
                spatialIndex = vector.spatialindex(pointLayer)
              File "C:/PROGRA~1/QGIS2~1.18/apps/qgis/./python/plugins\processing\tools\vector.py",
  line 220, in spatialindex
                idx = QgsSpatialIndex(layer.getFeatures(request))           Exception: unknown


texts: Seeking reference/tutorial for ArcGIS Spatial Analyst?

texts:
Does anyone have a recommendation for a good reference/tutorial(s) on Spatial Analyst?

The most recent I can find are from early to mid 2000s so perhaps SA has not changed much since then.


texts: Converting data stored in Excel (grid/matrix) to a raster?

texts:
I have some grid data that is stored in a xls file. The rows are latitude and the columns are longitude values in decimal degrees (see image).

What would be a quick and efficient way to convert these data to a raster with 2° grid cells? I am using ArcMap 10.2 but also have access/experience with R or QGIS.


texts: Find nearest polygon to point using Python, OGR and Rtree

texts:
I have a dataset with points and one with polygons. When the points intersect with one or more polygons they get assigned an attribute value inhereted from the polygon, but not all points intersect with a polygon. In these latter cases I want to assign a value based on which polygon is nearest (the code block below relates to these latter cases only). Since I am dealing with over 200,000 polygons I use an Rtree spatial index. However, with the use of index.nearest() it is not always the case that the indexed polygon is the nearest when comparing euclidian distances to the point with other near features.

What I would like to do is pick the polygon that is nearest after evaluating the distance of several nearest polygons (according to the index) to a point. However, I am unsure how to translate this into the code I use now, any tips?

import sys, rtree
from osgeo import ogr

areaFile = r'C:\data\categories\classesCombined.shp'
pointFile = r'C:\data\buildings\buildingPoints.shp'

areaData = driver.Open(areaFile, 0)
areaLayer = areaData.GetLayer()
pointData = driver.Open(pointFile, 1)
pointLayer = pointData.GetLayer()

index = rtree.index.Index(interleaved = False)
for i in range(0, areaLyer.GetFeatureCount()):
    areaFeat = areaLyaer.GetFeature(i)
    areaGeom = areaFeat.GetGeometryRef()
    xmin, xmax, ymin, ymax = areaGeom.GetEnvelope()
    index.insert(i, (xmin, xmax, ymin, ymax))

for point in range(0, pointLayer.GetFeatureCount()):
    pointFeat = pointLayer.GetFeauture(point)
    if pointFeat.GetField('CLASS') is None:
        pointGeom = pointFeat.GetGeometryRef()
        xmin, xmax, ymin, ymax = pointGeom.GetEnvelope()
        for area in list(index.nearest(xmin, xmax, ymin, ymax), 5)):
            areaFeat = areaLayer.GetFeature(area)
            areaGeom = areaFeat.GetGeometryRef()
            ...

areaData.Destroy()
pointData.Destroy()


texts: Trouble with GeoPackage - is SpatiaLite a solid alternative?

texts:
I like GeoPackage because it is very versatile, small and compact and works well with QGIS. But the more I use GeoPackage and recommend it to my customers, the more negative experiences I have: especially data loss! I can't reproduce it, but my customers tell me that geometry or tables are lost or damaged while editing. Sometimes the tables still exist, but there are no entries (but the number of rows) and they cannot be deleted. I cannot repair or reproduce the problem.
Does anyone else have such experience with GeoPackage in QGIS and what could be a good alternative? SpatiaLite?
EDIT: the GeoPackage is not stored in the network or the cloud, nor is used/edit by many people at the same time, but what all GeoPackages have in common are virtual fields!


texts: Best way to build my own simple web-based map server

texts:
I've got a requirement for a project where I need some relatively simple mapping functionality on a web page.  I'm familiar with the Google Maps APIs but unfortunately this is a project that will be hosted on a closed network with no internet access.

Is there a relatively straight-forward way to host my own map server with the following functionality (preferably on Windows but I would consider Linux):

I only need to cover a map area no bigger than the size of Texas
I don't need streets, terrain, or any other complex layers. At most I would needs outlines of the counties in the region with a few text markers for large cities.
I need to be able to display custom data driven markers on the map when it is rendered (I don't care if that is on the server side or if I need to do it in javascript on the client)
2 or 3 levels of zooming would be great, but no necessary

I've tried SUAS server v3 with XAMPP and ran into numerous issues, including low-level PHP errors.


texts: Using CSV file with coordinates and multiple attributes to extract interpolated values for different set of coordinates in ArcGIS Desktop?

texts:
I have a .cvs file of 501 rows that correspond to 500 coordinate points, with each row consisting of latitude, longitude for a single location, followed by a list of about 100 weather attributes (precipitation, temp, etc. for various time periods). The first row is a header specifying all these variable. Together, the 500 points make a fairly regular grid that covers Uganda.

I would like to use these 500 points to interpolate (natural neighbor) and extract values for all 100+ weather attributes for a list of 1000 points that represent the locations of surveyed households in Uganda (also in csv format, with a household identifier). The end product that I would like to generate is a new csv (or excel file) with each row as lat, long, household identifier, and the 100 weather variables (basically the original csv but for a different set of coordinates).

Which sequence of tools can I use to accomplish this in ArcGIS? I have access to Spatial Analyst.

So far I have tried a few things but haven't been able to complete the full process with success (I only mention the steps that in hindsight seem most promising):

First, I converted the original csv of weather attributes to an excel table and then used the Excel to Table tool. This successfully generated a table that preserves the lat, long, weather attributes, and generates an ObjectID. From here, I next tried to use the Point to Raster tool (thinking I might be able to make a multiband raster from the table, but this step might not be needed), but this only allows the conversion of a single field at a time. Then I thought I would skip conversion to a raster, and directly applied the Natural Neighbor tool to the table, and this does successfully generate the interpolated raster to cover the entire geographic area of interest, but it also only allows a single field/attribute at a time.

Is there a tool that will allow me to do more than one attribute at a time, or is a loop over all the attribute names my best option?

My hope is to get the data into a form from which I can apply the Extract Multi Values to Points tool. (I am happy to use Shift + Control + C to simply copy the resulting table at the end.)


texts: If any of B contains all of A, then count # of times A is in all of B

texts:
BACKGROUND: Long time reader, first time poster. Looking for help using QGIS to score data in a CSV table that was yielded from processing the contiguous neighbors of every polygon from a file with over 750k features.

VERSION: 2.2 or 3.0

GOAL: an equation to run in Field Calculator

DATA: Column A (small string), Column B (long string), Column C (output of function, a small integer)

NATURAL LANGUAGE ARTICULATION: If any part of B contains all of A, then count all instances of all of A in all of B, output number of instances as integer to Column C

EXAMPLE

COLUMN A     COLUMN B                                              COLUMN C
---------    --------------------------------------------------    --------
Doe, John    Doe, John; Smith, John; Johnson, Harry; Doe, John;    2


texts:
I am trying to add a text field to my attribute table in ArcMap, but every time I do it adds a string type field instead.I've tried turning off the Allow Null Values and shortening the Length, but no matter what it makes it string instead of text field. If I can't make it a text field, how do I fill in the records for a string field? The table is attached to a polygon layer, I'm trying to add in short character codes to ID each polygon.
Why is it doing this and how do I get around it?


texts: Appending data from Google spreadsheet into Leaflet LineString feature properties?

texts:
How would I add several properties into a set of leaflet LineStrings from google spreadsheet via Tabletop?

I've created LineStrings from the google spreadsheet linked above and need to add the data from columns B-D into the properties of the json using the key IDs in column A, which match the IDs in the JSON string in the geometry column. I've tried a loop (and several others) that I thought would add the data into a new properties field, but it doesn't work.

I believe the issue surrounds the following:

var values = data.toArray;// thought this would grab the data
improvementsHeaders = ["code", "include", "summary"];
for (var j = 0; j < improvementsHeaders.length; j++) {
        featureCollection.properties[improvementsHeaders[j]] = values[row][improvementsHeaders.indexOf(improvementsHeaders[j])];
      };

The current script is throwing flag, Uncaught TypeError: Cannot read property '51' of undefined
Uncaught TypeError: Cannot read property '4' of undefined, which suggests there's something wrong with interplay between headers and the indexing loop.

var polygonLayer;
var pointGroupLayer;

function addPolygons(data) {
  if (polygonLayer != null) {
    polygonLayer.remove()
  }

  var featureCollection = '{"type":"FeatureCollection", "features":[';
  var delim = '';

  for (var row in data) {
    featureCollection = delim + featureCollection + data[row].geometry + ',';
    delim = '';
  }

featureCollection = featureCollection.substring(0, featureCollection.length - 1);
featureCollection = featureCollection + ']}'

var values = data.toArray;// thought this would grab the data
improvementsHeaders = ["code", "include", "summary"];
for (var j = 0; j < improvementsHeaders.length; j++) {
        featureCollection.properties[improvementsHeaders[j]] = values[row][improvementsHeaders.indexOf(improvementsHeaders[j])];
      };
//index through retrieved data and attach to features

  var geojsonStates = JSON.parse(featureCollection);

  var polygonStyle = {"color": "#2ca25f", "fillColor": "#99d8c9", "weight": 3};
  var polygonHoverStyle = {"color": "green", "fillColor": "#2ca25f", "weight": 3};

  polygonLayer = L.geoJSON(geojsonStates, {
    onEachFeature: function (feature, layer) {
      layer.on({
        mouseout: function(e) {
          e.target.setStyle(polygonStyle);
        },
        mouseover: function(e) {
          e.target.setStyle(polygonHoverStyle);
        },
        click: function(e) {
          L.DomEvent.stopPropagation(e);
        }
      });
    },
    style: polygonStyle,
  }).addTo(map);
}


texts: Polygons which are completely overlapped by second polygon - QGIS

texts: I have multipolygon like here: . And I have another, second polygon. Anyway I want to get only polygons which are completely overlapped by that second polygon. Intersection and clip in QGIS only cut the edges and that is, but like I said - I want only pieces (polygons) of this multipolygon which are completely overlapped by another.


texts:
I am using Rodeo as my Python IDE for data analysis, and I'm trying to import the ArcPy module in Rodeo. However, when I attempt to import it, I receive the error "No module named 'arcpy.'"

I think it may be an issue of Rodeo wanting to use the built-in Python, but I'd like to have it use the Python executable that ships with ArcGIS (I'm using ArcGIS 10.5).

Can I use the ArcPy module in Rodeo?


texts: How To Change Extent of An Ungeoreferenced Raster in QGIS

texts:
After photoshopping a tiff in Photoshop, the tiff loads in QGIS as shown in the photo below. Looking at the raster properties, you can see that the extent of the image was changed somehow by Photoshop. So now it loads with the top lefthand corner being placed where the bottom lefthand corner should be. One of the attached images compares the raster properties for the original and processed images. Note the difference in the extents for the two images.

Task: I would like to change the extent data in the processed image so that it is the same as the original image. I have not been able to find a method to this via Googling.


texts: Converting Spatial Points to Neighbours List using R?

texts:
I am trying to convert some spatial points into a listw object to run the moran.mc function on my data. The data are not adjacent but the points cluster in spatial patterns. Here is how I have tried to write the code:

library(sp)
library(raster)
library(spdep)

latitude <- c(56.33, 48.12, 51.13, 60.83, 52.46, 60.42, 51.07, 52.02, 52.87, , 55.27)

longitude <- c(84.58, 88.35, -55.93, -101.55, -116.19, -150.90, -89.90,
-90.13, -89.93, -88.93)

core_locations <- cbind(longitude, latitude)

pts <- *SpatialPoints(core_locations)

I have tried two approaches to converting this into the listw; the first is to create a weights matrix, but I can't figure out how to pass this to moran.mc

dist <- pointDistance(pts, lonlat = TRUE)

weights_matrix <- 1 / dist
weights_matrix[!is.finite(weights_matrix)] <- NA
rtot <- rowSums(weights_matrix, na.rm = TRUE)
weights_matrix <- weights_matrix / rtot

The second is try to convert the points into polygons, but maybe this is just a bad idea....

p <- Polygon(pts)
ps <- Polygons(list(p),1)
sps <- SpatialPolygons(list(ps))
w <- poly2nb(sps)
ww <-  nb2listw(w)

*I have a spatial points data frame for these locations with associated attribute data for each point.


texts:
Is it possible to go to a specified row number (or perhaps the next/previous) row using PyQGIS?

This code gets the value of the first attribute in the first row of the table

features = aLayer.getFeatures()
attrs = feat.attributes()
attribute0 = attrs[0]
print attribute0

How can I get the value of the same column (0), for row 2, 22, 222...? Is there a way to do this without using a for loop?


texts: Ways to recover all files in GeoPackage that were accidentally overwritten in QGIS 3.4?

texts:
I was saving a new layer in my GeoPackage when a window showed that I was going to overwrite a file. I chose the "Overwrite File" button because it usually overwrites a single file, not the whole GeoPackage, in my experience.

How can I recover this? The whole working folder is backed up to Google Drive but the main gpkg-file was not there, only the shm and wal. The syncing was off for the whole weekend until the mishap so the main gpkg-file should be there.


texts:
I want to create a valid Shape object in shapely which has properties but has no geometry.

e.g. For creating features to be saved into a shapefile. All features have attributes/properties but some have null/empty geometry.

I am currently using the workaround of creating an empty GeometryCollection.

s = shapely.geometry.collection.GeometryCollection()
s.properties = {"prop1":"val1"}

Is there a more suitable/valid approach or a better shapely geometry class for storing null/empty geometry?


texts: Applying Raster Calculation to Image Collection in Google Earth Engine

texts:
I am trying to perform a raster calculation (EVI) on an imagery collection.  This operation is seen on Line 130 and is based on:
https://developers.google.com/earth-engine/image_math (which I know is for single images).  In addition to being an imageCollection, my version of this GEE provided code differs in that I've renamed my bands prior to calculating EVI.
I am having trouble finding the necessary format to apply a band math equation to an imagery collection.
Presently, I am receiving an error:
"Cloudfiltered.expression is not a function"
Here is my script:
https://code.earthengine.google.com/0e9ae4a031ddd2d46fe44d40cc03b05e
The Region of Interest: "ROI" was made with the Rectangle Tool in GEE and can be totally arbitrary; I don't know how to share it. Here is mine for context: https://imgur.com/a/86H9iJ1


texts: Adjusting Override Settings (Cellsize) for Optimized Hotspot Analysis using ArcGIS Pro

texts:
I am presently working on a component of a project that involves utilizing the Optimized Hotspot Analysis (i.e. GI* Statistic) tool in ArcGIS Pro. Specifically, I am investigating tornadogenesis hotspots in the states of Arkansas and Oklahoma and determining potentially statistically significant areas (i.e. hotspots and coldspots) of tornado activity.
I have my genesis point features as my input and my incident aggregation method chosen as count incidents by hexagon grid. Furthermore, I have set my bounding polygon as the Arkansas and Oklahoma. The hotspot analyses are being executed separate of eachother (i.e one for AR; followed by one for OK)
By default, ArcGIS computes the appropriate cell size. It is evident that it changes the output (having experimented with the settings).
How might adjusting the cell size in the override settings impact the overall output?


texts: Making two pen style in one vector file in QGIS?

texts:
I am working with QGIS and I have a vector file of the earth Latitude and longitude .
I need to show the horizontal lines as a dash lines and the vertical as a solid line. I am sure that if I make to layer, one horizontal and vertical, it will solve me problem.

However, I wonder is there a more elegant way to solve this problem?

Can it be done with the Attribute Table?


texts:
when I use gdawarp command it was failed - reproject with proj4. It was related with msvcr120.dll.
Surely I was already install  MS VC++2013 redistributable (x64). I was copy dll to gdal folder but same.
I can't find any questions related with msvcr120.dll and gdalwarp in history.
Has anyone encountered this problem?

windows 7 pro, x64

I was installed 'gdal-201-1800-x64-core'  and I was set system variables GDAL_DATA, GDAL_DRIVER_PATH

My command is:

gdalwarp -of GTiff -co TFW=YES -t_srs "+proj=tmerc +lat_0=38 +lon_0=127 +k=1 +x_0=200000 +y_0=600000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs" .\data\a102113.tif .\data\b5186.tif


texts: Recovering shapefile with green question marks?

texts:
I have a shapefile with a green question mark in the ArcCatalog.
The shp-extension in the shapefile has 0Kb size

Is there any way to recover the corrupted file?


texts: Splitting feature classes into different geodatabases using ModelBuilder

texts:
Like another post on here - I'm using ArcGIS Desktop 10.6.1 (ModelBuilder) to split feature classes in a geodatabase by an attribute field ("Country").  This process is fine but the issue is exporting these individual feature classes to separate geodatabases (using Feature Class Iterator).

I want to copy these split "country" feature classes, copy these separate features to a geodatabase with the same "Country" name - the result to be separate country geodatabases with a single country inside.

So far with the site model below - I can make a split the site feature classes by country into a geodatabase and can create separate "Country" geodatabases using the iterator - but fail to copy the feature classes into these.

Can someone provide me step-by-step directions to solve this issue?


texts: Maps hosted on my own server instead of ArcGIS servers?

texts:
I know ArcGIS provides maps for me to use off of their server for the price of their services. http://resources.arcgis.com/en/help/main/10.1/index.html#//0154000003vt000000

I want to share my map with other clients without requiring them to buy a license or anything.

Is it possible for me to use my own map hosted on my own server, maybe something like google maps for free?

The closest option I have seen is OpenStreetMap http://www.openstreetmap.org/. I want to add all the same layers as ArcGIS such as legend, shapes and editing.


texts: Best database for Geographical database for my case

texts:
I'd like help to make a decision on database change.

For now my case is:

MySQL on Cloud SQL (GCP): Storing GPS Tracking Data for now is 400Million Record (Around 400GB). And going bigger every day. For now is quite slow to processing the report like routing report, trip report

I have plan to move all of these to Google Cloud Datastore.

Is this is a good case for my situation to move to Google Cloud Datastore?
The reason why i'm choosing Google Cloud Datastore is:
Easy Scaling, Faster query.

But i had wondering that Google Cloud Datastore can write only 1 operation/second. I need to write the GPS Tracking every second for more than 10K Device

(I'm sure that i'm not understand this clear so can anyone explain me about transaction or the way can write more entity per second)

Does anyone have anything to recommend to let me know the best solution?


texts:
Using QGIS and other free tools, how can I use a Google Maps route to create a route animation in the style of the Indiana Jones movies (https://www.youtube.com/watch?v=5TY5Fp6O5iM)?

I have found two questions on this site (Show vehicle traveled route as animation and Interpolate points between coordinates for smooth animation in google maps or openmap) which deal with similar issues, but both of them are based on Google Maps or OpenLayers integrations and require coding.

I have found a very useful article on how to use the Time Manager QGIS plugin to map trajectories which fade out. In my case I would not need the path to fade over time, but that method could provide a useful basis for getting the laser effect that I want.

A further complication is that while I know the route that I want to animate, I do not have any actual data to use as inputs for the animation. I could get a route out of Google Maps, but that does not include timestamps so I'd either need a way to animate the steps in order, or otherwise add a timestamp somehow.

How can I put all of this together, or is there perhaps a more elegant way to achieve this effect using tools that are freely available?


texts: Querying Thousands of Points with ST_Value()?

texts:
Currently I have an application using the psychopg2 Python library to query a database for elevation data. My current python implementation looks like the following:

   def GetElevation(lat, lon, cur):
    point = "'SRID=4326;POINT({} {})'::geometry".format(lat,lon)
    cur.execute("SELECT ST_Value(rast, {}) FROM dted0 WHERE ST_Intersects(rast, {});".format(point, point))
    return cur.fetchone()[0]

This works, but I am was curious if I could pass in an array of latitutdes and an array of longitudes. I tried building a "point string" with several thousand queries, but I get an error saying that I can only pass 100 arguments to the function.

Is there a way to pass all the points I need in a single query?

The reason I am doing this is because I assumed doing only one transaction with the database would be faster rather than looping for each query.

-------------------EDIT 1----------------------

In an attempt to implement the top rated post, I've done the following:

import json
json_arr = []
for lat, lon in zip(lats, lons):
    json_arr.append({'lat:': lat, 'lon':lon})
json_str = json.dumps(json_arr)

Next I defined a function per his advice:

def GetElevationsBatch(json_elevations, cur):
    cur.execute(
    """SELECT ST_Value(rast,z.point)
       FROM dted0
       JOIN (
          SELECT ST_SetSRID(ST_MakePoint(lat, lon), 4326) as point
          FROM json_to_recordset(%s) AS z(lon double precision, lat double precision)
       ) AS z
       ON ST_Intersects(rast,z.point)""", (json_elevations,))
    return cur.fetchall()

However, when I call the function I get I don't get good things:

a = GetElevationsBatch(json_str, cur)
print("a = ", a)
# Result: a =  None

-----------------------EDIT 2--------------------------

I've left the python plugin out for now until I can find the best query, so here is my latest attempt:

SELECT (ST_Dump(gv.geom)).geom, gv.val
FROM srtm , LATERAL ST_Intersection(rast, 'SRID=4326;MULTIPOINT(-111.305048568 38.0601633931,-111.822991286 38.6025320796,-111.136977796 38.3631992596,-111.206470006 38.971228396)'::geometry)  AS gv
WHERE ST_Intersects(rast, 'SRID=4326;MULTIPOINT(-111.305048568 38.0601633931,-111.822991286 38.6025320796,-111.136977796 38.3631992596,-111.206470006 38.971228396)'::geometry);

In the above query, I am attempting to use the notion of MULTIPOINT geometries. However, I found that these calls "work", but they are much slower than my original method of simply querying for every point. I understand that I am not running on any amazing hardware, but should calls really take on the order of seconds for a simple elevation query? Seems to me there is something awry here. It is taking nearly 3 seconds to retrieve 4 points with the above code. To compare to my original method, it only takes about 100ms to retrieve 4 points.

Shouldn't a solution in which I only query the database once be quick than one where I have to query several times?


texts: NAD 1983 to NAD 1983 CSRS Transformation

texts:
I have data in NAD 1983 Zone 20 and want to transform it to NAD 1983 CSRS Zone 20.  I found this source: https://community.esri.com/thread/6772 which states:

NAD_1983_To_NAD_1983_CSRS_4  > 1950  Canada - Newfoundland, New Brunswick, & Nova Scotia offshore

Since I am not native in English, the word "offshore" bothers me. I hope somebody could confirm if this is the best transformation for data located in the Nova Scotia mainland or suggest a better transformation?


texts:
I have a number of MXD files with data sources linking to an Enterprise geodatabase connection string.
I would like to update these connection strings with ArcPy programmatically to a new data source, much the same as if I was using "Set Data Sources" in ArcCatalog manually.

This works when I hardcode the current data source like this:

mxd.findAndReplaceWorkspacePaths(r"path/to/sde", path2newcon, False)

However, it would be better if I do not have to hard code the current connection but get it from the MXD file with ArcPy.

I can get the sources of each layer like this:

for lyr in arcpy.mapping.ListLayers(mxd):
    if lyr.supports("DATASOURCE"):
        print "Layer: " + lyr.name + "  Source: " + lyr.dataSource

But I would like to get the connection string for the whole MXD file without any layer information.

How can I get the current data source with ArcPy?