texts: Error during a for loop in Python terminal from QGIS 2.18.17

texts:
When I run the script below in Python Terminal at QGIS 2.18.17 (Linux Mint 19) it always crash after doing 73 iterations at the 2nd for loop from the script. Could anyone explain why this happens and could suggest how can I fix it?

For now, the way that I found to deal with this problem is to close QGIS, open it and run the script again from the item that it stopped. But when you have more than 2000 itens in a list of a for loop this can be quite annoying.

import processing
import os
import pandas as pd
wd = "/home/tiago/Documentos/PRIM_Mineracao/mata_atlantica/up_centroides/split"
lc = list()
for file in os.listdir(wd):
    if file.endswith(".shp"):
        lc.append (wd+"/"+file)
drenagem = QgsVectorLayer('/home/tiago/Documentos/PRIM_Mineracao/mata_atlantica/redes_drenagem/rio_7.shp', '','ogr')
fim = QgsVectorLayer('/home/tiago/Documentos/PRIM_Mineracao/mata_atlantica/fozes/foz_1.shp', '','ogr')
ext = drenagem.extent()
xmin = ext.xMinimum()
xmax = ext.xMaximum()
ymin = ext.yMinimum()
ymax = ext.yMaximum()
coords = "%f,%f,%f,%f" %(xmin, xmax, ymin, ymax)
ups = "/home/tiago/Documentos/PRIM_Mineracao/mata_atlantica/up/up_mata.shp"
pasta_intermediarios = "/home/tiago/Documentos/PRIM_Mineracao/mata_atlantica/v_net_distance"
p = 1
## 2nd FOR LOOP ##
for pto in lc:
    nome = os.path.basename(pto)
    nome = nome.replace(".shp","")
    nome = nome.replace("up_centroide_diss_COD_","")
    v_dist= pasta_intermediarios+"/v_dist_"+nome+".shp"
    disol= pasta_intermediarios+"/v_disol_"+nome+".shp"
    explo= pasta_intermediarios+"/v_explo_"+nome+".shp"
    points= pasta_intermediarios+"/v_points_"+nome+".shp"
    inter = pasta_intermediarios+"/v_inter_"+nome+".shp"
    processing.runalg("grass7:v.net.distance",drenagem,pto,fim,50,2,"","",0,"","",None,None,None,False,coords,-1,0.0001,0,v_dist)
    processing.runalg("qgis:dissolve",v_dist,True,None,disol)
    processing.runalg("qgis:explodelines",disol,explo)
    processing.runalg("saga:convertlinestopoints",explo,True,0.004,points)
    processing.runalg("qgis:intersection",points,ups,False,inter)
    inter_shp = QgsVectorLayer(inter,'',"ogr")
    list_up = list()
    for feature in inter_shp.getFeatures():
        up = int(feature["UP_UC"])
        if up not in list_up:
            list_up.append (up)
    jusante = list_up[:]
    del jusante[0]
    jusante.append(0)
    df = pd.DataFrame(jusante,list_up)
    nome = os.path.basename(pto)
    nome = nome.replace(".shp",".csv")
    nome = nome.replace("up_centroide_diss_COD_","")
    df.to_csv('/home/tiago/Documentos/PRIM_Mineracao/mata_atlantica/tabelas_jusante/'+nome)
    path = os.path.dirname(v_dist)
    files_to_del = []
    for r, d, f in os.walk(path):
        for file in f:
            files_to_del.append(os.path.join(r, file))
    len(files_to_del)
    for file in files_to_del:
        os.remove(file)
    os.remove(pto)
    print p, "de", len(lc)
    print nome,"done"
    p = p+1

This is the message that I get.

Unable to execute algorithm
Wrong parameter value: /home/tiago/Documentos/PRIM_Mineracao/mata_atlantica/v_net_distance/v_dist_UP_764949.shp
Unable to execute algorithm
Wrong parameter value: /home/tiago/Documentos/PRIM_Mineracao/mata_atlantica/v_net_distance/v_disol_UP_764949.shp
Unable to execute algorithm
Wrong parameter value: /home/tiago/Documentos/PRIM_Mineracao/mata_atlantica/v_net_distance/v_explo_UP_764949.shp
Unable to execute algorithm
Wrong parameter value: /home/tiago/Documentos/PRIM_Mineracao/mata_atlantica/v_net_distance/v_points_UP_764949.shp
Traceback (most recent call last):
  File "<input>", line 1, in <module>
  File "/tmp/tmp2K40DP.py", line 42, in <module>
    del jusante[0]
IndexError: list assignment index out of range


texts: Catching WMS error message from log messages panel in QGIS with python?

texts:
I am using a WMS-service as background map in an automatized atlas print process.

Unfortunately the WMS-service seems to have a problem on server side, as it is sometimes not responding (see error message: map request failed ... server replied bad request (see attached image)). By just triggering a repaint after some seconds or minutes, the WMS loads correctly.

Now I would like to catch this behaviour via python/pyqgis because the atlas print process produces empty background pages when the WMS is not responding. So I would need to check if this particular error occurs, and force a triggerRepaint() until the WMS-server responds correctly.
I don´t know how, since the error appears only in the log message panel, under the tab "WMS".
I already tried with try:...except... or by checking wmsLayer.isValid(), but this did not do the trick...

Is there anybody who generaly knows a way to catch this particular error appearing in the log messages panel?

Unfortunately the WMS-service is licensed, so I´m afraid the problem might not be reproducable...But I hope there is anybody out there who knows a way to catch errors from the log message panel.

I am using QGIS 2.18.2

With advice from Germán I was able to catch the error in the following way:

# coding=utf-8

from qgis.core import *

wmsLayer_name="wms-dtk50_wgs"
url_with_params ='url=http://sg.geodatenzentrum.de/wms_dtk50?&crs=EPSG:25832&featureCount=10&format=image/png&layers=DTK50&styles='

wmsLayer = QgsRasterLayer(url_with_params, wmsLayer_name,'wms')
QgsMapLayerRegistry.instance().addMapLayer(wmsLayer)

def errorCatcher( msg, tag, level ):
    if tag == 'WMS' and level != 0:     #Warnings or Errors (0: Info, 1:Warning, 2:Error)
        print "WMS error detected!"
        myWMSLayer = QgsMapLayerRegistry.instance().mapLayersByName("wms-dtk50_wgs")[0]
        myWMSLayer.triggerRepaint()

# connect with messageReceived SIGNAL from QgsMessageLog to an errorCatcher custom function
# instantly reacts if error/warning occurs
QgsMessageLog.instance().messageReceived.connect( errorCatcher )

#after 100 times triggering a "wmsLayer.triggerRepaint()",
# I get following warning in log messages panel "WMS":
# "2017-01-17T07:17:52  1   Not logging more than 100 request errors."

Obviously the "WMS" provider seems to have a restriction of sending error requests to the messages log...opened a new question with this issue here: How to solve issue with log messages panel in QGIS: “Not logging more than 100 request errors.”?


texts: How do you automatically number cells in a grid with row numbers and column numbers (or letters)

texts:
I have a fishnet grid consisting of 330 rows and 330 columns.  Each cell is 500 feet x 500 feet.  The centroids of each cell is 500 feet from its adjacent centroid throughout the grid.  Also, the centroid points have been extracted for each polygon (i.e. cell) as well as adding an X and Y field, and calculating the geometry to get those values.
Now comes the fun part.  I would like to find a way to automatically number each row (top to bottom) and each column (left to right).  In other words, first row would all have the number "1" in each cell, the second row would all have the number "2" in each row, etc. without having to resort to selecting each row and Field Calculate a "1", "2", "3", etc. for 330 rows or "A", "B", "C", etc. for 330 columns.


texts: Using ui.select for Administrative levels dropdown in Google Earth Engine

texts:
How can I use ui.select for Administrative levels dropdown in Google Earth Engine? I would like to use ui.select to first select the state and then select the county.

Here is the GEE code:

var countries = ee.FeatureCollection("TIGER/2016/States")
var Counties = ee.FeatureCollection("TIGER/2016/Counties")

var countryNames = countries.aggregate_array('NAME')
var countyNames = Counties.aggregate_array('NAME')

var countryNames = countries.aggregate_array('NAME').getInfo()
function countryNames(names) {      return names}

var countyNames = Counties.aggregate_array('NAME').getInfo()
function countyNames(names) {      return names}

var state = {
country: ui.Select({
items: countryNames,
}),

county: ui.Select({
items: countyNames,
})
};

var update = function(NAME) {
var roi = Map.layers().get(0)
[countries.filterMetadata("NAME", "equals", NAME).getValue()]
[counties.filterMetadata("NAME", "equals", NAME).getValue()]

}

state.country.onChange(update);

state.county.onChange(update);

ui.root.insert(0, ui.Panel([state.country, state.county]));

https://code.earthengine.google.com/7d0286ea2eaf4de2a1a8915b9c0d71d6


texts: Buffer does not intersect or contain the point used for the center of the buffer when creating the buffer

texts:
I made the buffer below (see code).  I then tried to see if the buffer contains or intersects the center of the buffer, but False is returned.  Why?

from shapely.geometry import Point, Polygon

import pyproj

srcProj = pyproj.Proj(init="epsg:5070")
dstProj = pyproj.Proj(init='epsg:32612')

pt = (-94.6091, 39.1063)

x, y = pyproj.transform(srcProj, dstProj, pt[0], pt[1])

pt = Point(x,y)

buf = pt.buffer(100)

buffer_points = zip(buf.exterior.coords.xy[0],
buf.exterior.coords.xy[1])

proj_buffer_points = []

for point in buffer_points:
    x = point[0]
    y = point[1]
    x, y = pyproj.transform(dstProj, srcProj, x, y)
    proj_buffer_points.append((x, y))

print Polygon(proj_buffer_points).intersects(pt)
print Polygon(proj_buffer_points).contains(pt)

False


texts: How I can merge three rasters with sum cell

texts:
I am working by ArcGIS 10 and I have three raster
all of raster have value cell
and I want output raster have the sum of pixels values


texts: Stopping Torque from looping in CartoDB editor

texts:
I'm working on a torque map and would like it to just animate once and then stop. I've found loop=false, but am unsure if I can apply this in any way when I'm using the CartoDB editor.

Does anyone know if that is possible?


texts: Cache not generating in ArcGIS Server v10.1

texts:
Apologies for any initial ambiguity in the question, I will attempt to make it a little clearer.

I have published a map service in ArcGIS Server.
In publishing, I have specified this map service to be a cached map service (I have set the levels, cache location, tile format, etc. AGS has no problem with any of these settings).
The service publishes successfully and AGS tells me that I can monitor the  creation of the cache through ArcCatalog.
When I go to ArcCatalog and check on the Cache Status it tells me that the cache generation is not in progress, (see screen grabs).

The AGS log files are telling me:

  SEVERE Instance of the service 'System/CachingTools.GPServer' crashed.
  2014-06-12T10:03:25 Server ‘Server Name’  8252  11464 1

  SEVERE Instance of service 'Service Name.MapServer' failed to process
  a request. AutomationException: Unspecified error  2014-06-12T10:03:24
  Server Name  8001  3161900 17

  SEVERE IPersistStream::Save failed.  2014-06-12T10:03:24 Service
  Name.MapServer

Has anyone else experienced this, and if so, have you figured out what the issue is?


texts: How can I reproject a shp file in NAD_1983_UTM_Zone_15N to something more standard

