texts: Arcpy: Geometry object - distanceTo method in a for loop

texts:
I'm new to Python so sorry for this maybe trivial question.

I have 4 lines and 4 points, which belong together. I want to get the closest distance from point 4 to line 4, point 3 to line 3, ... . I want to use a SearchCursor and Geometry methods (distanceTo), no geoprocessing tool like "Near".

The code is used so far is:

    import arcpy

    point_feature = r"C:\Users\name\Documents\ArcGIS\Default.gdb\Gas_points"
    line_feature = r"C:\Users\name\Documents\ArcGIS\Default.gdb\Gas_lines"

    pointcursor = arcpy.da.SearchCursor(point_feature, ['SHAPE@'])
    linecursor = arcpy.da.SearchCursor(line_feature, ['SHAPE@'])

    for row in pointcursor:
        geometry = pointcursor[0]
        print geometry.type
        for row in linecursor:
            newgeometry = linecursor[0]
            print newgeometry.distanceTo(geometry)

The result i get:

point
6.47818694862
67.4681286231
66.0830276147
67.4681286231
point
point
point

Should be a logical mistake in the for loop, I guess.


texts: Graduated symbology shows borders between matching categories

texts:
I'm using graduated symbology on a vector layer in QGIS and i do not want lines between polygons of the same color. I tried making the border outline transparent, and also setting the border outline to match the fill color and set the border width to 1. I can still see a thin grey border regardless of the settings.

I should add, i'm trying to do this in python. This worked in QGIS but not python: Styling border to match fill in QGIS?

for actmin, actmax, myLabel, myColor in ranges:
    mySymbol = QgsSymbolV2.defaultSymbol(maplayer.geometryType())
    mySymbol.setColor(myColor)
    mySymbol.setAlpha(myOpacity) mySymbol.symbolLayer(0).setDataDefinedProperty('color_border','color_rgba'+str(rgbacolors[i]))
        myRange = QgsRendererRangeV2(actmin, actmax, mySymbol, myLabel)
        myRangeList.append(myRange)
        i+=1

myRenderer = QgsGraduatedSymbolRendererV2('', myRangeList)
myRenderer.setMode(QgsGraduatedSymbolRendererV2.Quantile)
myRenderer.setClassAttribute(mappingcol)
maplayer.setRendererV2(myRenderer)


texts: Multi layer interactivity not working with carto.js

texts: I've two sublayers with interaction set as true and the proper columns added in the interactivity method. The featureClick listener is attached to the layer. On the one hand, popups are working on both sublayers independently when they are alone. But when both sublayers are present only the top-sublayer is triggering the featureClick event. How can I launch both popups from the two sublayers?


texts: Merging neighbor polygons up to a specific area size in Python

texts:
I am looking into a way to merge neighboring polygons in Python into 10 groups so that each group has a specific area size. At the beginning I have 389 micro-zones stored in outmulti which I want to split into 10 neighboring (of equal size) zones. Assuming that the whole area of the 389 micro-zones is 19447273 m², I want to get 10 (being neighbors) bigger zones of size around 19447273/10 m² each.
My variable which stores the polygons is outmulti in the form:

Type
Value

geometry.polygon.Polygon
Polygon object of shapely.geometry.polygon module

geometry.polygon.Polygon
Polygon object of shapely.geometry.polygon module

geometry.polygon.Polygon
Polygon object of shapely.geometry.polygon module

geometry.polygon.Polygon
Polygon object of shapely.geometry.polygon module

I know I should use unary_union but it will group together into a single geometric feature of all my polygons (which I don't want). I do not know how to do that so that specify the neighboring polygons, as well as the limit of the area size of each of the 10 (bigger) zones.


texts: How to change dimension annotation style in ArcMap 10.1?

texts: For some reason, the option to add a "New" or "Import" style is grayed out. I tried creating a different dimension annotation feature class with the style I want, and then copy and pasting the records...but the option to copy was grayed out. I have a ton of dimension annotations in the feature class I am trying to change so re-doing it is not an option. Anyone know what is going on?


texts:
I am aware that a similar question has been asked as Sites with recreational shapefiles for United States? but I am looking for comprehensive sources of recreational data (specifically for the province of Ontario, Canada).

Land Information Ontario has some data (I have everything they can provide), but I am looking for something covering more categories (and fairly current if possible).

I'm hoping to find golf courses, campgrounds, hiking/biking trails, kayaking/canoeing routes and other data of this nature.

Other categories are welcome if you have any suggestions.


texts: If/Else Statement in Python ArcGIS Field Calculator

texts: I have a shapefile and want to calculate a field based on another field. The two fields are "DURATION" and "STATUS". I need "STATUS" to display three different outputs based on "DURATION". I want "STATUS" to display a "1" if "DURATION" is greater than 31. I want "STATUS" to display a "2" if "DURATION" is less than or equal to 31 but greater than 0. I need "STATUS" to display "3" if "DURATION" less than or equal to 0.


texts: ArcMap calculate distance (from and to) from a points shape to the lines shape table

texts:
I am trying to calculate the From and To fields in the line table.
I have got a points shape and i have got a line shape. ( see for example the figure bellow.)
The points( in red ) intersect with lines.
My plan is to get the distance from points for the begin and end of the each line segment. I have tried spatialjoin but it gives only the distance but not giving +/- form the points location. I need to update each line with the 'ID'of the point + begin and end distance from that 'ID'point ( see the figure : in table i have updated 2 records manually). I think i am clear about what i want to do.
What i have done is just measured the distance with measuring tool and updated the fields manually.
I have more than 500 line segemnts.
I think there must be a better way than this.
I am working with ArcMap 10.6

figure ( it's an example) :


texts: Retaining polygon feature attributes in points randomly generated within using QGIS

texts:
How do I  retain polygon feature attributes in random points generated within these polygons using QGIS?
Like Generating random points inside polygons that capture values from multiple attribute fields, I found that QGIS 3.14 "Random Points inside Polygons" tool does not retain any attributes from the polygon features.


texts:
I'm trying to understand the math behind this transformation. I should note that I've been able to perform the transformation using GeoPandas, but I'd like to understand what is going on under the hood.

I've got this dataset of ZIP code polygons from NYC open data, that appears to be using NAD83 / New York Long Island (ftUS) (EPSG:2263), and I'm trying to transform it to WGS84 (EPSG:4326).

How do I turn, for example, 986038.6614370644, 213051.06312088 tuple into -73.99354429600449, 40.75145227487899?


texts: specify coordinate system on openlayers Filter

texts:
Is there anyway to specify which coordinate system to use with a filter?

using the following filter:

new OpenLayers.Filter.Spatial({
    property: "boundsText",
    type: OpenLayers.Filter.Spatial.INTERSECTS,
    value: geom
})

Produces this HTTP POST to my GeoServer:

<ogc:intersects>
<ogc:propertyname>bounds</ogc:propertyname>
<gml:polygon xmlns:gml="http://www.opengis.net/gml" srsName="EPSG:4326">
    <gml:outerboundaryis>
        <gml:linearring>
            <gml:coordinates decimal="." cs="," ts=" ">-7670608.6611,2132898.8368875 -6261721.356,254382.4300875 -6653078.94075,6633511.0615125 -12758257.26285,6790054.0954125 -7670608.6611,2132898.8368875</gml:coordinates>
        </gml:linearring>
    </gml:outerboundaryis>
</gml:polygon>
</ogc:intersects>

and i get the following response error from GeoServer:

  -7670608.6611 outside of (-180.0,180.0)

How can i change whatever coordinate system is being sent in my filter to decimal degrees?


texts:
Does anyone have a bit of code to plot a multiline?  I just want to add some IPoints to a collection and display in a particular color and width, but I can't find a good way to do this.

Right now I'm using this after a dynamic draw event:    dynamicDisplay.DrawMultipleLines(pointcollection);
But there's no clear place to set the graphic attributes, and it plots incorrectly.


texts: how to open geotiff as base layer on openlayers

texts:
i want to add a geotiff image as base layer on openlayers :

copy my image.tiff in a directory and open it with openlayers directly without using a cartography server.


texts:
This is probably going to annoy the more experienced map people here, but I'm sort of new to all this, or at least I don't have a lot of experience.
I know this is a bit strange and not ideal, but I have a huge SVG file of geographical shapes and regions.  I want to be able to plot points and draw shapes in Openlayers.  The problem is that, of course, this file is over 50 MB, so I can't exactly expect users to load the whole thing then navigate around.
I have a less detailed map for when users are zoomed out, so I'd like to use this map as they zoom in, but then tiling would of course be best.  I don't need much detail beyond what I have so using just these two maps as starting points is fine.  Like I said, I was planning on using OpenLayers since I've used it in the past but I'm not tied to it.
How is the best way for me to proceed?
Thanks.


texts:
I noticed that I can't use map() to change the system:index of a FeatureCollection in GEE. I can change one system:index at a time. And I can use map() to edit other properties. Is there something special about 'system:index' that prevents me from using map() to edit?
var features = [
  ee.Feature(ee.Geometry.Rectangle(30.01, 59.80, 30.59, 60.15), {name: 'Voronoi'}),
  ee.Feature(ee.Geometry.Point(-73.96, 40.781), {name: 'Thiessen'}),
  ee.Feature(ee.Geometry.Point(6.4806, 50.8012), {name: 'Dirichlet'})
];

var fc = ee.FeatureCollection(features);

// Rename one feature
var newName = fc.first().get('name');
var renamed = fc.first().set('system:index', newName)
print('I can rename the system:index for one feature', renamed)

// Map across FC
var rename = function(f){
  var myNewName = ee.String(f.get('name'))
  return f.set('system:index', myNewName)
}
print('But cannot for a collection', fc.map(rename).first())

// Mapping a different function works ok
var addProperty = function(f){
  var newProperty = ee.String(f.get('name')).cat('Hello')
  return f.set('newProp', newProperty)
}
print('I can map set other properties', fc.map(addProperty).first())

I see this question found a workaround using ee.List. Is this the only solution for this?
https://code.earthengine.google.com/?scriptPath=users%2Fcaseyengstrom%2Freprex%3AsetSystemIndex


texts: Georeferencing points to an aerial image

texts:
I have a file containing a large number of points that were digitized from scanned plans in CAD.  The points represent pit lids at an industrial facility.

When the point file is laid over the aerial image the points are all out by varying amounts.

Around the perimeter of aerial I can see a decent number of the pits. Within the centre of the facility there is a lesser amount of pits visible.

Is there a technique or workflow to "transform" the pits to the aerial in QGIS?

I inherited the files in this manner and would do it differently if I had the chance


texts: MapServer reading SQL from 8-bits data type TIFF

texts:
I'm converting to SQL (PostGIS) some raster files (TIFF). Then I'm using those SQL files for queries with MapServer, and displaying it in my app with OpenLayers with no problems. Everything goes well. The TIFF files are like 600KB weight (by the way, the SQL files are like 22MB approx), and I can see with the QGIS explorer that they are Float64 as "Data Type" and -1.7e+308  as "No Data Values" (excuse me if the translations isn't OK, I'm using QGIS in Spanish). As I said, everything is OK.

In MapServer, I have a map file like this:

LAYER
    NAME zonificacion_cabernet_sauvignon_actual_rango5

    TYPE raster
    STATUS ON
    DATA "PG:host=192.168.1.133 port=5432 dbname='XXXXXX' user='postgres' password='XXXXXX' schema='raster' table='XXXXXX' mode='2' where='ano=2014 AND rango_5' "

    PROJECTION
        "init=epsg:32719"
    END
    PROCESSING "SCALE=AUTO"
    OPACITY 80
    CLASS
        NAME ""
        EXPRESSION ([pixel] = 1)
        COLOR "#3399cc"
    END
    CLASS
        NAME ""
        EXPRESSION ([pixel] = 2)
        COLOR "#99cc33"
    END
    CLASS
        NAME ""
        EXPRESSION ([pixel] = 3)
        COLOR "#cc3399"
    END
END

Everything goes well, I can see 3 colors (each for those cases) in the map of my app.

But now, my client, who provides this TIFF files, made a change, and he told me "I will use Integer instead Float". So, I started to replace the SQL files, doing all the previous actions. The TIFF files now are much lighter, they are around 60-70KB and the SQL files like 2MB, but the script on the map file is not working, but if I use this:

EXPRESSION it does:

CLASS
    NAME ""
    EXPRESSION ([pixel] = 85)
    COLOR "#3399cc"
END
CLASS
    NAME ""
    EXPRESSION ([pixel] = 170)
    COLOR "#99cc33"
END
CLASS
    NAME ""
    EXPRESSION ([pixel] = 255 )
    COLOR "#cc3399"

I checked with QGIS Explorer the new TIFFs, and I can see that "Data Type" is now "Byte - Natural Integer of 8 bits" and "No Data Values" are 255

What should I do to my client in order to have TIFF files with correct data? because I need to match the EXPRESSIONwith 1, 2 and 3 (the amount of pixels for each color) but not 85, 170 and 255.


texts: Getting Python error for raster analysis interpolation using nearest neighbor in QGIS

texts:
I am trying to interpolate temperature values into a raster using QGIS 2.18.15 in Windows. I need the output file as a GeoTIFF (.tif) file. I have my input for temperature already in a shapefile, but whenever I try to execute the command, it shows a Python error.

These are the procedures I have taken

But it shows me this

When I click Stack trace it says this:

An error has occurred while executing Python code:

UnicodeDecodeError: 'ascii' codec can't decode byte 0x83 in position 32: ordinal not in range(128)
Traceback (most recent call last):
  File "C:/PROGRA~1/QGIS2~1.18/apps/qgis/./python/plugins\GdalTools\tools\widgetPluginBase.py", line 88, in onFinished
    self.base.onFinished(exitCode, status)
  File "C:/PROGRA~1/QGIS2~1.18/apps/qgis/./python/plugins\GdalTools\tools\dialogBase.py", line 208, in onFinished
    msg = unicode(self.process.readAllStandardError())
UnicodeDecodeError: 'ascii' codec can't decode byte 0x83 in position 32: ordinal not in range(128)

Python version: 2.7.5 (default, May 15 2013, 22:44:16) [MSC v.1500 64 bit (AMD64)]
QGIS version: 2.18.15 Las Palmas, ac2d111ed1

The message log says this:

2018-01-18T19:04:55 1   Traceback (most recent call last):
              File "C:/PROGRA~1/QGIS2~1.18/apps/qgis/./python/plugins\GdalTools\tools\widgetPluginBase.py", line 88, in onFinished
                self.base.onFinished(exitCode, status)
              File "C:/PROGRA~1/QGIS2~1.18/apps/qgis/./python/plugins\GdalTools\tools\dialogBase.py", line 208, in onFinished
                msg = unicode(self.process.readAllStandardError())
            UnicodeDecodeError: 'ascii' codec can't decode byte 0x83 in position 32: ordinal not in range(128)

What I should do?

Here is a link to my shapefile
https://drive.google.com/file/d/17At2G7yBLRxSrWhS9yhIKIJBs8R0iHd_/view?usp=sharing


texts:
I created a field to display the whole route of several photos in a point layer in order to be able to visualize them on the map. I try to do so by using this expression:

concat('C:\Users\DIDAC\Documents\CONSULT\EIA OLIVERES CLARAMUNT\treball camp\fotos\IGN\',"file_name")

"file_name" being the name of every photo.

Here's the error message: "syntax error, unexpected Unknown_CHARACTER"


texts: Displaying raster tiles from OSM data stored in a postgis server

texts:
I have been recently tasked with setting up a raster tile server that serves tiles from data pulled from a postgis store.

I have managed to get a vector tile server up and running but I'm stuck on the raster tile server bit. This is what the raw data looks like (an example row, it's basically a list of all the streets in Switzerland):

{ id: 1537021,
  osm_id: '317400638',
  osm_name: null,
  osm_meta: null,
  osm_source_id: '3236568222',
  osm_target_id: '3236568210',
  clazz: 72,
  flags: 6,
  source: 1216745,
  target: 1216746,
  km: 0.048946153,
  kmh: 10,
  cost: 0.004894615,
  reverse_cost: 0.004894615,
  x1: 8.1708598,
  y1: 47.4290024,
  x2: 8.1707774,
  y2: 47.4285814,
  geom_way: '0102000020E610000003000000F7668AEF7A572040E8F3F68CE9B64740AEC2C1397E572040DBAC5516E0B64740B5CAA7227057204042565CC1DBB64740' }

What would be the best approach to generate raster tiles given this data?

I managed to cook up vector geoJSON data via the following:

SELECT osm_name As properties, ST_ASGEOJSON(geom_way) as geometry
FROM ch_2po_4pgr
WHERE ST_Intersects(ch_2po_4pgr.geom_way,
ST_SetSRID(ST_MakeBox2D(ST_Point($1, $2), ST_Point($3, $4)), 4326))

Most tutorials seem to deal with using Mapnik and a shapefile datasource but what I have is this postgis table :(


texts: Layer Stacking Sentinel-2 image on Google Earth Engine

texts:
How can I stack sentinel 2 Image?
My complete code is :
 var image = ee.ImageCollection('COPERNICUS/S2')
.filterBounds(table)
.filterDate('2019-06-01','2019-07-01')
.median()
var subset = image.clip(table)


texts: How to set up GeoDjango with PostGIS 2.0?

texts:
I would like to run GeoDjango (Django 1.4) with Postgis 2.0.

Whenever I try to run syncdb I get the following error for geographic fields:

Failed to install index for stores.Store model:
operator class "gist_geometry_ops" does not exist for access method "gist"

I have found various online discussions about this, but no clear answers.

It sounds as though I need to create some indexes by hand... or perhaps there is a patch that can help, but I don't know how to add it.

Please could someone explain simply how to fix this error?


texts:
I'm trying to plot several polygons on a raster. Several posts suggest to use matplotlib and descartes like this one: Plot shapefile on top of raster using plot and imshow from matplotlib

My polygons do not appear on the map although they have the same coordinate system and that the polygon coordinates are within the bbox of my raster. The plots below show the raster and the polygons as separate plots (a and b) and one can see that a) has "window coordinates" in pixels and b) in EPSG: 4326.

How can I plot my polygons on the rasterio window?

from descartes import PolygonPatch
import matplotlib as mpl
import rasterio

# Import .tif image (RGB)
src = rasterio.open(base_tif)  # src.crs = CRS.from_epsg(4326)
# Select window based on 4 coordinates
rst = src.read(window=windows.from_bounds(minx, miny, maxx, maxy, src.transform))  # 7.493877446690874, 7.545662625573129, 45.95562555558752, 45.99882089053693, Affine(0.00011273472018020788, 0.0,6.7805624095004235,0.0, -0.00011273472018020788, 46.66875202450472)

# My polygons
lakes_for_gl = [{'type': 'Polygon',
  'coordinates': [[[7.559165, 45.965006],
    [7.558827, 45.965006],
    [7.558827, 45.965118],
    [7.559277, 45.965231],
    [7.559165, 45.965006],
    [7.559165, 45.965006]]]},
 {'type': 'Polygon',
  'coordinates': [[[7.51362, 45.958918],
    [7.513056, 45.958918],
    [7.513056, 45.959256],
    [7.513733, 45.960271],
    [7.514296, 45.960496],
    [7.514409, 45.959594],
    [7.51362, 45.958918],
    [7.51362, 45.958918]]]}]

fig, ax = plt.subplots()
rasterio.plot.show(rst)
ax = plt.gca()

patches = [PolygonPatch(feature, edgecolor="red", facecolor="red", linewidth=2) for feature in lakes_for_gl]
ax.add_collection(mpl.collections.PatchCollection(patches, match_original=True))
ax.axis('scaled')
plt.show()


texts: Identify lat long points falling out of cluster

texts:
Problem Statement

I have a set of lat-long locations, say 30 in number. Ideally they should be just a few meters from each other on ground. But sometimes there's that odd junk data point that's 20kms away or so, which messes up the whole cluster. I want to programmatically identify these "odd ones out". What would be a good way of going about this? Will mention here: I have to repeat this thousands of times, so seeking a programmed way.

Strategies coming to mind

Use Sklearn package's AgglomerativeClustering as done in the Kaggle example below.
Draw distance buffers and dissolve, then find largest polygon and all points outside of it, based on PyQGIS answer linked below.
Convert lat-longs to OpenLocationCode / plus.code, and find the most frequent one. This is simplistic and crude as the plus code grid is predefined (what if two points are right next to each other? what if one  grid size resolution is too large and next level one that's 400x smaller is too small?)

Specifics

Max cluster area : 2 sq km
Or, max distance between any two points: 0.5 km.
Input would be one dataframe having one point per row
Output would be two dataframes or lists telling which is 'Out' and which is 'In'.

Good leads to pursue

Kaggle - Basic Visualization and Clustering in Python
Filtering out spatial-outliers from cluster of points using PyQGIS? - answer is promising but I can't use it directly : Will have find largest polygon and then find all the points in that. If would be great to adapt this to work directly with pandas/lists instead of having to convert into a shapefile layer and all.
How to cluster points and plot - It's in R so a python adaptation would make for a good answer here.

I'll update here if I find more leads.

Clarification notes

No weighing, scoring business etc. They're just lat-long points.
The data is ok enough that we can trust the majority; most points are near the right place and there may be just 1 to 3 outlaws in each set.
Preferred programming language is python3, but I'm open to seeing other solutions; will adapt.
I'm not too concerned about lat-long contorting (the work area is in the tropics so a lat-long grid is pretty squarish), so it's possible to boil this down to a simple x-y scatterplot, and finding the odd one out.
I'm aware that the algorithms can even generate multiple clusters. In that scenario I'll just take the most populated cluster and all the other points will be branded outlaws. In case of 50-50, I'll flip a coin :)
No need of visualising! This has to be a sub-routine at back end.

Related questions on this forum

Creat clusters using long and lat - basic question, could be great as a starting point, unfortunately no answers yet.
How to calculate/measure the spread of coordinates in an area using python - no solid answer
Clustering geographical data based on point location and associated point values - went into more complicated territory
Clustering of Spatial Data in R or Python - closed for being unclear


texts:
I need some clarification. I have a map the is set to the scale of 1:6000. In Settings>Project Properties, under the General tab I have the Layer Units set to feet. Therfore I assume that means the scale is 1 foot equals 6000 feet. Does this seem rights? Or is it 1 inch equals 6000 inches?

I am ultimately attempting to convert this to traditional engineering scale of 1 inch equal __ feet.

Found this handy description for anyone that interested. http://groundwater.fullerton.edu/groundwater/Maps,_Scale,_GIS_and_GPS/Guide_to_Map_Scale.html


texts: Same data displayed in different locations in ArcGIS Explorer and QGIS

texts:
I've been working on some GIS data in middle eastern countries.
A client gave me the data in .gdb format. So I thought it could be loaded in QGIS.

However, when I load the data with QGIS, they are located in a totally different place (well, to be more precise, the longitude seems correct but the latitude comes around the equator).

The strange part is when I load the same data with ArcGIS Explorer, they are on the correct location. So it's not the problem with the original data, but with the way I load them.

I don't have ArcDesktop, so my destination is to be able to deal with the data on my QGIS.

I tried ogr2ogr conversion from gdb to shp but it just gives the same result, which means, when I load the resulted shapefile it goes to the wrong location.

The above trial made me think there might be special georeference method or expression which can be loaded by ESRI system but not by QGIS/GDAL, though at the moment I have no idea what it is. Any clue?


texts:
I have this WKT file:

PROJCS["PRS92 / Philippines zone 3"
,GEOGCS["PRS92"
,DATUM["D_Phillipine_Reference_System_1992"
,SPHEROID["Clarke_1866",6378206.4,294.9786982138982]]
,PRIMEM["Greenwich",0]
,UNIT["Degree",0.017453292519943295]]
,PROJECTION["Transverse_Mercator"]
,PARAMETER["latitude_of_origin",0]
,PARAMETER["central_meridian",121]
,PARAMETER["scale_factor",0.99995]
,PARAMETER["false_easting",500000]
,PARAMETER["false_northing",0]
,UNIT["Meter",1]]

How can I use it in mapinfo .prj file so that i can use the CRS?
I'm very new in this kind of matter.


texts: Cannot merge several input datasets in ModelBuilder?

texts:
I am the new one in ModelBuilder, so have faced some troubles with merging my resulted feature class with my set of input dataset.

Could you please help me to find the mistake?


texts: ST_ClosestPoint(line, point) does not intersect line

texts:
In my PostGIS database (PostGIS 1.5 on Postgres 8.4.1), I have two tables: roads (composed of linestrings) and crashes (composed of points). I've been trying to correlate each crash to a road, but I'm having problems getting the following to work:

SELECT ST_ClosestPoint(road.the_geom, crash.the_geom),
    ST_Intersects(ST_ClosestPoint(road.the_geom, crash.the_geom), road.the_geom)
    ST_Distance(ST_ClosestPoint(road.the_geom, crash.the_geom), crash.the_geom) AS distance
    FROM
        --Table crashes is already in SRID 4326
        (SELECT the_geom FROM crashes WHERE gid = 360) as crash,
        (SELECT ST_SetSrid(the_geom, 4326) as the_geom from roads) as road
    ORDER BY distance;

This queury should return the closest point to the crash with gid 360 on each road, but the ST_Intersects function returns false for the first result (the true closest point over all roads). Am I doing something wrong? Is there another way to link the crash to the closest road?


texts: Why the values are NULL if I join the CSV table to a point shapefile

texts:
So my problem is all values are NULL if I join my CSV table to a point shapefile, why? (QGIS)

I tried so many things but I cannot find and handle the problem. I have customer names in the CSV, account numbers as well, digits and letters.

OK, I put a screenshot here, that is maybe better to see what I mean.


texts: QGIS problems executing 'Filter' inside and outside of attribute table with ESRI FileGeodatabase

texts:
I'm working with a file geodatabase in QGIS and am having some issues with the 'Filter'.  These issues occur whether I'm applying the filter to the attribute table or directly to the layer (from the layers pane).

When I go to apply the filter and start writing my expression, all of my fields show up correctly and I can select the one I need, but when I try to view 'All' or 'Sample' of the available values, I'm receiving values from the first field that exists in the feature class.  So, for example, I want to filter on street name, but when I try to view all the values from the street name field, I believe I'm getting the values from the first field that exists in my feature class which is my RelateID field.

I have also noticed when say I'm trying to filter on say a street prefix (I have 7 potential values in my feature class), that is the same number of values that appear in the 'value' window on the side of the query box.  I haven't checked to see if they match up, but i'm guessing the 'value' I'm receiving is the first occurance of that 'Prefix' in my feature class.  Interestingly, when I manually enter my filter expression, "StreetName" = 'MINNESOTA' without getting my value from the value window, the expression executes correctly (selects the appropriate 129 records).

I'm curious if anyone else has encountered issues similar to this before and if so, how did you handle it?

When I save the feature class as a shapefile, the expression builder works fine and all of the correct available values appear in the window as they are supposed to.  Unfortunatly I'm going to be working with 100+ different datasets, most of which will be in the File Geodatabase format and I really don't want to have to convert all feature classes I need to a shapefile and then after processing, etc have to convert them back to File GeoDatabase Feature classes.


texts:
I was just wondering what's the best way to organise QGIS projects. My problem is that I work across different platforms, Mac, Win and Linux and the project contents are stored with links to the files.

I have organised my files in a DATA partition, which gets mounted under /Volumes on Mac and /media etc. on Linux.
Hence, a QGIS project created on a Mac will not be able to load files when you open it on Linux due to the different file paths.

I am not enough of a computer geek to change mount points or anything like that.

I was just wondering if there is a relatively easy way to solve this problem?


texts: Is there a publicly available instance of ArcGIS Server

texts: Are there any publicly available instances of ArcGIS Server for learning? I have searched Amazon EC2 but I don't see any public instances available (unless you already have an ArcGIS Server license).


texts: Add image from my computer for background in Print Composer?

texts:
Add image either by dropdown or icon does not produce any kind of dialog box where I could upload an image. Item properties is checked and showing in the right hand panel, but it's blank (no options showing at all).

Do I need to set properties for this somewhere else (don't see how to do it inside Print Composer)?

QGIS is new to me. It looks to be a powerful suite of tools. Just can't get this working.


texts: Is there any tool/plugin for QGIS to easily select feature and edit attributes

texts: I need to select polygons and change some attributes in a shapefile. I'm using Quantum GIS - I select the polygon, open attribute table and edit the attributes I want. I know there must be an easier and faster way to do this. I was thinking in some tool which when the feature is selectd pop-up a form with the respective attributes and lets you edit them.


texts: Coordinates of the lowest point in polygon

texts:
I am working in QGIS 2.6. I have a file of several polygons and I want for each one, to find the lowest point (where is it located and its elevation).
I am using "Grid statistics for polygons" which gives me the elevation but not the coordinates of this lowest point.
How can I find these coordinates? Is there an existing tool in QGIS to perform this?


texts:
I have a model within a model. Each model has an iterator. The sub model has a iterate feature class instance and the main model has an iterate row selection. I need for each feature class a field calculated based on the individual polygon locations within the area of the selected row.

I have tried to pass the sub model to a Collect Values tool and use that as a parameter for the main model but "make feature layer" does not allow that as an input.

How do I get each instance of the iterate sub model into the make feature layer command to create my selection?

I will elaborate a bit about what is happening in the proxy_background that is causing me trouble before I try the suggestions. The iterate inside of the sub-model is running completely before passing each individual feature class to the main model. Essentially, if my file location has 20 feature classes, it iterates 20 times, not once, before it moves to the main model. This is why I put the Collect Values tool thinking I will at least have the paths to move forward. Alas, it did not want to connect to my main model.


texts: Convert Cassini-Soldner Projection coordinates from Kenya to UTM

texts:
I have a list of coordinates from Kenya. The coordinates were recorded in old survey maps using Cassini-Soldner Projection, with Clarke 1858 Datum. The approximate geographic coordinates (WGS84) of the area where the coordinates are apparently located are lat: -1.46 and lon : 36.9.

As an example, one set of the coordinates is (X: 150052.09 ,  Y: 4690.62). My interest is to convert these set of coordinates to UTM projection.


texts: Selecting Landsat bands to download from Google Earth Engine?

texts:
I am trying to download a specific Landsat scene (since the USGS EROS site is down from the government shutdown) using the imageToCOGeoTiff task to save to my Google Drive. I am using the USGS Landsat 8 Collection 1 Tier 1 TOA Reflectance data.
I get the error message:

  Error: Exported bands must have compatible data types; found inconsistent types: Float32 and UInt16.

I understand that all bands are "float" except for the BQA band which is "int16". I only need bands 1 - 11, not the BQA. Is there any way to select the bands I want (or to remove the BQA band)?

This is what I have so far (I am very new to GEE and most of this code was derived from elsewhere):

var lsat = ee.Image(landsat.filterDate('2016-06-27','2016-06-29').filterBounds(point).first());

// Print the image to the console
print(lsat);

// Define vis parameter
var trueColor = {
  bands: ['B4','B3','B2'],
  min: 6000,
  max: 13000,
};

// Add the image to the map
Map.addLayer(lsat,trueColor,'True Color Image');

// Export a cloud-optimized GeoTIFF.
Export.image.toDrive({
  image: lsat,
  description: 'imageToCOGeoTiffExample',
  scale: 30,
  fileFormat: 'GeoTIFF',
  maxPixels: 11629187501,
  formatOptions: {
    cloudOptimized: true
  }
});


texts: How to create a simple OSM map of San Francisco with a bounding box?

texts:
My objective is to create map of San Francisco, specifically [-122.5145, 37.8092, -122.3716, 37.7085].

Using my (older, 1.8) version of qgis and the open layers plugin, I can find the San Francisco area. I simply wish to export the area bounded by the lat/lon coordinates above to an image file!

Its worth noting that on the openstreetmap site, I'm able to export the bounding box (via the Overpass API only), but the file is nearly 400MB.

Clicking the share icon on openstreetmap does allow me to download an image file.  This is exactly what I am trying to do.  HOWEVER, clicking on the manual selection box only allows the user to draw the area.  I cannot manually specify the bounding box coordinates as I can when exporting the .osm file.

How can I get only the data within the bounding box ?


texts: QGIS: Extracting Raster LAYER Property Information (Source, name, etc.) and property name into a seperate table

texts:
Currently we are trying to extract the layer property information from our vector layers and raster layers (inside some geopackages). More specifically the property information. However though it is very easy to do this with vector layer with Rasters it is a different story. While a work around would be to "vectorise" and do the attribute table part, we are looking for another way.
I know the information in the Raster is there because when I click on Properties on my raster i get this:
[![enter image description here][1]][1]
This is something we also want to run as batch process to process many layers at the same time. Note we are not trying to store the information inside the raster layer but outside of it in a table without geometry.


texts:
Is there a possibility to get the centroid of each line?
When using VECTOR ==> ANALYSIS TOOLS ==> MEAN COORDINATS I always get a centroid of the whole lines. What I need is the centroid of each line without having to select individual lines.


texts:
My company receives GNSS in state plane / IMU in lat/long, but the control is in local. I need to be able to complete AT within Inpho and export to Summit. The Inpho support guy says that I have to convert my GNSS to local, but the Summit guy is certain that's not necessary, although his suggested technique bombs out.

Is there anyone out there who knows how to run the AT from local coordinates only?


texts:
I need to create a Python code for validating domains compared to another set of domains. The initial premise here is taking Collector data, downloading the attributes, and comparing certain domains against a set of parameters. This is based off of an educational geology project where students collect field data and list which geologic regions they're in and identify any feature/rock/mineral. There are only 2 domains I need/can compare to a set list. (also, I can not change how the map has been set up)  I have planned on having the attribute table exported to CSV and comparing that to another, but I need to use Python to automate it (or maybe model builder and extract the code from that.

import arcpy

# Local variables:
Geology2015New = "GeologyofGa2015New\\Geology2015New"
Regions = "Regions\\Regions"
Prov_correct = "C:\\Users\\Joshua\\Desktop\\School Work\\Summer 2015\\Project\\Prov_correct"

# Process: Feature Compare
arcpy.FeatureCompare_management(Geology2015New, Regions, "Province", "ALL", "", "0.003280833333 Meters", "0.001", "0.001", "", "SedimentaryStructure;Fossils;IgneousRox1;SedimentaryRox1;MetamorphicRox1;IgneousRox2;SedimentaryRox2;MetamorphicRox2;GeologicStructure;Strike;Dip;Description;EditDate;Editor;Shape;FID_physio;FID_phys_1;AREA;PERIMETER;PHYSIO_;PHYSIO_ID;FCODE;FENCODE;DIVISION;Shape_Leng;FillTransparency;OutlineTransparency;SHAPE_Length;SHAPE_Area", "NO_CONTINUE_COMPARE", Prov_correct)


texts:
I'm following Mike Bostock's tutorial, let's make a map, but for a map of Mexico. I want to display state boundaries.

I was able to get the map of Mexico from Natural Earth Data's world maps. My source for the state boundaries is here:

http://duopixel.com/LIMITE_ESTATAL.shp

NOTE: this should be the same file as here, but INEGI's server seems to be down for days. I was able to convert it from ogr2ogr into geojson. When you convert from geojson into topojson, this notice comes up:

duopixel$ topojson -s 5 -o estados_simp.json -- estados.json
quantization: maximum error 78.607km (0.707°)
simplification: retained 1130 / 454450 points (0%)

It still seems to generate the file and it looks like it's correct. This is the file it generates. This is my javascript code:

d3.json("estados.json", function(error, mex) {

  var estados = mex.objects.estados.geometries

  svg.selectAll(".estado")
    .data(estados)
    .enter().append("path")
    .attr("class", "estado")
    .attr("d", path)
});

However, I get an error when running this code...

Uncaught TypeError: Cannot read property 'length' of undefined

Now, since I'm not even sure what I'm doing, I have no idea where the failure point is. I'm stuck at this point and any help is appreciated.


texts: How to create a custom coordinate grid in the QGIS 2.14.2?

texts:
I want to create a custom grid different from latitude / longitude, with letters and numbers on the axes.

X = A, B, C, D, D, ...
Y = 1, 2, 3, 4, 5, ...

A  B  C  D  E

1

2

3

A  B  C  D  E

I also intend to serve to use the data in an overview in the print the composer to give the location of the map, "A - 1, or C - 5, etc)


texts: How to quickly fill multiple gaps in a series of adjacent polygons

texts:
I've been given a "background" dataset which is designed to draw as the first layer in a complicated map. For some reason the data provider has previously erased this dataset with roads, parks, etc. The white parts of this image are the areas which have been cut from the  dataset, which is comprised of many individual polygons.

The problem is that the datasets which draw above this layer have changed, leading to many small gaps visible when the datasets are overlaid.

What is the best approach to filling in all the gaps in the background polygon? The full dataset looks like the image below - I want the outer edges to remain, with the white gaps in the middle filled in:

I would normally try Buffer and Dissolve, but I don't have a dataset showing the outside edges of this shape, so I can't clip the resultant polygon back to the original shape. Typically, this dataset doesn't align with other "authoritative" coastline datasets.

I only have a Basic license so don't have access to fancy Topology tools. Am I stuck with manually digitising the outer boundary of this shape?


texts: Converting Multiple Polygons to Raster and ASCII

texts:
I've 27084 polygons representing protected areas. I need to convert to raster (where the protected area will be value 1 and the wider landscape 0), then convert the rasters to ASCII. all of which need the same extent.

Any ideas how to do this other than individually converting each one?

I've access to ArcMap 10.5 and R.


texts: Write to QGIS python console from plugin with sys.stdout

texts:
From my own QGIS plugin, I want to print text to the python console (see motivation below). Most simply, I use

import sys
sys.stdout.write('Outputtext')

And get:

  AttributeError: 'NoneType' object has no attribute 'write'

When I open the python console, close it again and start my plugin, it works. Some kind of initialisation routine seems to be missing when used by the plugin.

The background for this - feel free to suggest alternatives

I'd like to use the console as an output for my plugin unit tests:

unittest.TextTestRunner(verbosity=3, stream=sys.stdout).run(suite)

As the console is inside QGIS, it is non-blocking when running the tests (I need to inspect some results in the map, as I will not implement automated checks for all test results that maybe need to be checked. In that sense, it is rather combined unit and integration testing).

Compared to the MessageLog the console provides a better readability (Syntax highlighting.)


texts: Adding shapefiles to specific .mxd project using ModelBuilder

texts:
I am new to automating geoprocesses.
In ModelBuilder, I conducted a geoprocess of transforming some points in excel to a shapefile, and then to make a buffer of 4 kms to these.
As a result, two shapefiles are produced, like the ones I mentioned above.
I leave you a screenshot :

In the same ModelBuilder, is there a way to automatically add these two shapefiles to a given .mxd project?


texts: Adding graphics from different spatialReference using ArcGIS API for JavaScript?

texts:
I want to know if I'm approaching this incorrectly. We are attempting to create an esri-based map where we add graphics to a map thru the javascript API. The points that we're adding are lat/lng coords(wkid:4326). Part of our requirements are to be able to show these items regardless of the spatialreference of the map. So we tried to use the geometry service that esri provides to project the 4326 coords to the map spatialreference.

That all works fine, unless we have a lot of points to project. If the geometry service is fed too many points at a time it returns Error:"timeout exceeded". We've split out the amount we give the service and called it iteratively, but it also errors out if we call it too many times too quickly. So we've tried to spread the interval between calls to limit that. The more I do this, the more I feel like I'm digging myself into a deep hole.

Should we be approaching this problem in another way? If so, any suggestions?
Should we homebrew our own geometry projection service to bypass the limitations we notice with the esri implementation?

If so, any resources that would help with the required algorithms?


texts: ArcGIS tools for hiking distance calculations

texts:
Objective: Calculate best route (hiking distance/difficulty) from point A to point B using mostly trails and non-trail terrain, but not walking in a big river. B is not necessarily on a trail. In real life: Point A: trailhead in a park. Point: Benchmark or control point monument somewhere in the park. Need to calculate most logical way to get to that point, honoring rules of the park in terms of trying to stay on trails.

What tools do I need to learn to use to be able to perform this analysis in ArcGIS Desktop?

I know about least cost path, but that does not take trails into account. There could be multiple routes using the trail, so the analysis needs to figure out the least distance, and possibly the easiest (some logic describing slope).


texts: How to create a grid whose vertex coordinates are shown on an attribute table in Arcmap?

texts:
I have a specific area to study. Before starting on my field work, I have to create a 310m X 310m grid of the area. Later on I have to get the coordinates of each point of the grid in an attribute table (each angle of each square!). I tried fishnet; but first, it does not cover the area exactly (there are a lot of rectangles outside of the study area); second, when I tried to get the x, y coordinates, it gave an error message (the calculated value is invalid for the row with ObjectId=0...); and third, the output shows the centers of each rectangular polygon, which I do not need.

I am a newcomer to the field, so any help regarding my problem will be deeply appreciated. (I hope I have narrated my problem well.)


texts: Problems using a "custom" projection with Leaflet

texts:
I'm currently trying to display a dataset on a webapp using Leaflet.
All the dataset is projected into French projection Lambert II ESPG 27572.

Instead of converting coordinates, I'm trying to setup Leaflet with this ESPG:27572.
I gave it a first try using proj4leaflet :

var crs = new L.Proj.CRS('ESPG:27572',
    '+proj=lcc +lat_1=46.8 +lat_0=46.8 +lon_0=0 +k_0=0.99987742 +x_0=600000 +y_0=2200000 +a=6378249.2 +b=6356515 +towgs84=-168,-60,320,0,0,0,0 +pm=paris +units=m +no_defs',
    {
        resolutions: [
            8192, 4096, 2048, 1024, 512, 256, 128,
            64, 32, 16, 8, 4, 2, 1, 0.5
        ],
        origin: [0, 0],
        bounds: L.bounds([5168.43, 1730142.53], [1013247.20, 2698564.20])
    }),
    map = new L.Map('map', {
        crs: crs,
    });

L.tileLayer('https://api.geosition.com/tile/osm-bright-3006/{z}/{x}/{y}.png', {
    maxZoom: 14,
    minZoom: 0,
    attribution: 'Map data &copy; OpenStreetMap contributors, Imagery &copy; Kartena'
}).addTo(map);

map.setView([0,0], 2);

It returns a white map with no console message. Commenting the bounds line "did the trick", except the map is centered on South Africa instead of France.

//bounds: L.bounds([5168.43, 1730142.53], [1013247.20, 2698564.20])

I think that I'm not far from the solution, as when I insert a marker which is supposed to be in France :

L.marker([656513, 2626522]).addTo(map);

I'm getting the marker, in the sea near South Africa, but correctly positionned if France was over South Africa.

Expectation :

Result :


texts:
My qgis2web Plugin is not working, as I receive the error as you see below:

Basically, it's able to display layers in Leaflet, but not correctly enough. The problem is similar to these:

QGIS2Web won't display OpenLayers but does display Leaflet
Python error when using qgis2web

As I see on Github there is a lot of issues

https://github.com/tomchadwin/qgis2web/issues

reading all the hints here it looks like the Plugin is no longer supported.
Is there any other counterpart, which I can use?
I would like to populate the GeoPackage layer as an interactive map (i.e. Leaflet or similar) with keeping all the colouration, labels, etc the same.
I tried to convert the GeoPackage layer into GeoJSON, but finally I got a separate layers, which in my geopackage were bounded in one group. When I uploaded it on http://geojson.io I saw everything in the greyscale.
How can I show my geopackage layer interactively?


texts:
Need to import Symbology from ArcGIS PRO to QGIS 3.10 over ubuntu.
Will .SLD do the work?


texts: Need to Set the Map in Layout to fit completely to a layer in QGIS

texts:
So, I know that there is the "Set to map canvas extend" but it just doesn't do the job for me. I used ArcGIS more than QGIS so I am having issues accomodating to the map layout options in QGIS. I want to set the map layout completely to a specific polygon layer that encapsulates the project area. I use QGIS 3.10.

However, when I zoom to layer and "Set to map canvas extend", I get a different outcome:

Now, I am looking to have the lines of the map extend completely fit the polygon extend. How can I accomplish this? So, at the end, I want the grid (blue-white) to be as a border.


texts: Projection of CSV after converting geographic to planar coordinate system

texts:
I have a conceptual question. I have lat-long points which I have used to extract environmental variables from. My coordinates are in WGS1984. For most of the variables, I converted the environmental data into WGS1984 before extracting the data into a matrix/vector (Lat, Long + env variables).

For a few of the variables though, I increased efficiency by converting both my lat-long points and the environmental data into a planar coordinate system. Once I extracted the value for such environmental covariate, I added the values back to this data matrix.

I ran a regression on this data and now that I'm displaying the end result, I'm unsure whether it would be proper to label my data as WGS1984. I was planning on creating a raster out of this XYZ matrix and I have to list a projection in R. Would I have to retransform the extracted data in the planar projected data back in WGS1984 before I can add it into the data matrix? If I understand correctly, these are the environmental values that would correspond to data for the lat-long in either WGS or planar?


texts: 'Exceptions.NameError' > "Name 'sheet_count' is not defined" Error in Python Code

texts:
I know this is something simple but for some reason I am having trouble remembering. In my code I am getting a name is not define for sheet_count. This means I would have to define at top correct. For Example: Is it sheet_count = true? My code is below:

import arcpy, os, sys, string
from os import path as p
from arcpy import mapping as m

#Make parameters for people to choose mxd and folder of PDFs to compare

mxdList = string.split(arcpy.GetParameterAsText(0), ";")
pdf_path = arcpy.GetParameterAsText(1)

#Use Search Cursor to go through Attribuite Table to get Sheet number info

for mapDoc in mxdList:
    arcpy.AddMessage(mapDoc)
    mxd = arcpy.mapping.MapDocument(mapDoc)
    for lyr in m.ListLayers(mxd, "*, Grid"):
        if lyr.description == "*, Grid":
            max_list = []
            rows = arcpy.SearchCursor(lyr)
            for row in rows:
                max_list.append(row.Sheet_ID)
            sheet_count = max(max_list)
            arcpy.AddMessage('Sheet count: %s'%sheet_count)
            print 'Sheet count: %s'%sheet_count
        else:
            arcpy.AddError('No Layers in %s match data source'%mapDoc)


texts:
Is it possible to create custom symbols in mapinfo and use it in thematic shading?How?

  Mapinfo - 10.5


texts:
Following on from How can I add points to a LineString in shapely?, I'm hitting a problem where floating point representations mean objects aren't precisely where I expect them be after they are moved using affinity.translate.

This would be fixed if there was a way of rounding all the coordinates to the nearest cm (or other insignificantly-small value).

Does this exist?

For what it's worth, I'm doing all my intersection and union operations first before the translation and that seems to be working.


texts: Plotting map in R that looks old (yellowish and stained)?

texts:
I want to make a map in R that appear old (yellowish and stained). For example:

Is it possible to achieve this effect entirely within R rather than bringing a map into an image processing program such as Illustrator or Photoshop?


texts: Removing intersect or overlap of points in the same vector layer

texts:
Is there anyway to remove the intercept part of the vector points?

P.S. Please ignore the red circle, it's not a polygon, it is a hand drawn by me to show the error part
So far I've done: "Select by location" and Geometries Checking. But there isn't a way to work.


texts: Saving Bookmarks between sessions using ArcGIS API for JavaScript?

texts:
The ArcGIS Server JavaScript API contains a Bookmark dijit which allows users to add and delete bookmarks. These are not persisted between sessions, but are reset each time the page refreshes.

Has anyone written code which saves the bookmarks between sessions, and which they're willing to share?

Failing that, does anyone have any pointers as to how you'd go about writing this functionality?

Presumably you would need to save each BookmarkItem to a cookie and read them in again when the page loaded - are there any samples I can use as a starting point?

Ideally this would be tied to a user, but more realistically it could be tied to a browser only.


texts: How to create a list and append objects to it through a loop

texts:
Still struggling with this script.  Now I need to create an arpy.AddMessage that says "Projected: with the list of feature classes that were projected in the loop"

i.e.
Projected: filename, filename, filename etc.

I've been at this for 3 days and know I am missing something stupid.  Any help would be appreciated.  Here is my code again:

# Get the spatial reference of the Feature Class from folder A
# Examine all Feature Classes in folder B and report if their
# spatial references match the Feature Class in folder A

import arcpy
arcpy.env.overwriteOutput = True

#try:
#set up the paths
#targetFC = r"C:\MGIS\geog485\Lesson2\Lesson2\CountyLines.shp"
#folderToExamine = r"C:\MGIS\geog485\Lesson2\Lesson2"

targetFC = arcpy.GetParameterAsText(0)
#print "target FC path is "+targetFC
folderToExamine = arcpy.GetParameterAsText(1)
#print "folder to examine"+folderToExamine

#get spatial reference for the target feature class
targetDescribe = arcpy.Describe(targetFC)
#print targetDescribe
targetSR = targetDescribe.SpatialReference
#print targetSR
targetSRName = targetSR.Name
#print "Target feature class is "+targetSRName

# Get a list of my feature classes
arcpy.env.workspace = folderToExamine
listOfFCs = arcpy.ListFeatureClasses()
#print listOfFCs

#Loop through the list of FCs
for currentFC in listOfFCs:
    #print currentFC
    #Read the spatial reference of the current one
    currentFCDescribe = arcpy.Describe(currentFC)
    currentFCSR = currentFCDescribe.SpatialReference
    currentFCSRName = currentFCSR.Name
    #print currentFCSRName

    if currentFCSRName != targetSRName:
        print "Spatial references don't match"
    else:
        print "Spatial references do match"
    if currentFCSRName == targetSRName:
        #skip
       continue
    else:
       # Determine the new output feature class path and name
        outCS = currentFC[:-4] +"_projected.shp"
        #Reproject datasets that are different from target file
        arcpy.Project_management(currentFC, outCS, targetSR)

        #Create a message stating which files were projected
        projFCs = outCS[:-14] +".shp"
        #print projFCs
        projFCList = arcpy.ListFeatureClasses(projFCs)
        x = ", ".join(projFCList) + ", "
        arcpy.AddMessage("Projected " + str(x))


texts:
I am trying to build a REST API using Flask and Postgis.

I have a PostGIS table with a GEOM column. Querying the table using ST_AsGeoJSON gives me a MULTIPOLYGON with lat, lon in WGS84.

[('Reykjavíkurborg', '{"type":"MultiPolygon","coordinates":[[[[-21.9793904176809,64.1449965355369],[-21.9793723947542,64.1449939711492],[-21.9793486086784,64.1449948273916],[-21.979313323496,64.1450030433758], ...

Using the TopoJSON Python port I can manage to render the output in D3, which in turn creates an SVG like so:

What I am looking to do is to export a PNG of that polygon. Effectively rendering the image (rasterizing?) but doing that server side rather than client side.

I've come across what seems to be various different ways of doing this.

How to rasterize simple geometry from GeoJSON file using GDAL?

This method uses the GDAL library. I can load the GEOM column directly into
ogr.Open() but I just get a black image outputted.
Convert geojson to png
Original question is asked in Java but at the bottom there's a Python solution. The Python Pillow imaging library is manually used to draw an image but I'm not seeing any transformation of the coordinate system. This seems to be a rather manual way of drawing the image.
Convert ST_Union polygon to PNG file in PostGIS? Using the ST_AsPNG PostGIS function seems intuitively to be the most practical way but I'm having a bit of trouble building that function in Flask. To output the GeoJSON I do this:

for r in db_session.query(geom.ST_AsGeoJSON()).all():
    ...

but conversely:

for r in db_session.query(geom.ST_AsRaster()).all():

gives me an error: No function matches the given name and argument types. You might need to add explicit type casts. (Googling ST_AsRaster + the error result indicates I've got some catching up to do)
Rasterizing shapefiles with GDAL and Python? Here I find unclear the purpose of the ndsm value as it seems to be neither the input nor the output.

So there seems to be more than one way going forward. And my question is: What is the most efficient way of extracting a PostGIS geom and rendering it as a PNG?

Using either builtin PostGIS functions or Python GDAL bindings seems to make the most sense.

Lots of things I don't really understand such as what bands are in raster files.


texts: Merge vector layers problem with german names (ä, ü, ö...)

texts:
When I'm trying to merge (saga merge layers tool) some shape files I get an error log

  'ascii' codec can't decode byte 0xfc in position 91: ordinal not in
  range(128) See log for more details

Log:

  Uncaught error while executing algorithm
  Traceback (most recent call last):
  Traceback (most recent call last):
    File "C:/Users/user/.qgis2/python/plugins\processing\core\GeoAlgorithm.py", line 232, in execute
      self.processAlgorithm(progress)
    File "C:/Users/user/.qgis2/python/plugins\processing\algs\saga\SagaAlgorithm212.py", line 152, in processAlgorithm
      filename = dataobjects.exportVectorLayer(layer)
    File "C:/Users/user/.qgis2/python/plugins\processing\tools\dataobjects.py", line 302, in exportVectorLayer
      unicode(layer.source()).decode('ascii')
  UnicodeDecodeError: 'ascii' codec can't decode byte 0xe4 in position 86: ordinal not in range(128)

I guess its caused by the names of the vector files (Umlaut like ä or ü). In SAGA it works but not in QGis. Is that a known issue?


texts: Symbolizing multiple occurrences at the same address

texts:
I would like to have the size of the point on my map change to represent multiple occurrences on an event at the same address.  I imported the data from a CSV file.

I'm using QGIS 2.16


texts: errors after installing 64 bit background geoprocessing on ArcGIS desktop 10.2.2

texts:
I installed 64 bit background geoprocessing on my ArcGIS 10.2.2 and then changed my pythonpath to "C:\Python27\ArcGISx6410.2" so that my code can benefit from the 64bit. then I had to install a 64b Scipy, this lead into an error as follows:

RuntimeError: module compiled against API version 9 but this version of numpy is 7
Traceback (most recent call last):

File "C:\Users\Ehsan.abdolmajidi\workspace\ArcPy_matching\NVDB_Matching_v18_H_4\Runner_v18_H.py", line 6, in <module>
from roundabout import *
    File "C:\Users\Ehsan.abdolmajidi\workspace\ArcPy_matching\NVDB_Matching_v18_H_4\roundabout.py", line 8, in <module>
from scipy import spatial
    File "C:\Python27\ArcGISx6410.2\Lib\site-packages\scipy\spatial\__init__.py", line 89, in <module>
from .kdtree import *
    File "C:\Python27\ArcGISx6410.2\Lib\site-packages\scipy\spatial\kdtree.py", line 8, in <module>
import scipy.sparse
    File "C:\Python27\ArcGISx6410.2\Lib\site-packages\scipy\sparse\__init__.py", line 206, in <module>
from .csr import *
    File "C:\Python27\ArcGISx6410.2\Lib\site-packages\scipy\sparse\csr.py", line 13, in <module>
from ._sparsetools import csr_tocsc, csr_tobsr, csr_count_blocks, \
ImportError: numpy.core.multiarray failed to import"

now I am trying to install a new 64 numpy but I can't find a .exe file. .whl file needs pip which when I try to install return new error: "no module named _socket "

do you know if I can find a version 9 of 64 bit Numpy which is .exe? or do you know how I can fix the first problem?

Answer:
I solved the problem by updating the numpy. the problem was that I installed the SciPy which was not compatible with the numpy embedded in the 64 background geoprocessing.


texts: What is some theory on the number of GCP when georeferencing?

texts: Just wondering, is it always better to have more than 3 GCPs? Should the GCPs be spread out evenly? Where can I find some general resources on the theory of this? Other than trial an error, which simply teaches me what works at that moment in time and not what should work in given instances. Not critical, but just interested.


texts: How to Automate Symmetric Difference using OGR in Python

texts:
I need to try and automate the process of creating mask layers by using the symmetric difference of two different polygon layer files (they are currently in .TAB format, however these can be converted to .shp if need be) using OGR in Python.

The following is the code:

outputFileName = new+"mask.TAB"
driver = ogr.GetDriverByName("MapInfo File")
f1 = driver.Open(new+boundname,0)
layer1 = f1.GetLayer(0)
feature1 = layer1.GetNextFeature()

if f1 is None:
  print "Could not open file ", f1
  sys.exit(1)

f2 = driver.Open(new+regionname,0)
layer2 = f2.GetLayer(0)

if f2 is None:
print "Could not open file ", f2

### Create output file ###
if os.path.exists(outputFileName):
  os.remove(outputFileName)
try:
  output = driver.CreateDataSource(outputFileName)
except:
  print 'Could not create output datasource ', outputFileName
  sys.exit(1)

newLayer = output.CreateLayer('SymmetricDifference',geom_type=ogr.wkbPolygon,srs=layer1.GetSpatialRef())

if newLayer is None:
  print "Could not create output layer"
  sys.exit(1)

newLayerDef = newLayer.GetLayerDefn()

featureID = 0

while feature1:

  layer2.ResetReading()
  geom1 = feature1.GetGeometryRef()
  feature2 = layer2.GetNextFeature()

  while feature2:

    geom2 = feature2.GetGeometryRef()

    if geom1.Overlaps(geom2) == 1:
        newgeom = geom1.SymDifference(geom1,geom2)
        newFeature = ogr.Feature(newLayerDef)
        newFeature.SetGeometry(newgeom)
        newFeature.SetFID(featureID)
        newLayer.CreateFeature(newFeature)
        featureID += 1
        newFeature.Destroy()

    else:
        newFeature1 = ogr.Feature(newLayerDef)
        newFeature1.SetGeometry(geom1)
        newFeature1.SetFID(featureID)
        newLayer.CreateFeature(newFeature1)

        featureID += 1
        newFeature2 = ogr.Feature(newLayerDef)
        newFeature2.SetGeometry(geom2)
        newFeature2.SetFID(featureID)
        newLayer.CreateFeature(newFeature2)
        featureID += 1

        newFeature1.Destroy()
        newFeature2.Destroy()

    feature2.Destroy()
    feature2 = layer2.GetNextFeature()

feature1.Destroy()
feature1 = layer1.GetNextFeature()

f1.Destroy()
f2.Destroy()

From this, I get the following output:

However, what I actually want this script to do is create a hole in the larger polygon as shown below

Any ideas where I am going wrong in the code?


texts: Saving to database on OpenLayers map changes?

texts:
I have put a map on website using OpenLayers. This map is fetched from geoserver by using wms request. and in the same way, it is drawn some shapes such as point, line etc. using open layers library in this simple project.

My Spatial data(tables) is held on Oracle spatial database. as I mentioned above, these tables are served by geoserver to the web page and are mapped using open layers library.

How can I save new changes on map to related to Database Table When I draw a new shape on this map?

When it has occurred changes, the second image was created.


texts: Is it possible to catch ServiceException in OpenLayers when ol.Protocol.WFS used?

texts:
I have openlayers-2. Can i intercept ServiceException when openlayers make request to GeoServer? I have found this answer, but for HTTP protocol. However in my code is used WFS protocol. Here is how WFS protocol is created:

function getWFSProtocol() {
  return new ol.Protocol.WFS({
    geometryName: "Obj",
    url: clientSettings.map.layers['wfs'].url,
    featureType: "USA",
    featureNS: "http://www.mapinfo.com/wfs"
  });
}

I have tried to use createCallback like this:

function getWFSProtocol() {
  var protocol = new ol.Protocol.WFS({
    geometryName: "Obj",
    url: clientSettings.map.layers['wfs'].url,
    featureType: "USA",
    featureNS: "http://www.mapinfo.com/wfs"
  });
  protocol.createCallback(function(method, response, options) {
    console.log('wfs protocol callback called');
  });
  return protocol;
}

But callback not called response with ServiceException.


texts:
When I am adding a point to my web map I want to give it a letter that represents the area where the point is, so I need to know what area the point is in.
In OpenLayers 2 I did it like this:
areas[i].geometry.containsPoint(point)   ---> returns true if it does.

Is there anything like this in OL3?

I searched in the OL3 API but didn't find it.


texts: Using related records added on Collector (for ArcGIS) to update original feature class values and symbology

texts:
I currently use Collector for ArcGIS to carry out tree inspections. The design consists of two parts - a feature class layer which contains a point for each tree including data such as its species, height and last inspection date and a table which contains previous tree inspections. The two parts are related using a 1:many relationship class based on a shared GlobalID field. I used ArcGIS 10.6.1 desktop to create to the relationship class which I then host on ArcGIS Online and subsequently update on Collector for ArcGIS (Based on https://www.esri.com/arcgis-blog/products/field-mobility/field-mobility/related-tables-exploring-new-ways-to-use-collector-for-arcgis/).

When carrying out tree inspections, I use Collector to add a related record to trees within the original feature class. Is there a way to use the data recorded in the related record to update attributes in the feature class such as 'next inspection date' without having to change this manually?

For example, a tree's last inspection date was 20/09/2018. If I carry out an inspection today and use Collector to add a related record, is there a way to get the feature class to update its 'last inspection date' to 22/09/2019 without having to do this manually for every tree?


texts: Using Layer Actions on overlapping polygons with PyQGIS

texts:
I have a file with many overlapping polygons that outlines different datasets. In the attribute table I have an absolute Windows path to mapinfo tab files. I would like to define an action on the layer that uses this field to present a list of files and when pressed the file gets loaded into QGIS 3.8.

How is this possible?

I have gotten as far as getting a popup but I get an error when I run the action and also it does not give me the names of the layer in the popup, only the rowid number (I think).

I have this Python script running fine in the console and it executes fine:

location = r"Z:\02 Data Library\Geology\AUSTRALIA\QUEENSLAND\GA 2.5M\GEOL_QLD_Geol_2500k_Ln.TAB"
iface.addVectorLayer(location,"","ogr")

But when I reformat it to use the path column in the attribute table for the layer and actions, I have inserted this code and get the below error:

location = r[%Path%]
iface.addVectorLayer(location,"","ogr")

Error:

An error occurred during execution of following code:
location = rZ:\02 Data Library\Geology\AUSTRALIA\QUEENSLAND\GA 2.5M\GEOL_QLD_Geol_2500k_Ln.TAB
iface.addVectorLayer(location,"","ogr")

traceback.print_exception() failed

Python version:

QGIS version:
3.8.3-Zanzibar 'Zanzibar', 685d8b15d2

Python path:

What am I doing wrong here?

I do not know how to show the "FileName" field when I get the list after clicking on the overlapping polygons in the layer (figure 2)


texts: View multiple tiles distributed over a larger area

texts:
I just started using GEE API to work with bigearthnet dataset. I want to view by highlighting as best as possible (maybe using a marker per image) a small list of images in map, without knowing a specific geographic area. I wrote a snippet below, but I can not determine visually the location of those images.
var tilesToSearch =
           ["S2B_MSIL2A_20180225T105019_59_8",
            "S2B_MSIL2A_20180204T94161_1_17",
            "S2B_MSIL2A_20170831T095029_28_80",
            "S2B_MSIL2A_20170924T93021_25_44",
            "S2B_MSIL2A_20180522T093029_2_75",
            "S2B_MSIL2A_20180525T94030_49_35",
            "S2A_MSIL2A_20170905T95032_27_10",
            "S2B_MSIL2A_20170924T093019_35_34",
            "S2B_MSIL2A_20180515T094029_35_82",
            "S2A_MSIL2A_20171101T94132_7_65",
            "S2A_MSIL2A_20171101T094131_32_59"
            ];

var ic = ee.ImageCollection('TUBerlin/BigEarthNet/v1')
            .filter(ee.Filter.inList('system:index',tilesToSearch))

var count = ic.size();
print('Count: ', count);

var visParams = {bands: ['B4', 'B3', 'B2'], min: 0, max: 3000};

Map.addLayer(ic, visParams);
Map.centerObject(ic, 11);

How can I visualize in the best way those images?


texts:
I have created cell_tower icon in .emf format and inserted in Arc Map 10.1. I want to apply this icon as symbol to point layer.

I am trying but I guess I am doing something wrong.

Tried following things:

1.Going though following links/articles :

Creating new styles , Adding custom north arrow

2.Created new style (with the help of customize tool) but unable to add icon as EMF file.

3.Searching on Google.

Any help will be great !!! Thanks in advance : )

Software : ArcGIS for Desktop 10.1

Screenshot :


texts:
I have a div.

layers

Content Here

 <

Inside this Div I want my LayerControl widget to be placed.

layerControl: {
            include: true,
            id: 'layerControl',
            type: 'domNode',
            srcNodeRef: 'layerTOC',
            path: 'gis/dijit/LayerControl',
            title: 'Layers',
            open: true,
            position: 0,
            options: {
            map: true,
            layerControlLayerInfos: true,
            separated: true,
            vectorReorder: true,
            overlayReorder: true,
                               layerInfos: [
                                      {
                            layer: cities,
                            type: 'dynamic',
                            title: 'i18n.viewer.operationalLayers.cities',
            controlOptions: {
            noZoom:false,
            swipe:false
 }
 }
            ]}
            },


texts: QGIS 3.10 buffer around polygon in kilometer units

texts:
I am trying to make the buffer in kilometers around my polygon.

When I run the "Buffer" option, then I am getting the degrees:

There is quite a few solutions, even similar to mine:

Create buffer in meters?

Understanding QGIS buffer tool units?

How to Change the buffer distance from degrees to meters in QGIS 2.14 Essen?

of which the common denominator is changing the CRS system from WGS to NAD 83.
I done it for my layer, which eventually disappeared from the map.#

In terms of this situation, I started looking for the solution, and I found this:

Shapefile disappeared after changing CRS

However, changing the CRS at the project stage didn't resolve my problem:

So even with the best will in the world, having al of the solutions above I cannot make it through.
Is there a way to get over this problem and set the buffer in kilometers instead of degrees?


texts:
I published simple polygon (shp) as a WMS and enabled caching. In geowebcache I created new task for caching as a test (zoom levels 0 and 1, EPSG 900913, type of operation: "generate missing tiles") and submitted. The task was submitted but no cache was generated! There was no specifed folder with tiles in the gwc directory.

What could be the reason for the absence of tiles? In other words, why generation of tiles doesn't work?

The log:
2017-01-17 16:29:29,782 INFO [seed.SeedTask] - GWC Seeder Thread-11 begins seeding layer : test:mende
2017-01-17 16:29:29,784 INFO [geoserver.wms] -
Request: getServiceInfo
2017-01-17 16:29:29,784 ERROR [geoserver.ows] -
org.geoserver.platform.ServiceException: org.geoserver.platform.ServiceException: Service WMS is disabled
    at org.geoserver.ows.DisabledServiceCheck.serviceDispatched(DisabledServiceCheck.java:69)
    at org.geoserver.ows.Dispatcher.fireServiceDispatchedCallback(Dispatcher.java:576)
    at org.geoserver.ows.Dispatcher.service(Dispatcher.java:571)
    at org.geoserver.ows.Dispatcher.handleRequestInternal(Dispatcher.java:257)
    at org.springframework.web.servlet.mvc.AbstractController.handleRequest(AbstractController.java:153)
    at org.geoserver.gwc.GWC.dispatchOwsRequest(GWC.java:1207)
    at org.geoserver.gwc.layer.GeoServerTileLayer.dispatchGetMap(GeoServerTileLayer.java:618)
    at org.geoserver.gwc.layer.GeoServerTileLayer.getMetatilingReponse(GeoServerTileLayer.java:563)
    at org.geoserver.gwc.layer.GeoServerTileLayer.seedTile(GeoServerTileLayer.java:757)
    at org.geowebcache.seed.SeedTask.doActionInternal(SeedTask.java:136)
    at org.geowebcache.seed.GWCTask.doAction(GWCTask.java:76)
    at org.geowebcache.seed.MTSeeder.call(MTSeeder.java:36)
    at org.geowebcache.seed.MTSeeder.call(MTSeeder.java:25)
    at java.util.concurrent.FutureTask.run(FutureTask.java:262)


texts: Speed up ArcPy CopyRows_management process on large feature layer?

texts:
I am attempting to automate an ArcMap process with Python and I am having difficulty with a particular part of a process taking too long when performed on larger datasets.

The process is as follows:

Step 1: Begin with 4 Feature Layers of geospatial data.

Step 2: Use Spatial Joins to add together for comparison.

arcpy.SpatialJoin_analysis(input_table_1, input_table_2, output_table)

Step 3: Change these into temporary feature layers and use AddJoin_management to join them together

arcpy.MakeFeatureLayer_management(output_table, temp_table)

Step 4: Use AddJoin_management to add these new feature layers together to make a complete master table with all of the previous spatially joined tables.

arcpy.AddJoin_management(temp_table_1, "field", temp_table_2, "field")

Now the further steps require summary statistics performed on this master table but joined tables may not be used as input for summary stats.

I attempted to work around this using the CopyRows_management function and performing summary statistics on the identical table output:

arcpy.CopyRows_management(joined_table, copy_rows_output_table)

However, while this works fine on a smaller table of approx 100k rows, it is now needed to be performed on a feature layer with over 5m entries and it takes virtually the entire night.

Are there any known ways to either:

Convert temporary feature layers to a permanent table faster, or;
Perform any form of editing on a joined table?

I have already attempted:

Running from the command line
Writing master joined table to csv
Using compression on gdb with Compact_management


texts: Patching a large number of Landsat 8 scenes together

texts:
I need to extract land cover information (vegetation, water bodies, etc.) using Landsat 8 scenes.  I have to extract the information on a country scale which means that I have about 15 scenes (different Rows and Paths and different days of the year --> 2015). I am using the Semi-Automatic Classification plugin in QGIS in order to download and pre-process the scenes (conversion from DN to TOA and DOS1 correction). What I thought I should do is pre-process the scenes separately and then merge the scenes for the different bands (i.e. in the end I will end up with one band 1 where I have merged all pre-processed bands 1 of the separate scenes, then one band 2 with all pre-processed bands 2 of the separate scenes, etc....). The image below shows an example of what happens when I merge the first 3 pre-processed scenes of band 2:

You can see clearly the boundaries between the three scenes, but what I find troubling is that when I check the raster values of the same object at the interface between two scenes, the raster values differ significantly. My questions are:

Is it even correct to patch Landsat images in the way I am doing it (since the scenes are from different days of the year). If this is not correct is there a way in which I can accomplish what I am trying to do?
Is it possible to end up with smoother patched up image of the different scenes?
Why are the raster values belonging to the same object in two different scenes so significantly different? Is it due to the fact that images are taken in different days?


texts:
I have recently upgraded to 2.14.11 from 2.8.2. I opened an old project which hasn't been opened since before I upgraded and came across an issue.

I started editing a vector layer which was already present in the project (i.e. not a new layer), created a new polygon and then went into the attributes table to enter the information. I then noticed that there were two rows for the single polygon I had created, and that you couldn't select one without selecting the other.

When I attempted to copy the polygon and then delete the original one causing the issue, one of the original 2 rows was deleted but the second gave the error 'Problem deleting feature'.

Has anyone else had this problem and is it an artefact of moving from 2.8.2 to 2.14.11. I haven't had this problem with any of my other projects from 2.8.2 that I have opened, but I think that I have been creating new layers in them instead of amending old ones.


texts: QGIS2Leaf fails with AttributeError: 'QgsMarkerLineSymbolLayerV2' object has no attribute 'penStyle'

texts:
Can anyone tell me what this means? and what I have to do to get it done?

I'd like to export my map from QGIS to use as a webmap. Data is from a postGIS database.

  Fehler bei der Ausführung von Python-Code:

  Traceback (most recent call last):   File
  "/Users/kevin_hh/.qgis2/python/plugins/qgis2leaf/qgis2leafdialog.py",
  line 297, in export2leaf
      qgis2leaf_exec(self.outFileName, self.basemapName, self.basemapMeta, self.basemapAddress, self.width, self.height,
  self.extent, self.full_screen, self.layer_list, self.visible,
  self.opacity, self.encode2JSON,self.createcluster, self.webpage_name,
  self.webmap_head,self.webmap_subhead,
  self.legend,self.locate,self.address, self.precision, self.labels,
  self.labelshover, self.matchCRS, self.selected)   File
  "/Users/kevin_hh/.qgis2/python/plugins/qgis2leaf/qgis2leaf_exec.py",
  line 575, in qgis2leaf_exec
      opacity: '""" + str(symbol.alpha()) + """', AttributeError: 'QgsMarkerLineSymbolLayerV2' object has no attribute 'penStyle'

  Python-Version:
  2.7.10 (default, Jul 14 2015, 19:46:27)  [GCC 4.2.1 Compatible Apple LLVM 6.0 (clang-600.0.39)]

  QGIS-Version:
  2.8.3-Wien Wien,

  Python-Pfad:
  ['/Applications/QGIS_Wien.app/Contents/MacOS/../Resources/python/plugins/processing',
  '/Applications/QGIS_Wien.app/Contents/MacOS/../Resources/python',
  u'/Users/kevin_hh/.qgis2/python',
  u'/Users/kevin_hh/.qgis2/python/plugins',
  '/Applications/QGIS_Wien.app/Contents/MacOS/../Resources/python/plugins',
  '/Library/Frameworks/SQLite3.framework/Versions/C/Python/2.7',
  '/Library/Frameworks/GEOS.framework/Versions/3/Python/2.7/site-packages',
  '/Library/Python/2.7/site-packages/numpy-override',
  '/Library/Python/2.7/site-packages/matplotlib-override',
  '/Library/Frameworks/GDAL.framework/Versions/1.11/Python/2.7/site-packages',
  '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python27.zip',
  '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7',
  '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/plat-darwin',
  '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/plat-mac',
  '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/plat-mac/lib-scriptpackages',
  '/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python',
  '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/lib-tk',
  '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/lib-old',
  '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/lib-dynload',
  '/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python/PyObjC',
  '/Library/Python/2.7/site-packages',
  u'/Users/kevin_hh/.qgis2//python',
  '/Users/kevin_hh/.qgis2/python/plugins',
  '/Applications/QGIS_Wien.app/Contents/Resources/python/plugins/fTools/tools',
  '/Users/kevin_hh/Desktop/Projekt slowUp Basel-Dreiland']


texts: How to perform (sync or async) network request in PyQGIS

texts:
Following a question on the mailing list about what is the best way to perform network requests in PyQGIS both blocking and non-blocking and with or without support for authentication (QGIS' authcfg)

see https://lists.osgeo.org/pipermail/qgis-user/2019-November/044363.html


texts: How can I create polygons that are on the edge of point data and angle the polygons based on numerical angle data generated from flow direction?

texts:
ArcGIS 3D Analyst 10.5.0.

I have a set of well locations and I would like to define the groundwater capture zone based on the groundwater velocity and the direction of flow. I calculated flow direction using Flow Direction in Hydrology, under Spatial Analyst. Now I want to use this direction data to "draw" polygons at the edge of my well points in the opposite direction of flow.

How can I create a set of rectangle polygons that have an edge centered on point locations and are angled according to a numerical angle I have added to the point data?


texts: GeoPy API error: GeocoderInsufficientPrivileges

texts:
Continually receiving error from standard address geolocation request using Nominatim.

Code below:

from geopy.geocoders import Nominatim
address = "90, Park Avenue, New York City, New York, 10016"
geolocator = Nominatim(address, timeout=10, exactly_one=False)
location = geolocator.geocode(address, timeout=10, exactly_one=False)

Resulting error below:

Traceback (most recent call last):
    File "<input>", line 1, in <module>
    File ".../lib/python2.7/site-packages/geopy/geocoders/osm.py", line 176, in geocode
    self._call_geocoder(url, timeout=timeout), exactly_one
    File ".../lib/python2.7/site-packages/geopy/geocoders/base.py", line 147, in call_geocoder
    raise ERROR_CODE_MAP[code](message)
    GeocoderInsufficientPrivileges: HTTP Error 403: Forbidden

Question: The docs only acknowledge this occurs but does not prevent it (https://geopy.readthedocs.org/en/1.9.1/#geopy.exc.GeocoderInsufficientPrivileges). They wrote,

  "The remote geocoding service refused to fulfill a request using the
  account credentials given."

What can I do to avoid this?


texts: Bulk printing\exporting from ArcGIS for Desktop to PDF or JPG for 300 layers using same extent?

texts:
I am using ArcGIS 10.2.2 for Desktop and I need to print\export individual maps of some 300 point Shape\layers (plant occurrences) to PDF (or JPG). Basemap extent remains fixed. Naming each Map with the Shape Name (spp. name) and text box of spp. name.

Are there ESRI solutions to this, or will this need a Python script?

At first was thinking Data Driven Pages would help; not sure now (looks like that'll do multiple regions of the same layout).


texts:
If anyone is experiencing issues with ArcGIS Online, it may be down.

How can this be confirmed?


texts: What causes 'syntax error, unexpect $end' on iface.showAttributeTable?

texts:
I have a plugin that loads some tables on another dockwidget.
On QGIS 2.8 that's all ok, but on QGIS 2.14 the tables is loaded as layer correctly on a new dockwidget but when I request to open the attribute table I get the error message above.

The code is:

class MyMenuProvider(QgsLayerTreeViewMenuProvider):

      def __init__(self, view, iface):
          QgsLayerTreeViewMenuProvider.__init__(self)
          self.view = view
          self.iface = iface

      def createContextMenu(self):
          if  not self.view.currentLayer():
              return None
          m = QMenu()
          m.addAction("Abrir tabela", self.abreTab)
          return m

      def abreTab(self):
          self.iface.showAttributeTable(self.view.currentLayer())

Inside the plugin´s main class I have:

        lays = []

    #---[Conecta tabelas]
    for i, t in enumerate(lisBOG_t):
        uri_db.setDataSource('', t,'','','')
        lay = QgsVectorLayer(uri_db.uri(), lisBOG_l[i], 'spatialite')
        lays.append(lay)

    #---[Gera painel BOG]
    self.bogtg = QgsLayerTreeGroup()
    self.model =QgsLayerTreeModel(self.bogtg)
    self.model.setFlag(QgsLayerTreeModel.ShowLegend, False)
    self.model.setFlag(QgsLayerTreeModel.ShowLegendAsTree, False)
    self.model.setFlag(QgsLayerTreeModel.AllowNodeReorder,False)
    self.model.setFlag(QgsLayerTreeModel.AllowNodeChangeVisibility, False)
    self.model.setFlag(QgsLayerTreeModel.AllowLegendChangeState, False)
    self.model.setFlag(QgsLayerTreeModel.AllowNodeRename,False)
    self.bogtv = QgsLayerTreeView()
    self.bogtv.setModel(self.model)

    self.bogDock = QDockWidget( "BOG" )
    self.bogDock.setObjectName( "Tab_bog" )
    self.bogDock.setAllowedAreas( Qt.LeftDockWidgetArea | Qt.RightDockWidgetArea )
    self.bogDock.setWidget( self.bogtv )
    self.bogDock.setContentsMargins ( 9, 9, 9, 9 )
    self.iface.addDockWidget( Qt.LeftDockWidgetArea, self.bogDock )

    #---[Carrega painel]
    for l in lays:
        QgsMapLayerRegistry.instance().addMapLayer(l, False)
        self.no_lay = self.bogtg.addLayer(l)

    #---[T_BOG não editável]
    l.setReadOnly(True)

    #---[Menu da camada]
    provmenu = MyMenuProvider(self.bogtv, self.iface)
    self.bogtv.setMenuProvider(provmenu)


texts:
I work on several layers in arcmap 9.3.

My question is about creating binds to display theses layers directly.
like

ctrl+1 ( display / hide  layer 1 )
ctrl+2 ( display / hide  layer 2 )
ctrl+3 ( display / hide  layer 3 )
ctrl+4 ( display / hide  layer 4 )

I Didn't find it in the keyboard menu.

Is there a way of doing it ?


texts:
I would like to know if there is a satellite map for OpenLayers.
I made the following code:
new ol.layer.Tile({
  // A layer must have a title to appear in the layerswitcher
  title: 'Stamen',
  // Again set this layer as a base layer
  type: 'base',
  visible: false,
  source: new ol.source.Stamen({
      layer: 'watercolor'
  })
}),
new ol.layer.Tile({
  // A layer must have a title to appear in the layerswitcher
  title: 'Satelite',
  // Again set this layer as a base layer
  type: 'base',
  visible: false,
  source: new ol.source.Stamen({
      layer: 'terrain'
  })
}),
new ol.layer.Tile({
  // A layer must have a title to appear in the layerswitcher
  title: 'OSM',
  // Again set this layer as a base layer
  type: 'base',
  visible: true,
  source: new ol.source.OSM()
})

And now I need an openlayers satellite map or something compatible with it.


texts: Call Identify Click Function on Dropdown Change

texts:
I am a web programmer with no GIS experience and I am working with a GIS specialist with no web service experience. I am using the JavaScript API 3.7 but I am not very familiar with the API.

I have a select list containing unique identifiers for line segments within a layer. I can currently zoom to the location of the segment on change of the select but I would also like to highlight the selected segment on the map. I have set up a executeIdentifyTask on click of the map. I am not sure if I should call the click method and pass it the coordinates of the segment or if there is a different identify method that I should be using instead.


texts: Writing advanced label expression for ArcMap?

texts:
I have a label that I need to format using label expression. The label looks like this:

  EVA09C

and I want the outcome to be this:

  EA09C

Can I tell the label expression to skip the V in there?


texts: Making a GeoTIFF using Google Earth Data

texts:
I'm trying to make a GeoTIFF using geotifcp, which I'll then upload to GeoServer.

I've got an image overlay (part of a .kml) in Google Earth, so I tried to use its bounding box to create a world file, which I could then feed to geotifcp.

In GeoServer, I've set the declared coordinate reference system to EPGS:4326, which I understand to match Google Earth. To upload the GeoTIFF to GeoServer, though, I need to also specify the native coordinate reference system of the GeoTIFF, which I should be able to do with geotifcp using a proj4 string.

From here, I got the proj4 string that should match Google Earth:

+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs

But when I use that with geotifcp,

geotifcp -e image.tfw -4 "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs" image.tif geotiff_image.tif

it returns this error:

mode=w
Failure in GTIFSetFromProj4

What am I doing wrong?


texts: How to export a map in lower size (11*17 inch) with all the parameters set up for higher map size but still looking good?

texts: I am working with ESRI ArcMap (10.2) and have set up a map page size for 40*60 inch. I created a pdf and printed, it looks good. Next I printed with lower size 11*17 but the place labels does not appear good on the map. Do I need to go again wit ArcMap and set page size for 11*17 and do all those works that I did for bigger size, create pdf and print it. Or is there anyway that can minimize my work load either at ArcMap or pdf environment?


texts: Issue with r.watershed run by QGIS-Python-processing.runalg

texts:
I would like to extract flow accumulation and flow direction layers from a DEM calling the grass r.watershed module from a QGIS-python script using the processing class.
I obtain just an empty map, while, using the processing GUI I get the map I was expecting.
Here is my code:

from PyQt4.QtCore import *
import processing

# "mydem" is the loaded dem layer
dem=QgsMapLayerRegistry.instance().mapLayersByName("mydem")[0]
processing.runalg("grass:r.watershed", dem, None, None, None, None, 10, None, None, 300, None, 1, None, None, None, 40, "C:\\Users\\flc\\Desktop\\QGIS_Plugins\\layers\\acc.tif", "C:\\Users\\flc\\Desktop\\QGIS_Plugins\\layers\\dir.tif", None, None, None, None, None, None)

The 6th argument is the threshold. It is not needed for the computation of flow accumulation and flow direction, but in my trials with the GUI (the processing GUI, not the GRASS GUI) I found that the algorithm works when it is specified.
However, I tried a lot of combinations of parameters.
Do you have any idea about how could it work?


texts: Treating a Vector layer as a file parameter or making sure a vector layer only handles .json files

texts:
This question is related to this one: https://gis.stackexchange.com/review/suggested-edits/192945.
In that Question I had trouble. Exporting my model as a Python script because I was using a workaround to make sure a certain file was treated as a "vector layer" by using an expression. This JSON file is a feature collection and this has to be modified and filtered out. Which the algorithm is for.
To kind of "solve this problem" I am again using a vector layer input, and if I use the ... button and browse to it the file works fine... kind of I noticed there is not output where there should be.
However to make sure no mistakes are made, I have to clarify that the Vector layer input parameter only as a GEOJSON newline Delimited JSON. How do I do this in this line?
self.addParameter(QgsProcessingParameterVectorLayer('Vectorlaag', 'PMKL', types=[QgsProcessing.TypeVectorAnyGeometry], defaultValue=None))


texts: Using LIKE in a For Loop with UpdateCursor

texts:
I have a classified raster with a class value that is a very long concatenated string name (Class_Name below). I'm trying to use arcpy to generate integer labels ranging from 1-24 (which are the third concatenated sub-group in these long strings). Here's an example of a previous (manually manipulated) similar raster attribute table.

So far I've been manually using search criteria in ERDAS to find classes belonging to the same group such as:

$"Class_Names" contains "-14.01-" or $"Class_Names" contains "-14.02-" or $"Class_Names" contains "-14.03-"

I've tried using the following code (as a trial for the first two classes) with LIKE in a for loop with UpdateCursor, but I get a syntax error for my LIKE statement.

in_table = 'l12wm_agsup02.img'
field_names = ['Class_Name', 'crop']
with arcpy.da.UpdateCursor(in_table, field_names) as cursor:
    for row in cursor:
        if row[0] 'LIKE \'1.%\'':
            row[1] = 1
        elif row[0] 'LIKE \'2.%\'':
            row[1] = 2
        cursor.updateRow(row)

I've looked at Using LIKE clause in update cursor gives syntax error? and https://community.esri.com/thread/87848 but can't figure out how to use LIKE in the if else statements instead as a where clause in UpdateCursor.

I'm using ArcMap 10.5, Python 2.7


texts: Reading GPS data from iphone and save data into computer?

texts:
I would like to use my iphone for navigation. I want to read GPS data (latitude and longitude) from iphone and save the data into my computer to use it for navigation.

I am using Matlab in my computer.

Any one can give suggestion/ solution?


texts: Saving user-imported shapefile to geodatabase

texts:
Is there any API or service out there that can do this? We have a SQL Server managed geodatabase (ArcSDE).

It would be ideal to be able to create a web application where someone can upload a shapefile and we save it to our geodatabase (we'll do further processing on it as the person works their way through the web app form). Is that possible from any platform? Javascript API? ArcGIS Online? .NET with ArcObjects?

I cannot find any documentation at all for any platform regarding this capability which makes me think it's not possible. I've seen documentation for displaying user-imported shapefiles on a map, but not actually saving them to a geodatabase. We have no need to do any mapping with them.


texts: Sequential numbering from North to South using ArcMap

texts:
Is it possible to label points from North to South sequentially?
I have about 2000 points and the table is not in order. I know how to number sequential in the table, but I need it to be sequential from north to south. (like below)


texts: How to use multiple field calculator in QGIS graphical modeler

texts:
I find myself struggling to understand how to chain field calculator process in the graphical modeler in QGIS.
I have build a model that take 3 layers input and 3 vector field parameters and try to make some calculation with the field calculator, 3 time one after another. I stuck in the first one and the thirdץ
The first take pre-calculte value (choose the field name from the parcels layer in the start of the model), and calculate the concatenation of the numeric fields "GUSH" and "PARCEL", with the "_" sign between them = "1000_10" for example.
concat(to_string(@SelectGushField) || '_' || to_string(@SelectParcelField))

the second calculator calculate the area, as name "Calc_Area" (which not show in the next field calculator)
$area

the third calculate the ratio between the legal area of the parcel and the calculate area of the same feature
  @SelectLegalArea / @add_calc_area_field_OUTPUT

this calculate is also not working

using QGIS 3.16.1 on Windows
Can someone explain to me how to solve this?


texts: make a matrix from rasterLayer / plotting custom coloured raster in ggmap

texts:
I'm trying to add a raster layer to a ggmap object. The basic data is a .tif and as I get strange errors by loading it by the package tiff or rtiff, I decided to make first a rasterLayer (it works) and than add it to the map.
The ggmap() + inset_raster() works for me, however, I cannot assign custom colours (which are predefined) to the raster layer (see this post).
That's why I ended up with ggimage, which requires a matrix, array, imagematrix or something which understands as.raster.
When I convert my RasterLayer to one of these formats and put it in ggimage, I get the error message:

Error in if (max(mat) > 1 || min(mat) < 0) { :
missing value where TRUE/FALSE needed

The .tif is quite big (107Mb), which may be the reason why I cannot use tiff/rtiff (converting to polygons seems not to work for the same reason), but how could I handle the error of ggimage?
Or do you have maybe a better idea to get the raster into the map?
P.S.: It does not have to be ggmap, but I need to use a dynamic change of layer extent later, which is only possible with coord_map without weird displaying (see this post).


texts:
I have a point. I want to add to a Point layer in QGIS 3.16, but I can't get the Z-dimension working on my layer. The point is created but has no z-value...
coord = qgis.utils.iface.activeLayer().crs().authid()
r_layer = QgsVectorLayer("Point?crs=" + coord, "MyPoints", "memory")
r_layer.setProviderEncoding(u'UTF-8')
pr = r_layer.dataProvider()
#CODE that handles fields#

zPoint = QgsPoint(55, 34, -2)
seg.setGeometry(QgsGeometry(zPoint))
pr.addFeatures([seg])
r_layer.updateExtents()
QgsProject.instance().addMapLayer(r_layer)


texts: Counting all features in all layers of ArcGIS map to get grand total?

texts:
I have 20-50 polygon layers in an mxd and I want to find out the aggregated amount of all features (global count). The shapefiles which are the source of those layers are stored in different folders. Just knowing the count would fulfill my requirements because it is just a quick check before I go on processing these files. Thus, I don't want to merge the files or count separately. I just need the simple amount.

Up to now, I can only think of a solution by coding a script in ArcPy, but I would prefer to find an easier way.
I use ArcGIS 10.1.


texts: CDO without proj and hdf libraries (Regrid HDF file AMSR-2)

texts:
I am trying to use CDO for AMSR-2 data analysis, but when I installed CDO 1.9.9, it does not include Proj and HDF library. Even when I installed Proj and HDF libraries in Ubuntu, the CDO still unable regrid HDF file (cdo mygrid, infile outfile).
How can I include Proj and HDF libraries in CDO?
I am using Ubuntu 20.04 and CDO 1.9.9.