texts:
I have a public shapefile data set which is in a format that doesn't contain actual latitudes and longitudes (prj file below).  I'd like to load it into SQL and be able to utilize the latitude and longitude of the geometry directly, but I don't know how to reproject the shapefile.  I don't have ArcGIS.  I sounds like ogr2ogr is the best route to do the reprojection.  However, I don't know what parameters to enter to a) tell it what to convert from (based on the prj file below) or b) what I should be converting it to (what is a good standard?).

Here is all the information from the prj file:

  PROJCS["NAD_1983_UTM_Zone_15N",GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-93.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]


texts:
Not sure what I am doing wrong. Keep getting this error:

Exception has occurred: Execute Error ERROR 000622: Parameters are not
valid ERROR 000628: Cannot set input into parameter Coordinate_System

Code:
import arcpy
from arcpy import env

env.workspace = r"C:\PSP\WeeklyWork\Week10\ArcPyAutomation\AutomationData\CANADA"

in_features = "lakes.shp"

properties = "PERIMETER_LENGTH_GEODESIC"

length_unit = "KILOMETERS"

coordinate_system = arcpy.SpatialReference("Projected Coordinate
System/Continental/North America/Canada_Lambert_Conformal_Conic")

arcpy.AddGeometryAttributes_management(in_features, properties, length_unit, "", coordinate_system)


texts: Finding values that are different between two tables in MapInfo/MapBasic based on a unique ID

texts:
I am comparing 2 tables which are mostly the same but there have been some amendments, in effect comparing a new version to the old one and displaying all the records with changes in a selection table based on a unique ID.

I run into problems with the following SQL part of the script:

Select * From selectedTableNew, selectedTableOld
Where (selectedColumnNew = selectedColumnOld)
And (aColNew <> aColOld)
Into AmendedRecords

I receive the error

  No join specified between selectedTableNew & selectedTableOld. Invalid join condition in WHERE clause."


texts:
I've been trying to import all the ways for Romania into my PostgreSQL database using Geofabrik and osm2po.
However, I noticed that some chunks of a county road are missing by using the following query:
SELECT id, geom_way FROM ro_2po_4pgr,
(
    SELECT ST_Expand(ST_Extent(ST_SetSRID(ST_MakePoint(23.2161769, 45.6347162), 4326)), 0.2) as box
) as bbox
WHERE ro_2po_4pgr.geom_way && bbox.box

Which resulted in:

The missing chunks for road DJ705A that I have found on openstreetmap.org are:
254048490
258364773
661526471
258367634
254112343
778597709
778324715
You can find them by replacing the last param at this url: https://www.openstreetmap.org/way/254048490
With the help of a gentleman from Geofabrik (I was thinking that there was the issue) I received the confirmation that these ways do exist, plus he gave me the following command to run for myself:
osmium cat -tway romania-latest.osm.pbf -fopl | grep -e w254048490 -e w258364773 -e w661526471 -e w258367634 -e w254112343 -e w778597709 -e w778324715

Indeed, these ways do exist in the .pbf file, but they are not present in the generated .sql file.
I tried increasing the memory, the tile size, even importing service and living_street highway types.
I tried parsing the data from .bz2 with the comand:
curl -s -L http://download.geofabrik.de/europe/romania-latest.osm.bz2 | bzcat -c | java -Xmx12g -jar osm2po-core-5.3.2-signed.jar prefix=ro cmd=c tileSize=x postp.0.pipeOut=true log.0.to=err postp.0.class=de.cm.osm2po.plugins.postp.PgRoutingWriter > ro.sql

(here it spits the output straight to an .sql file)
... it still did not work.
This is the command that retrieved the maximum amount of ways (740073)
java -Xmx32g -jar osm2po-core-5.3.2-signed.jar cmd=tjsp prefix=ro tileSize=x http://download.geofabrik.de/europe/romania-latest.osm.pbf postp.0.class=de.cm.osm2po.plugins.postp.PgRoutingWriter

(whether I used 6GB, 10GB, 12GB or 32GB it did not make any difference)
I could really use some help, since I ran out of options and all of this is new to me.


texts: Select stream network within a fixed distance from a point

texts:
Is there a tool out there that would allow me to select a portion of stream network that is within a fixed distance from a specific point on that network?

In the example below, there is an acclimation site on Palmer Creek (green dot).  I want to be able to select all portions of my stream network that are within 4 steam miles (not euclidean) of that point.  In the example below I manually selected the stream segments within 4 miles.

Is there a tool to automate this process, as I have upward of 25 points?


texts: Using UNC pathname with Calculate Field tool in ArcPy?

texts:
This question is different to Using Calculate Field to fix hyperlinks in ArcPy gives Cannot use VB for services?
which was asked previously because it is using the Python parser specifically, rather than VB.

Syntax calculates the field correctly within tool (see raw below).  Copy the python snippet directly from the results window to migrate into main script, and the syntax returns as invalid.  Not sure how to fix.

Raw (working):

"\\\\ftp\\raw\\" + "20160206" + "\\Win\\Charles\\" + !NAME! + ".shp"

Copied snippet verbatim (Invalid):

arcpy.CalculateField_management(outShp, "Link", """"////ftp//raw//" + "20160206" + "//Win//Charles//" + !NAME! + ".shp"""", "PYTHON_9.3","#")

Having trouble getting it to run properly, have tried adding 'r' and using different combinations of escapes and it just won't run.


texts:
How do I find duplicates in a certain field?

Plus, for every duplicate, I would like to obtain corresponding row value from another field.

I am new to python programming.


texts:
I'm new at making scripts and I´m having some issues in my first arcpy script.

I've created a code that automates maps creation, and it's running ok but gradually increases the memory consumption until the memory become full, the process stops and arcmap exhibits an "out of memory" message.

I've used the arcpy.Delete_management() to try to clean up the memory, but is not affecting so much the memory cleaning...

Is there another way to clean the memory?

The link to the code is posted below...

I posted the link because the code is too large to paste here in gis.stackexchange

https://drive.google.com/file/d/0B7468mS6n7QJUHdLUUl1Ym5lcUk/view?usp=sharing

Please, remember that I´m a scripting rookie


texts: directions (path find) given kml supplied bounds

texts:
Im trying to generate directions using google.maps.DirectionsService
and I really need those directions bounded by the paths defined from
within a pre-made kml file. I have reviewed the API docs for google maps and doesnt seem they have this functionality.

This is for a custom app where areas are
already defined inside a kml file that is layered over google maps.

Any ideas?

Thanks


texts: Splitting a Polygon by a MultiPolygon is unsupported. Both geom columns are MultiPolygon

texts:
I'm trying To st_split() on two tables. the output is like this:
Splitting a Polygon by a MultiPolygon is unsupported
 . Although both fields are Multipolygons. Created new columns where the geometry type is not specified. try:

create table processing.split_20 as (
  select
  inventory_elements.ogc_fid as inv_id, seg_20.ogc_fid as seg_id,
  ST_Split(inventory_elements.geom, seg_20.geom)
  from input_d.inventory_elements, input_d.seg_20
  )

-nothing changes - Splitting a Polygon by a MultiPolygon is unsupported and the types have not changed again both multipolygons

although after creating the columns and set the data I initiated alter to geometrii

ALTER TABLE input_d.inventory_elements ALTER COLUMN geom type geometry(Geometry, 4326)


texts: Counting number of points of interests within polygon in ArcGIS ModelBuilder?

texts:
I already checked all similar questions concerning this topic,
but somehow it's not working for me.

I use ArcMap 10.4.1 and want to build a script with ModelBuilder to count the number of points within my polygon. In addition the points display differnt POI'S (each poi has a unique ID) and the polygon are districts (each district with a unique ID).

Now first of all I want to count just the quantity of points within each district.

I already tried it with -

Intersect tool (output type = POINT)
Summary statistics (COUNT FID of Polygon)
Join Field

But all I get is 0 as a quantity ...

I checked it and it seems that until summary statistics everything works as it should, but as soon I want to join the fields there is the problem that some parts of the polygon don't have any points within and I guess arcgis has some issues with that - is there a way where I can still join my table to the attributes of the polygon shape?


texts: Is there a Field Calculator or GUI based way to move points by angle and distance in QGIS?

texts:
In QGIS what is the best way of moving an existing point (or digitising a new point) when you have:

the starting x,y (UTM) coordinate;
the angle from compass bearing towards the unknown point, and;
a distance to the new point?

Hopefully this could be achieved using the field calculator by adding specific formulas in a new field, the QGIS GUI or a pre existing plugin?

(Would like to do this for singular points as well as for a list of coordinates or point shapefiles)

QGIS 2.16.2
Using UTM coordinates in Australia

EDIT:
Ideally the solution would be able to be incorporated into/accessed by or calculated on a point shapefile in the QGIS graphical modeler to work with other algorithms.

(NB: the usual originating points will be coming from someone's Garmin GPS ... so my initial steps are to smooth that transition in the best way possible to convert to .shp (converting to UTM format) as well as adding the bearing and distance fields from field notes but hopefully I can also automate that in to the data capture stage.)

Here is some sample data that would be worked on
Shapefile:
http://www.filedropper.com/pointstomoveshpfile

Spreadsheet:
http://www.filedropper.com/pointstomove


texts: Why does my python script run out of memory from IDLE but runs Ok from Python Window in ArcMap?

texts:
Why does my code work in ArcMap but not in IDLE?

It is similar to:
ArcPy function works in ArcMap Python window but not from standalone script?
except that the answer to that question does not apply to me.

I also am not just running a function from my script in the python window. I am running the exact same script.

There is also Script runs in ArcMap Python window but not in standalone PythonWin or IDLE? but no great answer there.

I suspect it has something to do with Geoprocessing Options (my specific hunch is "background processing"?) that may not apply when I'm using IDLE (see last screenshot) but I'd really like to know what specifically causes this so I can run scripts like this from IDLE (and not from ArcMap every time). Is there a way to set options like this in the script? I am using 10.3.1 so I assumed what I've read about needing to install background processing no longer applies.

I wrote this script to process contours from a series of smoothed DEMs. I was surprised when it crashed when running it in IDLE after successfully testing on a small subsample. On the full-sized project data, it crashed before successfully creating the first set of contours (1/17) and the result was unfinished (big gaps of contour lines missing).

Just for grins, I turned to the Python window in an open ArcMap session and ran it with the "execfile" command and lo and behold it is running as I write this up (appears now to be almost done). This screenshot shows that it crashed at the shell and also shows two successful sets of contours in the data frame in the background.

Here is the script:

import arcpy
from arcpy import env
from arcpy.sa import *

# Check out the extension
arcpy.CheckOutExtension("Spatial")

contourInterval = 2
baseContour = 0
#inputSmthDEM = r"L:\2014\FULL_DEM\WORK\garfmn9"
inWorkspace = r"L:\2014\FULL_DEM\WORK\Contours.gdb"
outWorkspace = r"L:\2014\FULL_DEM\WORK\Contours.gdb\CONTOURS"

env.workspace = inWorkspace
rds = arcpy.ListRasters()

contourBaseName = "contours_"
for rd in rds:
    clen = len(rd)
    #cnew = clen - 7 # minus 7 because base name output_ has 7 characters in it
    # MODIFICATION
    cnew = clen - 4 # minus 4 because FMN_ has 4 characters in it or minus 2 because G_ has 2 characters etc.
    # print contourBaseName + str(fc[-cnew:][:-4])
    outputName = contourBaseName + str(rd[-cnew:]) # [:-4]) # add a "minus 4" if input are shapefiles but these are rasters
    Contour(rd, outWorkspace + "\\" + outputName, contourInterval, baseContour)

and here is the error:

    Traceback (most recent call last):
  File "C:\dev\batchContours.py", line 25, in <module>
    Contour(rd, outWorkspace + "\\" + outputName, contourInterval, baseContour)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.3\ArcPy\arcpy\sa\Functions.py", line 5722, in Contour
    z_factor)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.3\ArcPy\arcpy\sa\Utils.py", line 53, in swapper
    result = wrapper(*args, **kwargs)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.3\ArcPy\arcpy\sa\Functions.py", line 5715, in Wrapper
    z_factor)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.3\ArcPy\arcpy\geoprocessing\_base.py", line 504, in <lambda>
    return lambda *args: val(*gp_fixargs(args, True))
ExecuteError: ERROR 999999: Error executing function.
Out of memory [contours_1]
ERROR 010005: Unable to allocate memory.
ERROR 010067: Error in executing grid expression.
Failed to execute (Contour).

My Geoprocessing Options:


texts: QGIS Plugin to add point to map without creating a file

texts: I am looking for a QGIS Plugin to add a point to the map by entering it's coordinates. The point should only be there temporarily (as long as the QGIs project is open) and not be saved in a file.


texts:
I have added this add on with g.extension on my grass 6.4.3RC3. But it doesnt work... this is the error I can see:

Traceback (most recent call last): File "C:\Program Files (x86)\GRASS GIS 6.4.3RC3\etc\wxpython\gui_core\prompt.py", line 839, in OnItemSelected

self.cmdDesc = gtask.parse_interface(GetRealCmd(cmd)) File "C:\Program Files (x86)\GRASS GIS 6.4.3RC3\etc\python\grass\script\task.py", line 470, in parse_interface

tree = etree.fromstring(get_interface_description(name)) File "C:\Program Files (x86)\GRASS GIS 6.4.3RC3\etc\python\grass\script\task.py", line 448, in get_interface_description

"\n\nDetails: %(det)s") % { 'cmd' : cmd, 'det' : decode(cmderr) } grass.script.core . ScriptError : Unable to fetch interface description for command 'r.basin.py.bat'. Details: C:/Program Files (x86)/GRASS GIS 6.4.3RC3/scripts/r.basin.py: line 74: import: command not found C:/Program Files (x86)/GRASS GIS 6.4.3RC3/scripts/r.basin.py: line 75: import: command not found C:/Program Files (x86)/GRASS GIS 6.4.3RC3/scripts/r.basin.py: line 76: import: command not found C:/Program Files (x86)/GRASS GIS 6.4.3RC3/scripts/r.basin.py: line 77: import: command not found C:/Program Files (x86)/GRASS GIS 6.4.3RC3/scripts/r.basin.py: line 78: from: command not found C:/Program Files (x86)/GRASS GIS 6.4.3RC3/scripts/r.basin.py: line 79: from: command not found C:/Program Files (x86)/GRASS GIS 6.4.3RC3/scripts/r.basin.py: line 80: import: command not found C:/Program Files (x86)/GRASS GIS 6.4.3RC3/scripts/r.basin.py: line 82: syntax error near unexpected token (' C:/Program Files (x86)/GRASS GIS 6.4.3RC3/scripts/r.basin.py: line 82:if not os.environ.has_key("GISBASE"):'

can someone help me? thanks


texts: Unzipped LiDAR LAZ files now I have +3,000 folders each with its own LAS file

texts: I am new to using LiDAR. I recently downloaded a massive amount of LiDAR. I unzipped using the LASTools, but now I have +3,000 folders, each containing a LAS file. Is there a trick to getting all of these files in to one folder so a LAS Dataset can be created?


texts: Layer map plots on top of one another in R/ggplot2

texts:
I'm making a choropleth map in R/ggplot2. I have the colored part working just fine. This code

p <- ggplot() +
  geom_map(data=c_1860,
           aes(map_id = GISJOIN, fill = AHL002),
           map = shp_1860) +
  expand_limits(x = shp_1860$long, y=shp_1860$lat) +
  scale_fill_gradient2(low="white", high="red") +
  theme_tufte()
print(p)

produces this map:

I can also make an outline of the United States without a problem. This code

us     <- map_data("usa")
p <- ggplot() +
  geom_path(data=us,
            aes(x=long, y=lat, group=group),
            color="black", size=0.5) +
  theme_tufte()
print(p)

produces this map:

I want to draw the outline around the choropleth to give the map a shape, but that doesn't work. This code just produces the same output as the choropleth map.

p <- ggplot() +
  geom_map(data=c_1860,
           aes(map_id = GISJOIN, fill = AHL002),
           map = shp_1860) +
  expand_limits(x = shp_1860$long, y=shp_1860$lat) +
  scale_fill_gradient2(low="white", high="red") +
  geom_path(data=us,
            aes(x=long, y=lat, group=group),
            color="black", size=0.5) +
  theme_tufte()
print(p)

How can I put the outline of the US around the choropleth map?


texts: Is it possible to add OS Open Carto vector tiles to QGIS 3?

texts:
I am interested to learn whether it is possible to view Ordnance Surveys OS Open Carto vector tiles in QGIS - I am using QGIS 3.4.2.

The data I would like to add can be found here, and I have tried to use the Vector Tiles Reader plugin to connect using the URL provided on the ArcGIS tiles server.

When I try this I get an error saying "The field 'vector_layers' is required but not found. This is invalid TileJSON"


texts: Can enhanced vegetation index have values greater than 1? I am using the ENVI software to calculate enhanced vegetation index (EVI) of a hyperspectral image and my EVI value range was from 0 to 1.9. Is it wrong?


texts: Error in overlay - formula not vectorized

texts:
I've run into another error while trying to get my predictive model to work. My base model is basically done (extracting locational parameters from archaeological sites in an area, and calculating a predictive map for archaeological risk areas, containing slope, tpi, tpi_15, aspect and distance to water), but it is not as accurate as I'd like it to be. So I'm trying to assign different weights to the parameters (e.g. distance to water is more important than aspect), and that's where I'm stuck.

weights <- apply(X = sn_sf_bind, MARGIN = 2, FUN = function(x)
{sqrt(1/sd(x)*range(x)[2]-range(x)[1])})

sn_sf_bind is my table with the archaeological sites and assigned locational parameters, looks like this:

For this I don't get an error, but I don't know if it could be in there already.

The code that isn't working is this:

wba <- overlay(x = ter_par_big,
               fun=function(w,x,y,z){
                 return((v*weights[1])+
                          (w*weights[2])+
                          (x*weights[3])+
                          (y*weights[4])+
                          (z*weights[5]))},
               unstack = TRUE)

I've also tried this:

fw <- function(v,w,x,y,z){
  return((v*weights[1])+
           (w*weights[2])+
           (x*weights[3])+
           (y*weights[4])+
           (z*weights[5]))}

wba <- overlay(x = ter_par_big,
               fun=Vectorize(fw),
               unstack = TRUE)

But it didn't change anything, I'm still getting the same error:

Error in .overlayList(x, fun = fun, filename = filename, forcefun = forcefun,  :
  cannot use this formula, probably because it is not vectorized

ter_par_big is my Raster Brick, consisting of 5 layers with the locational parameters for the overall area, it looks like this:


texts: QgsMarkerLineSymbolLayerV2 change default marker

texts:
I am trying to use the QgsMarkerLineSymbolLayerV2 line symbol layer in my app. I can only make it work using default circle as marker as below:

layer1 = QgsMarkerLineSymbolLayerV2(True, 3)

how do you change the marker from the default circle to a say, star (*)? Also, can I use a svg marker at all?


texts: How to get the upstream average of a raster data layer that falls into a river network?

texts:
I have a raster layer of a network of rivers and a raster layer of data I want to extract.

For each pixel of the rivers, I need to know the upstream average of the data layer that falls into the river. Kind of like a mix between Flow Accumulation and Zonal Statistics.

I have access to ArcGIS 10.1 and QGIS 2.8.2.


texts:
In Postgis/Postgres given a set of linestrings with topological connected edges, how can one calculate a unique id (network id) grouped by a set of attributes.
Installed postgis/pgrouting and looking at possibilities but so far could not find a solution.


texts: How to retrieve PROJECTION value for a map?

texts:
I am trying to read the Projection from GeoTIFF/shapefile, using GDAL 1.9.0 by the following command in the command line...

gdalinfo.exe -proj4 map-path

but its retuning empty string for "PROJECTION" parameter...  when I check the same map in QGIS getting the correct projection for this map.

Actually I want to retrieve "PROJECTION" value for any map from the C/C++ code, could you please let me know the way to retrieve this?


texts: Parcel Map - How Many Acres of Feature X fall Within Each Parcel

texts:
Can someone recommend a good tool for ArcGIS that would help summarize how many acres of a polygon feature occur within each parcel in a separate layer.

I have a layer containing a few hundred parcels stored as polygon features. Each parcel feature has an APN, calculated acreage, owner name/contact info, and a few other attributes. In a second layer, I have approximately 5000 smaller polygons with a five column attribute table. I need to figure out how many acres of each of those attribute values occur within each parcel. The issue is that the smaller parcels can potentially overlap multiple parcels.

Is there a tool that would be able to easily summarize this data for me? I'm envisioning something that slices those 5000 into separate polygons where they overlap a parcel boundary.

I've been learning on ArcGIS 10.2, but have QGIS installed on the computer if this would be easier in that program.


texts: Problem with import qgis.core when writing a stand-alone PyQGIS script

texts:
I am trying to write a stand alone python script in windows 7 that uses qgis.core and various other qgis libraries.

I have found an old question on this subject here: Writing standalone Python scripts using PyQGIS? and was following its advice, by using the guide here: http://docs.qgis.org/testing/en/docs/pyqgis_developer_cookbook/intro.html#python-applications

The code suggested there  is:

set PYTHONPATH=c:\qgispath\python

followed by:

set PATH=C:\qgispath;%PATH%

I found that IDLE did not like the "set at the beggining, so my code looks like this:

PYTHONPATH="C:\OSGeo4W64\apps\qgis\python"
PATH="C:\OSGeo4W64\apps\qgis\python\qgis\core;%PATH%"
import qgis.core

and I get:

Traceback (most recent call last):
  File "<pyshell#15>", line 1, in <module>
    import qgis.core
ImportError: No module named qgis.core

I have tried every variant of the paths I can think of including removing the ;%PATH%" at the end (that looks a bit superfluous) but have not been able to find any more information. Does anyone have any ideas where I might be going wrong?

EDIT:

After seeing Luigi Pirelli's comment I have now tried:

>>> PYTHONPATH="C:/OSGEO4~1/apps/qgis/./python\qgis"
>>> import qgis.core

Traceback (most recent call last):
  File "<pyshell#1>", line 1, in <module>
    import qgis.core
ImportError: No module named qgis.core
>>> PATH="C:\OSGeo4W64\apps\qgis\python\qgis\core;%PATH%"
>>> import qgis.core

Traceback (most recent call last):
  File "<pyshell#3>", line 1, in <module>
    import qgis.core
ImportError: No module named qgis.core
>>> PATH="C:/OSGEO4~1/apps/qgis/./python\qgis\core;%PATH%"
>>> import qgis.core

Traceback (most recent call last):
  File "<pyshell#6>", line 1, in <module>
    import qgis.core
ImportError: No module named qgis.core
>>>


texts: Getting coordinates from mouse click in QGIS 3 (python plugin)

texts:
I am trying to get coordinates from a mouse click in QGIS 3 and python 3 using the plugin builder.

I have tried using the response from this post, and i can manage to get the coordinates printed to the python console in QGIS but at the same time i get the message that "'TypeError: 'QgsMapToolEmitPoint' object is not iterable" or "'QgsMapToolEmitPoint' object does not support indexing"

I can't remember i ever seen a python program causing an error but still manage to compute and output the wanted result. The "not iterable" error traces back to the line with the for loop "for point in ... "
the "does not support indexing" error is caused if i do e.g. x = pointTool[0]

Here is a version of the code:

from qgis.gui import QgsMapToolEmitPoint

def display_point(pointTool):
    #print(pointTool)
    cood = []
    for point in pointTool:
        cood.append(point)
    print(cood)

# a reference to our map canvas
canvas = iface.mapCanvas()

# this QGIS tool emits as QgsPoint after each click on the map canvas
pointTool = QgsMapToolEmitPoint(canvas)

pointTool.canvasClicked.connect(display_point)

canvas.setMapTool(pointTool)

display_point(pointTool)


texts: Duplicating data for use in new ArcMap project

texts:
I am using ArcMap 10.8.1.
I wanted to make a new project based on the same data as an existing project, so I duplicated the project file and started editing one of the layers in Editor Mode. I then went back to my original project and realised that my edits had also affected it, because as I've now realised, I didn't duplicate the data, just the container for it.
I've done a fair bit of work aside from those silly edits on the new project already so I'd rather not start again.
Is there a straightforward way of duplicating the data and pointing the new project to that new set rather than editing the originals or do I have to go into each layer and edit the data source individually?
I had a search of the site but couldn't find anything specific for my question.


texts: How to delete openmaptiles (OSM Mbtiles) layers

texts:
I am struggling with converting PBF to Postgresql Table (importing only certain layers) or editing them, for the volume of the data. (world scale of PBF is about 80Gbs)

meanwhile, I have found that openmaptiles offers a very light mbtiles format of OSM data, though it is only up to zoom level 14

https://openmaptiles.com/downloads/tileset/osm/

From there mbtiles, I wanted to remove certain layers such as ocean polygons and other water boundaries (and maybe a few more as I wish).

I know mbtiles is a type of sqlite but I have realized that sqlite (mbtiles) viewer is not like QGIS layer like.

Will there be any GUI tools that I can handle mbtiles?


texts: Should I use a projection or not, and if so what SRID would be useful?

texts:
First off, I realize that this is perhaps not the most stackexchange-esque question. I am however having trouble finding a good source of information for this. So if there is a better place for me to look for this information please redirect me there.

I am creating an application for which the persistence layer is PostGIS, and I am unsure what if any projection I should use for storing my data.

My spatial data system is working on the following general requirements:

The data that I am storing pertain to events, and I would like to be
able to find events that are nearby another event. This is to say,
given event 'foo' I would like to find events within a certain radius
of 'foo.'
The data points are distributed at a national, and later at a
potentially global level, but the distance comparisons would only
occur at a relatively local level < 100k, so I only really care about
fidelity of distance up until that point.
Accuracy is much less important than speed in my use case, and single
digit percentage error is acceptable as long as the distance lookup
is quick.

I have looked up the geography type which is not tied to a particular position, however the PostGIS documentation appears to take a rather condescending view of the geography type:

  If you don't understand projections, and you don't want to learn about them, and you're prepared to accept the limitations in functionality available in GEOGRAPHY, then it might be easier for you to use GEOGRAPHY than GEOMETRY. Simply load your data up as longitude/latitude and go from there.

I am relatively comfortable with the mathematical notion of projection and the topology of manifolds, and I do care about performance, but I have little GIS specific knowledge and I have yet to find a place where there is a mapping between SRIDs and usecases. Is there a specific type of projection that I should be using, or a type that I should be looking for?


texts:
I built up a polygon layer from a road network for finding roundabouts using the compactness factor of a polygon. Then I have to intersect the extracted roundabout polygons with the road network so that I can find the composing links. This process can take ages to be done using for-loops for each roundabout and intersecting it with all links.

I am wondering if there is a way that I can add the composing links' ID to polygon in the very first step with out intersecting them.

Current algorithm takes three following steps:

FeatureToPolygon function in ArcPy to create the polygons.
Iterating over each polygon and computing the compactness,
And if it's over the threshold then select links which "SHARE_A_LINE_SEGMENT_WITH" the possible roundabout (the polygon in hand with compactness higher than threshold).

The dataset has more than 250000 line features but the number of detected roundabouts are depending on the threshold. it should not be more than 1000.

There is an optional argument about attributes. The ArcGIS (or ArcPy) instruction says not to use it though.

I tried the following code but I get an error.

polyoutputaddress = outputaddress+output_polygon_name
arcpy.FeatureToPolygon_management(inputaddress,polyoutputaddress,"", "NO_ATTRIBUTES", "")

targetFeatures = polyoutputaddress
joinFeatures = inputaddress

outfc = outputaddress+"fieldmapping"+output_polygon_name

fieldmappings = arcpy.FieldMappings()
fieldmappings.addTable(targetFeatures)
fieldmappings.addTable(joinFeatures)

inputFIDIndex = fieldmappings.findFieldMapIndex("JOIN_FID")
inputfieldmap = fieldmappings.getFieldMap(inputFIDIndex)
inputfieldmap
print inputFIDIndex

inputfieldmap.mergeRule = 'Join'
inputfieldmap.joinDelimiter = ','
fieldmappings.replaceFieldMap(inputFIDIndex, inputfieldmap)

arcpy.SpatialJoin_analysis(targetFeatures, joinFeatures, outfc, "#", "#", fieldmappings,"SHARE_A_LINE_SEGMENT_WITH")

here is the error:
    arcgisscripting.ExecuteError: ERROR 001156: Failed on input OID 0, could not write value '-1,-1,-1' to output field JOIN_FID Failed to execute (SpatialJoin).

as I checked this is because the FID field is numeric, but Join rule should be used for text values.

Any work around this problem?


texts: Why QGIS crashes when constructing QgsGraduatedSymbolRenderer object by QgsGraduatedSymbolRenderer.createRenderer() method?

texts:
Today I want to render a vector layer in QGIS by using PyQGIS QgsGraduatedSymbolRenderer renderer. When I directly construct it by using the constructor, it works. But when I want to use the createRenderer() method, QGIS then crashed unexpectedly. I have executed my code in QGIS 3.14.0 and QGIS 3.16.3 python console and QGIS crashed all. I want to know if there is a problem in my code or a bug with QGIS.
Here is the code:
LayerPath = r"C:\Users\15239\Desktop\QGIS\dln.gpkg"
DelaunayLayer = QgsVectorLayer(LayerPath,"Delaunay","ogr")
QgsProject.instance().addMapLayer(DelaunayLayer)
CpyRamp = QgsCptCityColorRamp("cb/div/BrBG_","05")
SimpleSymbol = QgsFillSymbol([QgsSimpleFillSymbolLayer(QColor(0,50,70))])
DelaunayLayer.setRenderer(QgsGraduatedSymbolRenderer.createRenderer(DelaunayLayer,'AREA',5,1,SimpleSymbol,CpyRamp))
DelaunayLayer.triggerRepaint()

Here are the test data and its attribute table:


texts: How to create a custom transverse mercator projection?

texts:
I'm working on a script that performs some simple geometric manipulations such as buffering vector data using shapely or geopandas. The input (i.e. original) shapefile could be geographically located anywhere on Earth and be in any spatial reference system. I need to ensure the coordinate system is standardized and projected prior to running the geometric operations because preserving distance is important to my analysis.

How can I create a custom transverse mercator projected coordinate system using geopandas and/or shapely?

The following is what I'm trying to accomplish:

import geopandas as gpd

in_shp = r'C:\path\to\input_shapefile.shp'
out_shp = r'C:\path\to\output_shapefile.shp'

# GeoPandas
geometry = gpd.read_file(in_shp)

# Project to a custom transverse mercator projected coordinate system
sr = {'a custom transverse mercator proj4 string'} # Here's the main issue
projected_geometry = geometry.to_crs(sr)

# Buffer 100m
buffered_geometry = projected_geometry.buffer(100)

# Convert to WGS84
final_product = buffered_geometry.to_crs({'init': 'epsg:4326'})

# Write to disk
final_product.to_file(out_shp)


texts: Split polygon with holes to get several hole-less polygons

texts:
Unlike this thread I actually need to split the polygons into several hole-less polygons. I'm not sure on how to approach to this problem.

I have polygons with holes of many diferent shapes, like buildings, esplanades, rivers, etc. And I need to split them in a coherent way (if the polygon is a river it should be splitted perpendiculary to his axis) and with less splits as possible.

I need to do it programatically. I'm using ArcObjects 9.3, but if there is no other way I could try to use a specific algorithm to solve this problem.


texts: Updating features from my active layers in QGIS

texts:
I'm trying to write a Python script with QGIS 3.4 but I have troubles when I want to change attributes values from a field.

With the code below, new values are calculated and the python console returns me True, but when I'm charging the layers in my QGIS Project, the attributes of all my features have a value = 0.000.

My Python code looks like:

champs = couche.fields()
nom_champs = [champ.name() for champ in champs]
feature=couche.getFeatures()
for elem in feature:
      elem.setAttribute('STD_2D', np.sqrt(elem[5]*elem[5]+elem[6]*elem[6]))
      elem.setAttribute('STD_3D', np.sqrt(elem[5]*elem[5]+elem[6]*elem[6]+ elem[7]*elem[7]))

      couche.changeAttributeValue(elem.id(), 11, elem.attributes()[11], skipDefaultValues=True)
      couche.changeAttributeValue(elem.id(), 12, elem.attributes()[12], skipDefaultValues=True)

      couche.dataProvider().updateExtents()

I think that they are not updated despite my command: couche.dataProvider().updateExtents(). Do I need the other command like this one? Is it a problem due to how I'm writing my script?

I tried to find help with subjects like this one but without success.

Edit :
I tried to use 'couche.startEditing()' or 'couche.commitChanges()' but when I'm doing that, the attributes of the field I have to change take no values. Attributes cases are blank.

If It can help, here are previous lines of my code:

if couche.isEditable() == False :
   couche.startEditing()

res=couche.dataProvider().addAttributes([QgsField("STD_2D", QVariant.Double,len=5,prec=4), QgsField("STD_3D",QVariant.Double,len=5,prec=4)])
couche.updateFields()

Edit 2:
For those who have the same probleme, I solved it by using the library "csv" with python, without PyGis. Then I import the new csv with PyGis in my current project.

import csv

entree = csv.reader(open("Export_QGIS_bon.csv","r"),delimiter=',')
sortie = csv.writer(open("Pour_shp.csv", "w"),delimiter=',')
row = next(entree)
row.insert(13,'STD_2D')
row.insert(14,'STD_3D')
sortie.writerow(row)
for row in entree:
   row.insert(13,round(float(np.sqrt(float(row[5])*float(row[5])+
   float(row[6])*float(row[6]))),3))
   row.insert(14,round(float(np.sqrt(float(row[5])*float(row[5])+
   float(row[6])*float(row[6])+float(row[7])*float(row[7]))),3))
   sortie.writerow(row)


texts: Loop to check multiple polygons overlap in r

texts:
I have a set of polygons (available here: http://ulozto.cz/x39ELbFt/polygons-zip) created in QGis and processed in R. I would like to know if they are overlapping or not, and if yes, where is the overlapped area (and between which polygons)? I know that to figure this out I can use tools:

over{sp} - return data.frame not a spatial object!
intersect{raster}
gIntersection{rgeos}

I can find intersection by comparing 1 by 1 polygon each time step. However as my real data are huge, I prefer to make a loop and include the iteration of polygons:

intersect(p1, p2)
intersect(p1, p3)
intersect(p1, p4)
intersect(p2, p3)
etc...

For the moment, I made up only loop for one polygon intersecting another 3. From resulting list of overlapping polygons I can't however really say which polygons are overlapped. And if they are not (p1 and p4) this NULL value is not included in my resulting list!

library(sp)
library(raster)
library(rgeos)
library(spatstat)
library(rgdal)

setwd("D:/...")
p1<- readOGR(dsn=getwd(), layer="p1")
p2<- readOGR(dsn=getwd(), layer="p2")
p3<- readOGR(dsn=getwd(), layer="p3")
p4<- readOGR(dsn=getwd(), layer="p4")

ab<-list(p2, p3, p4)

# my loop only for 1 polygon, not for combinations of polygons
by two intersection methods

# intersection by gIntersection{rgeos}
int1<-list()
for (i in 1:length(ab)){
int11<-gIntersection(p1, ab[[i]], byid = FALSE)
int1[[i]]<-int11
}

# intersection by intersect{raster}
rst.int<-list()
for (i in 1:length(ab)){
rst.int1<-intersect(p1, ab[[i]])
rst.int[[i]]<-rst.int1
}


texts: Unable to locate "Package for QField" in the drop down menu. (See pic below)


texts: Showing postive and negative range values separately in legend of ArcMap

texts: The legend here shows it's positive and negative ranges separately with a decreasing and increasing order with circles. I know it is related to the graduated symbology and multiple attributes but which part could I do for it I am getting no idea. I can't show ranges together.


texts: I've already used Sentinel-2 images with the resolution of 10m/pxl for urban heat island analysis. Now I need higher resolution image (1 to 5m per pixel) for Ekaterinburg in Russia. Searching in the web I found some commercial sources like KOMPSAT 2 or sentinel-hub, but can't understand wich source has LST band that I need for analysis? Need help.


texts: When you create a map in Google Fusion Tables, there's an option to weight your geographical data points against data in another column. What is the math behind this weighting?


texts: How do I create a dissolved buffer using PostGIS with EPSG 4326?

texts:
I'm trying to buffer a polyline (WGS84/EPSG 4326) by 50km, and then dissolve all the overlaps in the buffer. So far I've tried a combination of:

casting the EPSG 4326 geometry as geography (to get accurate buffer) then
creating a dissolved buffer using ST_union(ST_buffer(geo.my_polyline)).

However, I'm getting a syntax error around my CAST close-parenthesis: ERROR: syntax error at or near ")"

Here is my code:

create table geo.fifty_km_buffer as (
select id, ST_union(ST_buffer(
    CAST(geo.polyline.geom) as geography,
    50000
    )) as geom
from geo.polyline
)


texts: How can I use Oracle views with Geoserver and WFS-T?

texts:
What goes on when Geoserver handles updates to my table using WFS-T?

I have a table

AREA (
ID NUMBER(10) NOT NULL PRIMARY KEY,
SOME_OTHER_ID NUMBER(10) NOT NULL FOREIGN KEY REFERENCES FOO(ID),
GEOMETRY MDSYS.SDO_GEOMETRY NOT NULL
);

and a sequence SEQ_AREA for the primary key values.

When I publish this table as a Geoserver layer, it seems that it automatically detects that ID is the primary key and does not show it in the layer's properties.

When a new row is inserted, Geoserver does not use the index - how does it determine the next value? Does it do SELECT MAX(id) FROM AREA or something like that? Could I make it use the sequence?

And - if I create a view and a layer from that - the view does not have primary key as such - how can I make Geoserver understand it so that it is updatable? When I create a view such as:

CREATE VIEW TESTVIEW AS
  SELECT a.ID, foo.somefield, a.some_other_id, a.GEOMETRY
  FROM FOO foo, AREA a
  WHERE foo.ID = a.SOME_OTHER_ID;

insert into USER_SDO_GEOM_METADATA (table_name, column_name, diminfo, srid)
values ('TESTVIEW ',
'GEOMETRY',
SDO_DIM_ARRAY
(
  SDO_DIM_ELEMENT('X',50199, 761274,0.5),
  SDO_DIM_ELEMENT('Y',6582464, 7799839,0.5)
),
3047);

The layer ends up being viewable OK, but cannot be updated / inserted. When updates are attempted through wfs-t I get:

2013-02-26 19:00:33,636 ERROR [org.geoserver.ows] -
org.geoserver.wfs.WFSTransactionException: {http://www.opengeospatial.net/cite}TESTVIEW is read-only


texts: Count polygons that intersects polygon in Spatialite

texts:
I am trying to do a very simple query, but keep getting an error. I am looking to get the count of polygons that intersect with another polygon. I've used the following query:

SELECT feature__1 AS type, Count(feature__1) AS type_cnt
FROM cityprj_munparcel
WHERE (SELECT ST_Intersects(cityprj_munparcel.geom, to_districts.geom))
GROUP BY type

But keep getting the following error: no such column: to_districts.geom

I have looked at the table and have performed other queries with the layer and the the column is there. Unsure what is wrong.


texts: How do I install an external Python package in QGIS?

texts:
Is it possible to install an external package in QGIS?
For example: Can I install ReportLab in order to create a plugin to make pdf reports?

Or Pandas. How can I install Pandas in QGIS Python library?


texts: Transforming a single-band GeoTIFF into an RGBA PNG file using GeoTools

texts:
I'm fairly new with GeoTools.
I'm using it in order to generate GeoTIFF files that have a single band, with values between -1 and 1, or NaN.

An example of such a file is available here: https://m-innvoincy.box.com/s/ow1fentcgddgdeisfwj4i44yraeqkh4l

If you load this into QGIS or something similar, it should look like this (by default):

The next step, which I'm struggling with, is to generate a multi-color PNG file, that will be a visual representation of the raw data in the GeoTIFF (I'm ignoring the geography for now).
I would obviously need to specify a mapping between ranges of values in the GeoTIFF and RGBA values. NaN values should be mapped into the alpha channel ('transparent' pixels).

I would like to do it with GeoTools, but haven't found a way to do it yet.
It's important to note that the code should run on the server-side, there is no Map object and no GUI, Swing or otherwise, should be involved.


texts:
OK, I know that this is Python 101, but I can't seem to figure it out.

I have a function setup that will return a string, specifically an azimuth in DMS. I need to take the return values and pull out only the value that appears prior to an ":"

So I may get a value
x = 123:34:12.4

I need x to then equal 123 and use that to updateRow in an attribute table.

As usual... THANKS!


texts: gdal_merge.py just copies first input image

texts:
I am working with Sentinel-2 data, which I want to collect for a ROI that unfortunately is covered by two different (neighboring) granules. After downloading the granules, I cut them and convert them to to GeoTiff using

gdalwarp -of GTIff -crop_to_cutline -cutline cutline <10m_input_bands> datasetname_10m.tif

This gives the corresponding 4-band GeoTiffs (which contain only the 10m bands 4,3,2 and 8), including lots of no-data pixels where the ROI is actually located in the neighboring granule.

After that, I want to merge them by

python gdal_merge.py -of GTiff -o mergedImage.tif inputImage1.tif InputImage2.tif

The result, however, seems to be just a copy of inputImage1.tif, i.e. the no-data values remain and no mosaic is created.

What am I missing? Or is there another way to create coherent Sentinel-2 ROI images extracted from neighboring granules? At best, this would be automated for batch processing, i.e. search for products of similar name where only the granule indicator (e.h. UPU, UQU etc.) changes, and then performs the merging automatically.

EDIT: Maybe this was clear already, but to be sure: Since the output of gdalwarp is already corresponding to the ROI, the images are of the same size and follow the same geocoding reference (i.e. they reflect the same spatial subset).


texts:
I have a graph extracted from an American city: vertices represent road intersections, edges/arcs represent road segments connecting the intersections. For each road segment, I have the length and allowed driving speed. Furthermore, roads belong to different classes such as highway/residential/primary roads/....

I am looking for a routing algorithm which does the following:

Compute a pairwise travel time matrix for each intersection (i.e how long it would take to travel from intersection i to intersection j.
Query the quickest route from intersection i to j and returns a path (sequence of edges/arcs).
The algorithm should be able to take updates into consideration. For instance, if there is a congestion on road segment (x,y) which increases its anticipated travel time, or if road segment (x,y) becomes blocked, the data structures should be updated incrementally, thereby updating all affected shortest paths (i.e. update all shortest paths through (x,y)).

The distance matrix will be queried very often so this should be precomputed, but the actual shortest paths can be computed at a later stage if that's more efficient. The number of nodes in the graph is roughly 5000-10000.

I was hoping that people could refer me to literature about this, preferably something that doesn't take weeks to implement. Obviously, simply running Dijkstra for each intersection takes too long. Running Floyd–Warshall all-pairs shortest path algorithm could be a possibility, but I'm not so sure how to update the result if the cost of one of the edges/arcs changes. Furthermore, it seems to be a slow approach since the algorithm ignores that most of the routes go through a small subset of the road segments.

I'm not necessarily looking for an exact approach. A good approximation suffices. Any pointers to literature or algorithm descriptions are welcome. I presume that any kind of routing software needs a similar algorithm. I'm mainly interested in something that works; it shouldn't necessarily be state of the art.

If this isn't the right place to ask, please let me know (especially if you know a better place :) ).
Just to emphasize: I'm not looking for a software/database package, I'm looking for an algorithm.


texts: How to update raster color ramp with new min/max values using PyQGIS?

texts:
I'm trying to update a single band pseudo color renderer with new min/max values from the python console.  All the code samples i have been able to find deal with constructing a renderer from scratch, but in this case I have a color ramp setup for an existing layer and I want to be able to stretch to color ramp with new min/max values.

I can set the new min/max values for the renderer, I just haven't figured out how to reclassify the color ramp (I'm looking for the equivalent of hitting the "classify" button in the style properties dialog).

This is what I have so far:

l = iface.activeLayer()

provider = l.dataProvider()
extent = canvas.extent()
stats = provider.bandStatistics(1, QgsRasterBandStats.All, extent)

min = stats.minimumValue
max = stats.maximumValue

l.renderer().setClassificationMin(min)
l.renderer().setClassificationMax(max)


texts: Iterate the shortest path point to layer

texts: I have two layers. First layer defines the starting points and the second defines the end points. I just want to iterate the same algorithm "point to layer" by considering all the points of the first layer one by one and get the shortest point.


texts: How to get exact ISymbol image size in pixels

texts:
I found the next code to convert ISymbol into bitmap image,
but how to know the suitable image size? (suppose we have a ITextSymbol)

private static Image PointSymbolToImage(ISymbol symbol, Size imageSize)
{
    double x = imageSize.Width / 2;
    double y = (imageSize.Height - 1) / 2;

    IPoint point = new PointClass();
    point.PutCoords(x, y);

    Bitmap bmp = new Bitmap(imageSize.Width, imageSize.Height);
    Graphics g = Graphics.FromImage(bmp);

    symbol.SetupDC(g.GetHdc().ToInt32(), null);
    try { symbol.Draw(point); }
    catch { }
    symbol.ResetDC();

    g.Dispose();

    return (Image)bmp;
}

For more details:

I used QueryBoundary but the dimensions returned after transforming to screen coordinates were wrong.
Next is my code:

Here symbol is ITextSymbol and (x, y) are the lat, lon coordinates of the symbol.

IPolygon boundary = new PolygonClass();
boundary.SetEmpty();

IDisplay screenDisplay = axMapControl1.ActiveView.ScreenDisplay as IDisplay;
IDisplayTransformation displayTransformation = screenDisplay.DisplayTransformation;
IPoint point = new PointClass();
point.PutCoords(x, y);

symbol.QueryBoundary(screenDisplay.hDC, displayTransformation, point, boundary);

IPoint mapPoint = new PointClass();

int minX, minY;
mapPoint.PutCoords(boundary.Envelope.XMin, boundary.Envelope.YMin);
displayTransformation.FromMapPoint(mapPoint, out minX, out minY);

int maxX, maxY;
mapPoint.PutCoords(boundary.Envelope.XMax, boundary.Envelope.YMax);
displayTransformation.FromMapPoint(mapPoint, out maxX, out maxY);


texts:
I am using ArcView without any extensions.  Is there an alternate program out there which uses Cadastral Editing (Parcel Editing) or something similar?

I downloaded the ArcEditor 10.0 trial version but the parcel editing does not function as I don't know if its due to it being trial.


texts:
I'm new at Geoserver.

I've data in postgis tables and using QGIS OpenGeo plugin to publish data to Geoserver.
Postgis tables --> QGIS (Opengeo) --> Geoserver

I would like to give a different name (alias) to the columns of my tables and I did that at QGIS Layer Properties, but when I save the style in a .sld file but such information is not stored in the file  (in a .qml style file all the information is stored).

So, as anyone know if it's possible to give aliases to columns at geoserver (or any kind of workaround...) ?


texts: Python Create Dictionary with Shapely Points from Linestring coordinates

texts:
I have a Geodataframe with Linestring Z geometries as follows:

id_line
geometry

1
LINESTRING Z (0 0 0, 5 6 9, 1 2 3, ...)

2
LINESTRING Z (7 8 9, 9 8 7, 1 2 5, ...)

each of the Linestring is made up of n triples of x, y, z coordinates.
What I am looking for? Extract each of the triples that make up each line and convert it into Shapely POINT, then save each of these points in a dictionary like that:
{id_line : [points]}
points_by_line = {1 : [POINT Z(0,0,0), POINT Z(5,6,9), POINT Z(1,2,3)], 2 : [POINT Z(7,8,9), POINT Z(9,8,7), POINT Z(1,2,5)]

How am I trying to do this?
#Extracting the geometries from the Geodataframe
lines_segments = real_traj['geometry']
# saving the coordinates in a list
line_coords = []
for line in lines_segments:
    coords_line = list(line.coords)
    line_coords.append(coords_line)

Then I am trying to create a function that will extract each of the coordinates and convert them into Shapely Points, but here is where I am stuck
def linestring_to_points(line_coords):
    # creating a list for stock the points of each line
    line_points = []
    for coords in line_coords:
        for point_coord in coords:
            # access accessing the coordinates of each point on the line and creating a point from
            # each group of coordinates
            line_points.append(geometry.Point(point_coord))
    return {id_line:Point.coords}


texts: How to monitor changes in a PostGIS database in GeoDjango?

texts:
I have been working on a GeoDjango project which uses a PostGIS database.I have enabled the history in my application using the django-simple-history module and works all fine.
The issues comes in when I edit a record using QGIS by using the PostGIS tool in QGIS. The changes I make using QGIS are not recorded in the history of the record. These include operations such as split features and merge (editing the geometry of the record).
First, is it the most ideal way of doing things in such a scenario?
Secondly, how do I ensure the changes in a record made through QGIS are tracked in the historical_model which can be visible in my GeoDjango application?


texts: How to work with the service table layer_style in QGIS

texts:
Help me please.

Current versions:
QGIS 2.14
Sql Server 2012

How to save a layer style in a service table layer_styles?


texts: Rasterio does not execute the rasterize function

texts:
Rasterio's features.rasterize() function is not working as expected.

I have two files:

OpenStreetMap shapefile with land usages
Raster (.tif) with satellite imagery

Goal: Create a raster file from the OSM shapefile where attribute type = 'parking'.
Requirement: The raster file should have the same shape and projection as the .tif satellite imagery file.

Here is what they look like loaded into QGIS. You can see the feature polygons in green over the satellite background

Here are the projection and transform

using rasterio.open().meta.

Raster file

{'affine': Affine(1.0, 0.0, 543673.0, 0.0, -1.0, 4192469.0),
 'count': 3,
 'crs': CRS({'init': u'epsg:26910'}),
 'driver': u'GTiff',
 'dtype': 'uint8',
 'height': 21509,
 'nodata': None,
 'transform': (543673.0, 1.0, 0.0, 4192469.0, 0.0, -1.0),
 'width': 17245}

and using fiona.open().meta

Vector file

{'crs': {'init': u'epsg:4326'},
'crs_wkt': u'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]',
'driver': u'ESRI Shapefile',
'schema': {'geometry': 'Polygon',
'properties': OrderedDict([(u'id', 'int:11'),
    (u'osm_id', 'float:19'),
    (u'name', 'str:104'),
    (u'type', 'str:17'),
    (u'area', 'float:32.10'),
    (u'z_order', 'int:11')])
}}

Here's the code:

import rasterio
import numpy as np
from rasterio import features
import fiona

#Load vector file with fiona
vector_file = 'shapefiles/sf_osm/san-francisco_california_osm_landusages.shp'

shapefile = fiona.open(vector_file)
geom = [shapes['geometry'] for shapes in shapefile]
attrib = [shapes['properties'] for shapes in shapefile]

#filter to features where type = parking
    type_match_idx = [i for i, feature in enumerate(attrib) if feature['type'] == 'parking']
    feature_geom = [geom[i] for i in type_match_idx]
    feature_attrib = [attrib[i] for i in type_match_idx]

#get metadata for satellite imagery
rst = rasterio.open('satellite.tif')

Everything is working fine up until here. This is where it stops acting as expected.

#rasterize using the shape and transform of the satellite image
image = features.rasterize(feature_geom, out_shape=rst.shape, transform=rst.transform)

print(np.max(image)) #results in 0

This produced a numpy array with every value equal to zero. I would expect some to be non-zero where features exist.

#saving image
with rasterio.open(
        'rasterized-results2.tif', 'w',
        driver='GTiff',
        transform = rst.transform,
        dtype=rasterio.uint8,
        count=1,
        width=rst.width,
        height=rst.height) as dst:
    dst.write(image, indexes=1)

This produces a blank (all black) image.

Why is features.rasterize() producing no data?


texts:
I want to find out the direction from the 2 points from postGIS database. Please see the attached image. I want to find out the direction from the mainPt to any point. But want to get the result as NW, NE, SE or SW. Simply if I get the degree of the angle it can easily be solved.

Thanks


texts: Exporting Labels to Annotations using ArcPy?

texts:
I'm trying to automatically generate annotations using arcpy from a set of MXD files with preorganised labels (all settings done, ready to export the labels). I'm using ArcGIS 10.3 with an Advanced Licence. All I'm looking to do is scripting the process of right-clicking the layer selecting an output and exporting the labels to annotations (see screenshot below):

So far, I found various similar questions that all point to the arcpy.TilesLabelsToAnnotations_cartography() tool.

Labeling features and converting them to annotations with ArcPy?
Automate converting labels to annotation in ArcMap at multiple scales?
Exporting and editing annotation feature class using ArcPy?
Labels to Annotations in ModelBuilder --> this tool looks most promising, but is not available for ArcGIS 10.3...

So I tried using the tool as follows, which for some reason does not end in an error (green checkbox appears after only a few seconds) but simply produces nothing at all...

Where am I going wrong?

Any suggestions for changes?

Are there easier ways of achieving my goal?

Ideally I don't even want a tile for exporting, since I want to convert all labels, no matter of the extent. Sidenote: Some of the resulting annotation feature classes are quite large (300+ MB) when exported manually, I don't know whether this is of any relevance for the tool.

import arcpy
import sys

MXD_folder = r"C:\Users\sutter\Desktop\GIS_temp\stack_test"
outGDB = r"C:\Users\sutter\Desktop\GIS_temp\stack_test\testOut.gdb"

arcpy.env.workspace = MXD_folder
arcpy.env.scratchworkspace = MXD_folder
arcpy.env.parallelProcessingFactor = "100%"
arcpy.env.overwriteOutput = True

#-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

# 2. Create surrounding polygon:

spatRef = arcpy.SpatialReference(2056) # Swiss LV95
arcpy.CreateFeatureclass_management(MXD_folder, "tempTile", "POLYGON", "", "", "", spatRef)
tempTile = r"{}\{}.shp".format(MXD_folder, "tempTile")
# Setup Extent:
array = arcpy.Array([arcpy.Point(2600000, 1200000),
                     arcpy.Point(2600000, 1210000),
                     arcpy.Point(2610000, 1210000),
                     arcpy.Point(2610000, 1200000)
                     ])
polygon = arcpy.Polygon(array)
with arcpy.da.InsertCursor(tempTile, ['SHAPE@']) as cursor:
    cursor.insertRow([polygon])
# Create Feature Layer for following tools:
arcpy.MakeFeatureLayer_management(tempTile, "tile")

#-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

# 3. Find all files for processing:
mxdNameList = arcpy.ListFiles("*.mxd")

#-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

# 4. Export to annotations for all MXD files:

for i in range(len(mxdNameList)):
    mxd = r"{}\{}".format(MXD_folder, mxdNameList[i])

    # List Layers in MXD. If less or more than 4, exit. If 4, activate labels
    mxdPy = arcpy.mapping.MapDocument(mxd)
    layers = arcpy.mapping.ListLayers(mxdPy)
    if len(layers) != 4:
        sys.exit()
    else:
        for layer in layers:
            layer.showLabels = True
        arcpy.RefreshActiveView()
        mxdPy.save()
    del mxdPy

    # Get Reference Scale from MXD Name:
    refScale = int(mxdNameList[i].split("_")[-1].split(".")[0])

    # Define Dataframe name:
    df = "Layers"

    # Start conversion of all layers to annotations:
    arcpy.TiledLabelsToAnnotation_cartography(mxd, df, "tile", outGDB, "group", "anno", str(refScale), "", "Id", "", "", "STANDARD", "GENERATE_UNPLACED_ANNOTATION")


texts:
I made a new plugin for QGIS and I want to include it to the Processing toolbar. I can't find a good hint on how to do it cleanly (applying answer of this question on my algorithms seems dirty). So, how to do it?

I built the plugin using Plugin Builder and I already published it (FYI, it's this one). However, I would like to enable processing-specific options like using batch.


texts: QGIS SAGA Overland flow distance to channel network not working

texts:
I am trying to run the QGIS SAGA tool "Overland flow distance to channel network" with a DEM and the channel network raster developed from this same DEM using the SAGA tool "Channel network." However, I am getting an error:
"Error: Channel Network"
I'm using QGIS 3.4.2-Madeira.
What could the problem be? Has this tool worked for anyone?


texts: Identify intersections and assign them to one of the features in a single layer

texts:
I have a single SHP-layer (12,810 features) and want first of all to identify intersections (incl. duplicates --> features with identical geometries). In the next step I want to delete/clip these intersections and assign them to one of the features so no intersections are present anymore. It would be best to assign the overlaps to those features where the area fraction of the overlaps is greatest (in the best case).

I´ve tried the Intersection Tool but in this case all of the features will be automatically selected (cause working only with 1 layer).
Using the Topology Checker Plugin I get the table with overlap errors (using the rules: must not overlap + must not have duplicates). As a result all of the features with overlapping parts are highlighted but no layer has been created so I cannot create a spatial selection. But this also doesn´t help since I only want to identify the overlapping parts of polygons within a single layer and not select the whole polygons with overlapping parts (To select those features I can use the solution proposed by @Hast.
Identifying overlapping polygons in single layer using QGIS?

Using the Check Geometries Plugin (similar settings to the Topology Checker Plugin) I only got again a a new layer (idential to the original one).

I suppose for this task either a PyQGIS script or a model in QGIS modeler is necessary. But I have no idea which QGIS processing tools should I use or which workaround.

Do you have any ideas?

original layer with overlaps:

layer with intersection features (no overlaps):

final dissolved layer:

Update: How described below after running the Intersection tool I split the field with the FIDs with overlaps (e.g. 1|5|18 --> 1) via the field calculator and ran the Dissolvetool based on this field. Before that I ran the Fix geometrytool to repair the geometry (lots of editing errors were present). As you can see on the last screenshot some features haven´t been dissolved although for each feature at least 1 FID exists. (E.g. FID 27 in the original layer and 2 features which overlap with this one, so in the intersection layer 3 features with FID = 27 are present). After running the Dissolvetool, the Dissolved layer only contains the 2 intersection features but not that from the original layer (green))...

Do you have an idea why some features couldn´t been dissolved?

Using the Select by location tool (using the overlap parameter) (original layer vs. dissolved layer) selected much more features then it should which I also don´t understand. I manually selected the missing features (green) and merged them to the dissolved layer via Merge vector layers. After that I again ran the Dissolvedtool but again some intersections were still present.


texts:
I have multiple features that I display, and on click a popup appears, with a button to store it to favorites (with intention of show only favorites afterwards)

The thing is, when I click the button, I retrieve the selected feature through the select interaction, but the feature that retrieves is not a feature type, is an object like this:

featureAtClick Object { $a: Object, ta: Object, qa: Object, f: 1, ko: 568, I: Object, a: undefined, : "geometry", g: null, l: undefined, 2 more… }

Instead of:

Object { type: "Feature", geometry: Object, properties: Object }

The properties of my feature are all inside of the I key from the first object, but for a reason I can not access to it.

I tried this:

Object.keys(selected).forEach(function (prop){console.log(selected[prop])})
console.log('selected',selected['I'])

But it returns only empty objects.

How do I convert it to a feature type, so I can add it to a featurecollection and display it later?

EDIT

As requested I add the click event code. This is the click on the feature:

favorites = { "type": "FeatureCollection", "features": [] }
                map.on('click', function (evt) {
                    var coordinate = evt.coordinate;
                    map.forEachFeatureAtPixel(evt.pixel, function (feature, layer) {
                        if (feature) {
                            console.log('favorites', favorites)
                            console.log('featureatclick',feature)
                            if (favorites.features.indexOf(feature) === -1) {
                                $('#fav-select').text('☆')
                            } else {
                                $('#fav-select').text('★')
                            }
    /// So the idea is that if the feature is on the favorite list, display either one text or another
                            $('#popup-content').html('
Name: ' + feature.get('name') + '
');
                            $('#popup-content').append('
Address: ' + feature.get('address') + '
');
                            if (feature.get('rank') != undefined) {
                                $('#popup-content').append('
Rank: ' + feature.get('rank') + '
');
                                $('#popup-content').append('
Reviews: ' + feature.get('reviews') + '
');
                            }
                            overlay.setPosition(coordinate);
                            map.addOverlay(overlay)

                        }
                    });
                });

And this the click on the favorites button:

                $('#fav-select').on('click', function (t) {

                var selected = selectInteraction.getFeatures()
                Object.keys(selected).forEach(function (prop){console.log(selected[prop])})
                console.log('selected',selected['I'])
/// Here I retrieve the feature that is selected again, and check if it is on the list to either remove it or add it
                if (favorites.features.indexOf(selected) === -1) {
                    console.log('fav')
                    favorites.features.push(selected)
                     $('#fav-select').text('★')
                     console.log('favoritesfavselect', favorites)
                } else {
                    console.log('unfav')
                    $('#fav-select').text('☆')
                    favorites.features.splice(favorites.features.indexOf(selected), 1)
                    console.log('favoritesfavselect', favorites)
                }

                console.log('selected', selected)
                console.log()

            })


texts: Normality, Ordinary Least Squares (OLS) and Geographically Weighted Regression (GWR)

texts: Most of the publications I have reviewed, conduct both OLS and GWR for comparison, and prior to analysis, ensure dependent and explanatory variables are normally distributed. Can someone explain why this is necessary for GWR? Are both OLS and GWR parametric? Is there a non-parametric version of GWR?


texts: Spatial Join Using Python /ArcGIS Addin assistant

texts:
I have a point shapefile that contains all my field collected data, lat/long and all the field data collected during the season.  The one thing not imported was the name of the Lakes these points fell on.  I have a polygon data set that has the name of the lake.  What I want to accomplish without creating another data set is to have the Lake name populated for each point that falls within the lake.  Some of my lakes do not have names so some will have blank records if that makes a difference.  I was trying to accomplish this with python code that I found around the net trying to accomplish the same thing.  I tried importing it into the ArcGIS addin assistant that lets you create buttons with python code running them.  Here is the code I have so far that is not working.  I click the button and nothing happens.

  import arcpy
  import pythonaddins

  class Transfer(object):
"""Implementation for AIS_Transfer_addin.tool (Button)"""
def __init__(self):
    self.enabled = True
    self.checked = False
def onClick(self):
    pythonaddins.MessageBox ('Attributes have been Transfered!', 'Message', 0)
    arcpy.env.overwriteOutput = True
    arcpy.MakeFeatureLayer_management('E:\Violations\AIS_Registry\MISIN_Source_Data_All_20181120.shp','lyr_Source')
    arcpy.MakeFeatureLayer_management('E:\Violations\AIS_Registry\Michigan_Lakes_All_v3.shp','lyr_lakes')

    rows = arcpy.SearchCursor(lyr_lakes)
    for row in rows:
        arcpy.SelectLayerByAttribute_management('lyr_lakes', 'NEW_SELECTION', '\'FID\' = ' + str(row.getValue('FID')))
        arcpy.SelectLayerByLocation_management('lyr_source', INTERSECT, 'lyr_lakes', '', 'NEW_SELECTION')
        arcpy.CalculateField_management('lyr_source', 'Lake_Name', "'[0]'".format(str(row.getValue('Label'))), 'PYTHON_9.3', '')
    print 'Finished processing' + str(row.getValue('Label'))


texts: Turning part UTM, part degree coordinate into location on Google Maps?

texts:
I'm really new at this and not actually sure what kind of coordinate this is but I need to know where this coordinate leads to.
The cöordinates are part of a province wide "fox hunt" in gelderland, The Netherlands. I basically have to solve an obscure hint to find the exact location via coördinates. This hint involved a list of these weird coördinates. After a little research the only thing I found is that they might have to do something with UTM.

It reads

  49510.24298m; 65.72307667 degrees

It should lead to somewhere in the Netherlands in the western part of the province: Gelderland. Can any of you guys help me out with this?


texts: Joining two datasets based on matching values in multiple fields

texts:
I am trying to find a way to join CSV data to a vector layer only when they both share matching values in multiple fields, rather than just one field.

This is a problem I'm running into with multiple projects, but I'll just describe one here as an example:

The data I am working with is as follows:

a shapefile layer containing a list of all U.S. towns/cities with
population over 10,000.
49 CSV files (one for each of the lower 48 states + 1 for D.C.) from the Bureau of Justice Statistics (BJS), containing data on violent crime rates by city.

I am trying to join the data from the BJS .csv files onto the "U.S. Towns and Cities" layer so that each town/city feature contains data for violent crime in that city.

The problem is that the the "Joins" interface will only let me select one field to join on. So for instance, when go to "US Towns & Cities" --> Properties --> Joins and click "+" to create a new join, this is what I have (as of QGIS version 2.12):

So here I am joining the BJS data for Georgia to the "U.S Towns and Cities" layer by city name ... now assume that I've done the same for the BJS Florida data.

Now, there are several ways that this causes trouble. The most significant problem is that there are name collisions for cities in different states. For instance, the "U.S. towns & cities" map layer contains 8 cities named "Gainesville", ALL of which now have the BJS data for both Gainesville, Georgia AND Gainesville, Florida associated with them:

I need to only join the crime data for for a city when both the city name AND the state name coincide so that Gainesville, Florida's data only gets associated with Gainesville, Florida and not Gainesville, Georgia  ... How do I do this?

Also, another problem caused by the above approach, which isn't necessarily causing incorrect results, but seems "wrong" is that every single city now contains nearly 1000 NULL fields (~20 of them for each of the 49 CSV files). ... I feel like there must be a better way to do this.


texts: FieldPyculator Plugin - How to concatenate fields?

texts: How can i concatenate two fields using pyculator plugin? There is any document with all pyculator functions that i can use in "field expression"?


texts: Is dynamically changing a datasource possible in arcgis js api?

texts:
I've posted a related question is in the ArcGIS Desktop forums, but hasn't gotten any replies yet. In that case, I would like tackle this in another way, which is dynamically changing a datasource using arcgis js api. Is this possible?

To give some detail, the layers defined in my mapservice are joined with another view through oledb. Will the class QueryDataSource be able to do this, and how?


texts: Migration in PostgreSQL 11 from PostGIS 2.4 to 2.5 fails

texts:
This is my attempt to migrate the data:

ALTER EXTENSION postgis UPDATE;
WARNING:  'postgis.backend' is already set and cannot be changed until you reconnect
ERROR:  operator does not exist: gidx public.&& geography
LINE 1: SELECT $2 OPERATOR(public.&&) $1;
                  ^
HINT:  No operator matches the given name and argument types. You might need to add explicit type casts.
QUERY:  SELECT $2 OPERATOR(public.&&) $1;

More details:

CentOS 7
select postgis_full_version();

POSTGIS="2.4.7 r17330" PGSQL="93" (procs need upgrade for use with "110") GEOS="3.7.1-CAP
I-1.11.1 27a5e771" PROJ="Rel. 4.9.3, 15 August 2016" GDAL="GDAL 1.11.4, released 2016/01/2
5" LIBXML="2.9.1" LIBJSON="0.11" (core procs from "2.4.4 r16526" need upgrade) RASTER (ras
ter procs from "2.4.4 r16526" need upgrade)

I think maybe I'm affected by this bug.

UPDATE:

I think it is relevant that I've migrated from PostGIS 2.1.8 in Postgres 9.3 to version 2.4.7.  Then migrated to Postgres 11 and PostGIS 2.4.7, and finally tried to version 2.5.2.


texts: How to get perpendicular line script to work?

texts:
I'm trying to get this script to work by @gene but I'm not sure how to load the results into QGIS. I tried to make a memory layer but not sure how to add the geometries of the point, line and the intersection point:

point_layer = QgsVectorLayer("Point", "point_layer", "memory")
pr = point_layer.dataProvider()
fet = QgsFeature()
fet.setGeometry(QgsGeometry.fromPoint(QgsPoint(30,120)))

line_layer = QgsVectorLayer("LineString", "line_layer", "memory")
lr = line_layer.dataProvider()
seg = QgsFeature()
line_start = QgsPoint(50,50)
line_end = QgsPoint(100,150)
line = QgsGeometry.fromPolyline([line_start,line_end])

def intersect_point_to_line(point, line_start, line_end):
    magnitude2 = line_start.sqrDist(line_end)
    u = ((point.x() - line_start.x()) * (line_end.x() - line_start.x()) + (point.y() - line_start.y()) * (line_end.y() - line_start.y()))/(magnitude2)
    ix = line_start.x() + u * (line_end.x() - line_start.x())
    iy = line_start.y() + u * (line_end.y() - line_start.y())
    return QgsPoint(ix,iy)

line = QgsGeometry.fromPolyline([point,intersect_point_to_line(point, line_start, line_end)])

I am using QGIS 2.10.


texts: How to insert reversed map with Google Maps API?

texts: Is there a way to put map with Google Maps API to a website, but upside down? North must be down, and south up. Is there any possible way to do it? I really need to put map that way and I need really fast solution.


texts: GeoServer doesn't recognize imagery in GeoPackage

texts:
I'm currently trying to upload a GeoPackage containing raster images to a GeoServer (v2.13) instance with the GeoPackage extension. The file contains ~425 layers but after I add the file to the server's data directory and attempt to add it to the data stores (GeoPackage Mosaic), I get a vague error message:

Could not list layers for this store, an error occurred retrieving them: null

I'm wondering if this is an existing error somewhere but I can't seem to find anything current in the GitHub repo issues or elsewhere on the net. The images are perfectly in QGIS 3 and in Python.

See the screenshot below (ignore the locale setting weirdness):


texts: ArcMap Georefencing - sheet/image problem

texts:
I have a problem and I wanted to know if anyone else had experienced this. I think it is some type of bug as it doesn't occur all the time.

The problem is when I am georeferencing some 'tiff' images, when I zoom in to the image Arcmap seems to rotate the image without my consent. Even when the image has actually been georeferenced if I zoom in it will sort of flip the image, but if I zoom out of scale it will go back to the desired location.


texts: How to add features to geometry collection using .ShapeCopy

texts:
I'm having an issue with trying to complete the loop of adding the selected geometry to my geometry collection. The goal is finding streets that intersect the given street name (then locating the intersection, which is handled later in code).

Here is how I'm currently handling it and maybe someone can suggest an easier alternative or provide some help with the issues I'm having in my 'TODO:' statements.

// get the "Street" featureclass
IFeatureLayer pFeatureLayer = new ESRI.ArcGIS.Carto.FeatureLayer();
pFeatureLayer.FeatureClass = m_fcToQuery;

// query and get a feature cursor of selected streets
QueryFilter pQueryFilter = new ESRI.ArcGIS.Geodatabase.QueryFilter();
pQueryFilter.WhereClause = "TownCode = " + townCodeVal + " AND StreetName = '" + streetNameVal + "'";

// merge selected street segs into geometry collection 1
IGeometryCollection pGeomColl1 = new PolylineClass();
IGeometry pGeometry1 = new ESRI.ArcGIS.Geometry.IGeometry();
pGeometry1 = pGeomColl1 as IGeometry;

IGeoDataset pGeodataset = new ESRI.ArcGIS.Geodatabase.IGeoDataset();
pGeodataset = pFeatureLayer.FeatureClass as IGeoDataset;
pGeometry1.SpatialReference = pGeodataset.SpatialReference;
IFeature pFeature = new ESRI.ArcGIS.Geodatabase.IFeature();
IFeatureCursor pFeatureCursor = new ESRI.ArcGIS.Geodatabase.IFeatureCursor();
pFeature = pFeatureCursor.NextFeature();
if (pFeature == null) {
    return null;
}
while (!(pFeature == null)) {
    // TODO: Fix this
    pGeomColl1.AddGeometryCollection(pFeature.ShapeCopy); <--error
    pFeature = pFeatureCursor.NextFeature();
}

// get a feature cursor of the int streets (3rd dropdownlist)
pQueryFilter.WhereClause = "TownCode = " + townCodeVal + " AND StreetName = '" + streetsInt + "'";
pFeatureCursor = pFeatureLayer.FeatureClass.Search(pQueryFilter, false);

//Merge the selected streets into the Second Geometry Collection
IGeometryCollection pGeomColl2 = new ESRI.ArcGIS.Geometry.PolylineClass();
IGeometry pGeometry2 = pGeomColl2 as IGeometry;
pGeometry2.SpatialReference = pGeodataset.SpatialReference;
pFeature = pFeatureCursor.NextFeature();
if (pFeature == null) {
    return null;
}
while (!(pFeature == null)) {
    // TODO: fix this
    pGeomColl2.AddGeometryCollection(pFeature.ShapeCopy); <-- error
    pFeature = pFeatureCursor.NextFeature();
}

//Find the Intersection of the Two Streets
ITopologicalOperator2 pTopoOp = pGeometry1 as ITopologicalOperator2;
Polygon pBuffer1 = pTopoOp.Buffer(0.1) as Polygon;
pTopoOp = pGeometry2 as ITopologicalOperator2;
Polygon pBuffer2 = pTopoOp.Buffer(0.1) as Polygon;
pTopoOp = pBuffer2 as ITopologicalOperator2;

// TODO: fix this
IGeometry pintersection = pTopoOp.Intersect(pBuffer1, 4) as IGeometry; <-- error

if (pintersection.IsEmpty) {
    return null;
}


texts:
I am creating a Python script in QGIS 3 and I want to use the processing module. I want to transform my raster file to a vector file with the GRASS algorithm r.to.vect. My problem here is that the output I am getting is always empty. The way I am running the algorithm looks like this:

    param = {"input": path_to_some_raster_file,
             "type": 2,
             "output": some_tmp_file,
             'GRASS_OUTPUT_TYPE_PARAMETER': 0
             }

    vector = processing.run('grass7:r.to.vect', param)

But if I use the same parameters in the processing toolbox then the desired output is generated. What am I doing wrong with my Python attempt?


texts: ggmap openstreetmap SCALE, no zoom necessary?

texts:
I am relatively new to map projections but believe I have understood how one could translate lat/lon to the map tile coordinates X and Y, which is nicely explained here:

http://wiki.openstreetmap.org/wiki/Slippy_map_tilenames

When requesting Google static maps, we need to supply a zoom value which makes sense to me.
However for the function get_openstreetmap(), I am quite puzzled on the scale parameter and how one could possibly translate scale to zoom. That seems impossible/difficult since zoom values are discrete integers and scale more or less continuous.
If scale->zoom is not done, how does the package project lat/lon to Y/X ?
I dug through the sources but the only function I could find was LonLat2XY() which also requires a zoom level.

What function performs the translation for OSM maps?


texts: How to pass parameters to PostGIS queries containing geography type in Python using psycopg2

texts:
I have the following example (adapted from the PostGIS manual (version 2.2):

CREATE TABLE global_points (
    id SERIAL PRIMARY KEY,
    name VARCHAR(64),
    location GEOGRAPHY(POINT,4326)
  );

INSERT INTO global_points (name, location) VALUES ('London', ST_GeographyFromText('SRID=4326;POINT(0 49)') );
INSERT ..........
c = conn.cursor()
try:
  query ="SELECT name FROM global_points WHERE ST_DWithin(location, ST_GeographyFromText('SRID=4326;POINT(-110 29)'), 1000000);"
  c.execute(query)
  l = c.fetchall()

.................
How shall I modify the query to allow for specifiyng the variables lon and lat in POINT in order to be able to iterate over a set of geographic locations?


texts: How to de-stripe a Landsat 8 band with python?

texts:
I am working with Landsat 8 images and some of the bands have striping, e.g.:

I've been going down a rabbit hole by looking for a solution for this. Does anyone know a python library that could be used for this?


texts: Is it possible to edit geographic data over the web in QGIS cloud or QGIS server?

texts:
I'm looking for something similar to the ArcGis Server example but based on QGIS.

Could anybody tell me where I can get more information about functionalities like this?


texts:
I need to publish raster and vector data through GeoServer WMS service and I have already chosen PostGIS as datasource to keep vector data.

Now I would like to know what is the best practice to work with raster data.
I see two options:

store raster data in PostGis and then use ImageMosaicJDBC or Postgis Raster Plug In (I read about some limitations and instability)
store raster data somewhere on file system and setup GeoServer to read data from there.

I have little experience with GIS and I wonder if someone could suggest the best choice?


texts:
I need an approach to combine points close one to another into a common point.

I can take an average of points' {lat, long} to create one point.

How I can do this with ArcMap or Python? I have an Excel file for those points.

Should I combine according to {lat, long} information with Excel?


texts: Getting geometry from multiple tables using PostGIS?

texts:
I am very new to PostGIS, PostgreSQL, and SQL in general.

I have 44 tables in my PostGIS database, and each one represents a different layer of vector data. Each was loaded from a separate shapefile, and each has a column that describes the geometry for that layer, called wkb_geometry

I want to select a specific polygon on one layer, and then retrieve ALL the geometry from a subset of layers that overlap that polygon's bounding box. I'm not feeling picky about the order coming out, but it would be helpful if it was organized by the tables each group of geometry came from.

Here's a sample of my SQL statement:

SELECT
    ST_AsEWKT(wkb_geometry) /* Some of this data has z values, hence the EWKT */
FROM
    table2, table3, table4, table5
WHERE
    wkb_geometry &&
        (
        SELECT
            wkb_geometry
        FROM
            table1
        WHERE
            ogc_fid = 25
        );

which returns an error:

column reference "wkb_geometry" is ambiguous
/* note that every table that I am selecting from has a "wkb_geometry" column */

What is the right way to do this?


texts: Unable to retrieve custom extension state from a document in ArcEngine

texts:
I created a custom extension in ArcMap 10.1 which implements the IPersistVariant interface.  State related to the extension features is persisted and retrieved via the Load and Save methods.  Everything works fine in ArcMap.

I am now building an ArcEngine app which needs to load documents which have been modified with the custom extension in ArcMap.  I need to be able to load in the extension, so that when I load the document, the extension specific state is retrieved.  I have no problem loading the extension in ArcEngine as follows:

        Type extensionManagerType = Type.GetTypeFromProgID("esriSystem.ExtensionManager");
        System.Object extensionManagerObject = Activator.CreateInstance(extensionManagerType);
        this.comReleaser = new ComReleaser();
        comReleaser.ManageLifetime(extensionManagerObject);
        this.extensionManager = extensionManagerObject as IExtensionManager;
        IExtensionManagerAdmin extensionManagerAdmin = this.extensionManager as IExtensionManagerAdmin;

        // Add the PirsaExtension to the extension manager to load it.
        UID uID = new UIDClass();
        uID.Value = "PirsaX.Extension.PirsaToolsExtension";
        extensionManagerAdmin.AddExtension(uID, this as IBatchApplication);
        IExtension extension = this.extensionManager.FindExtension("PirsaMappingToolsExtension");
        // Remember the extension manager as a PirsaExtension so that later we can use it to load
        // a document.
        PirsaToolsExtension pirsaExtension = extension as PirsaToolsExtension;
        this.batchMapGenerator = pirsaExtension.BatchMapGenerator;
        // Enable the extension.
        IExtensionConfig extensionConfig = extension as IExtensionConfig;
        extensionConfig.State = esriExtensionState.esriESEnabled;

Doing this loads the extension properly.  In the debugger I see its Startup method being called and see the extension being enabled.

The problem arises when I open the document.  The extension's Load(IVariantStream stream) method never gets called.  Here is how I open the document.

        mapDocument = new MapDocumentClass();
                if (mapDocument.get_IsMapDocument(filePath))
                {
                    mapDocument.Open(filePath, String.Empty);
        }

There is a pause when the mapDocument loads, and then I can access its details, but none of the custom state is there.  How do I tie opening the document into the extension?

Any advice or insights would be greatly appreciated.

Tom Cuthill


texts:
I have been using mapinfo for some time now but I have not been able to complete the following task and I would very much welcome some help if anyone can offer it to me:

I have data which I need to link or join to railway lines or sections of adjacent railway lines. From this I would like to complete hot- spot lines of routes similar to thematic mapping to show troublesome routes.

I can do this using drivetime when there is a start X and Y and end X and Y point for the train data but the data I am working on does not have these. Instead I only have one geocoded data point which does fall on the railway line but I can not join it to the railway line. I would then need to thematically map these sections to show different widths for larger numbers.

I hope this makes sense to at least one of you.