texts: Equivalent function to shapely's envelope in ogr?

texts:
I wrote a very simple python function that gives me the envelope/bounding box/...based on two coordinate pairs (ower left and upper right corner):

# ...
def get_bounding_box(self, ll_corner, ur_corner):
    try:
        from shapely.geometry import MultiPoint
        mp_env = MultiPoint([ll_corner, ur_corner]).envelope
        return ogr.CreateGeometryFromWkt( mp_env.to_wkt() )

print mp_env
POLYGON ((4.8439699999999997 52.3661099999999990, 4.8469199999999999 52.3661099999999990, 4.8469199999999999 52.3706000000000031, 4.8439699999999997 52.3706000000000031, 4.8439699999999997 52.3661099999999990))

As you can see I'm returning an ogr geometry which I use, for instance, in this manner:

bbox = geo.get_bounding_box(lowerleft,upper_right)
...
p=ogr.Geometry(ogr.wkbPoint)
p.AddPoint(x,y)
if p.Within(bbox):
    ....

Isn't there an equivalent function to shapely's envelope in ogr? Both .ConvexHull() and GetEnvelope() return only the two corner coordinate pairs.


texts:
I want to know, can we use GeoExt freely? is there any need to purchase some license for using GeoExt ?

Actually i am already working with free Open source tools (OpenLayers+GeoServer) for my webMap applications. But i need some specific functionality which is possible with GeoExt. so can i use GeoExt freely with Openlayers and GeoServer ?


texts: Assigning attributes values to PointZM and convert to LineZM Geometries in QGIS?

texts:
I have elevation and measured values in attribute table of points layer.

Is there a way to convert this layer into PointZM (with elevation and measured) and assign elevation and measured values to PointZM geometry in QGIS/PostGIS?

How do I convert PointZM to lineZM geometry by preserving elevation and measured value of each vertex in lineZM in QGIS/PostGIS?


texts: WFS GetFeature control doesn't take 'viewparams' into account

texts:
I use GeoServer to get my WMS data. Everything works great for showing data in OpenLayers.

GeoServer has the ability to create SQL views and use parameters to return certain geometry. These parameters are called viewparams. This could be an example of a SQL view:

select *
from myTable
where age = %ageParam%

When using viewparams you are able to pass the value for the 'ageParam' parameter. If a parameter is not set, a predefined default value will be used. More information about this topic can be found here.

In OpenLayers however,  "getFeature" control does not take the 'viewparams' parameter into account. This means it always returns features with the default parameter values.

Is this a bug, a missing feature or something I just missed.

Edit 1

Here is a screencast I made to visualise the problem.

http://screencast.com/t/RgrQvaoyoqz

As you can see, the layer is redrawn with new points but the selectFeature control still returns the old features.

Edit 2

<wfs:GetFeature xmlns:wfs="http://www.opengis.net/wfs" service="WFS" version="1.1.0" xsi:schemaLocation="http://www.opengis.net/wfs http://schemas.opengis.net/wfs/1.1.0/wfs.xsd" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">
<wfs:Query typeName="test:point_query" xmlns:test="http://test.boedy.p:8080/wms">
<ogc:Filter xmlns:ogc="http://www.opengis.net/ogc">
<ogc:BBOX>
<ogc:PropertyName>geometry</ogc:PropertyName>
<gml:Envelope xmlns:gml="http://www.opengis.net/gml">
<gml:lowerCorner>490488.82657346 6800062.1058854</gml:lowerCorner>
<gml:upperCorner>491071.65891403 6800556.5579121</gml:upperCorner>
</gml:Envelope>
</ogc:BBOX>
</ogc:Filter>
</wfs:Query>
</wfs:GetFeature>

**edit **

This is the JavaScript code for the layer and controller

test = new OpenLayers.Layer.WMS("test",
    "http://test.boedy.p:8080/wms",
    {
        'layers': 'test:test_query',
        viewparams: 'x:'+myX+';y:'+myY,
        maxfeatures: 2000,
        ratio: 1,
        transparent: true,
        isBaseLayer: false}
    );

selectCtrl = new OpenLayers.Control.GetFeature(
        {
            protocol: OpenLayers.Protocol.WFS.fromWMSLayer(test),
            clickout: true, toggle: false,
            multiple: false, hover: false,
            box: true,
            toggleKey: "shiftKey", // ctrl key removes from selection
            multipleKey: "shiftKey",
            eventListeners: {
                featureselected: onFeatureSelect,
                featureunselected: onFeatureUnselect
            }
        }
    );


texts: Reading metadata from existing datasets in ArcCatalog?

texts: I would like to know which ArcPy functions can be used to read metadata from existing datasets in ArcCatalog?


texts: What is a simple python code to practice with in arcgis?

texts: I am trying to familiarize myself with python code. Could someone give me a simple code with an explanation of what it does so I can practice? Or a few codes with explanations if you have the time


texts: Locate features along routes: measurement unit output?

texts:
I am using ArcMap version 10.3.1:

I have a polyline (also converted to a route feature) and points layer. All points lie on the line. I want to calculate the distance each point is from either the start or end of the polyline.

I have used the Locate Feature Along Routes tool (with metres as the unit), and have not selected options to the 'keep only closest route location' and the 'use M direction offsetting'.

I get an output which makes sense for relative distances points should be from each other, but the units appear to be wrong (possibly in decimal degrees?). I am confused because the units of my map projection are in metres, i.e. when I use the 'measure' tool on the map it gives me back metre measurements.

Is there any other reason I have missed as to why the Locate Features Along Routes tool has not outputted measurements in metres?


texts:
I have a dataset of roads I am displaying in Arcgis Pro. I would like to label the roads, however, obviously there are large numbers of them so they do not display clearly and often overlap. I would like to create a buffer around the labels to reduce the amount being shown, something you could easily do in Arcmap.

I have found the buffer option in ArcGIS Pro however it is now limited to a 50% of the size of the label.

Is there anyway to increase the buffer of labels beyond this size?

Ultimately, I am trying to only display some road names (I'm not fussed on which ones get displayed) in a neat and legible manner for a map.


texts:
I have a large database of point data throughout all of Eastern Australia. Each point represents species presence/absence data and has attached about 18 predictor variables. To account for sampling intensity, I would like to make a 1km x 1km grid and count the frequency of the number of presence records (this would be my response variable) and create another column in the data frame that is the sum of all points within the cell (this would become an offset variable in a GLMM).

What is the best approach to this problem? I know of the Point to Raster tool, however I dont think it allows for keeping the other predictor variables. Aggregate Points is another option, however I am using ArcGIS 10.3 and I dont have access to the ArcGIS GeoAnalytics Server.


texts: Adding columns with easting and northing of points using ArcMap?

texts:
I've got LSOA dataset with centroids weighted by population. What I am trying to do is get coordinates of those points (centroids) so that I can then calculate distances to centroids outside of ArcMap (not the full licence).

Does anyone know how to do that?

Centroids are geocoded and do appear on a map.


texts:
I received “Hillshade” data that I'm trying to use in QGIS.

When I check the "Pseudocolor" option under "Color map", I get many colors, apparently corresponding to the terrain, though it doesn't jump out in an intuitive way.

"Freak out" also does that, w/ a different color scheme. (Don’t know if these options came with QGIS or with my data.) But "Colormap" doesn't show up at all, and Grayscale is just uniformly gray.

I expected grayscale to exhibit shades that would make the terrain stand out intuitively - but I don't really know what to expect.

Any idea what my problem might be?


texts: ArcServer 10.1- Get Original Document location (mxd) from published web service using python

texts:
I need to extract the document location of a published web service. Is there any way this can be accomplished in python?
I know where the maps are sitting on the server, and i know the name of the service and what mxd file they're based on, but i need to accomplish this programatically.

Thanks for your help!!


texts: How to calculate the area of an ellipse on the earth's surface using Python?

texts:
I have a dataset consisting of GPS coordinates. I use a clustering algorithm to generate a two dimensional Gaussian for me. The Gaussian takes the shape of an ellipse if visualized on the map.

However, in order to compute the area covered by the ellipse (the axis could be a slope, not necessarily parallel with latitude or longitude), I need to project this ellipse into a planar surface with equal scale, for ease of calculation. I need the area to be in square kilometers.

Could anyone point me to some keywords (for example, what system I shall project the ellipse to) and libraries I could use in Python?

I just got to know PyProj but am not sure how it could suits my needs (what system to use etc).


texts: Joining tables without losing original field names in ArcGIS for Desktop?

texts:
I'm wondering if there is any method to preserve the original field names in ArcGIS when I join csv. tables to a Shapefile.

I ask, because sometimes I  join tables with 10, 20 or even more columns and when exporting the joined data ArcGIS always renames them according to the name of the original table and the order of the column (e.g. table_xy_1,table_xy_2, etc.).

Since there is no easy way to rename the fields in the attribute table permanently (i.e. you have to create a new column and copy the old on into this new column AND delete the old column) this creates a bunch of work (despite from the fact that you have to orientate yourself in this mess of field names).


texts: How to set Marker Line symbol for QgsVectorLayer by using Python?

texts:
QgsVectorLayer Right click->Properties->Style->Single symbol->click on simple line in tree node ->symbol layer type -> Marker line

Kindly help me to set marker line symbol programmatically by using Python.


texts: QGIS how to set up buffer units to meters?

texts: I want to create buffer and set up distance in meters as units. The file is in WSG 84/UTM zone 18N. The project CRS is also set to WSG 84/UTM zone 18N. Units for distance measurement are set to meters (Project Properties/General settings) and Map tools preferred distance units are sets to meters. OTF is enabled. I use Fixed distance buffer tool and results are wrong it seems like a buffer is created in degrees not in meters. Is it a problem with a projection or the settings are wrong?


texts: Changing label next to symbol in Table of Contents in ArcMap using ArcPy?

texts:
I'm trying to write a script to change the label next to the symbol in the TOC so it appears in the legend. ArcPy has a way to change the name of the layer, but I can't seem to find a way to change the text one level below that. They also have a way to change unique symbol texts and descriptions, but that's not what I need either. Am I just missing it in the tutorials/help pages?

It's just the text box in the layer's Properties under the Symbology tab. There's a box that says "Legend" and then "Label appearing next to the symbol in table of contents."


texts: QGIS clip tool doesn't regenerate geometry column in attribute table

texts:
I have a lot of large datasets which all contain WKB defined geometry (point and polygon separate), which makes it easy to load in QGIS from a PostGIS-database. This format is also preferable when using Geoserver to display data through WMS later on.

When i clip the set of polygons (with a lot of other polygons) the result is good visually, but the new shapefile does not have a geometry column, it rather generates only a point_geometry, while the polygon_geometry disappears. The result is still shown in polygon-format in QGIS.

When trying to load this shapefile into PostGIS through the PostGIS Shapefile import/export manager (With "Create spatial index automatically after load") I don't get the polygon-geometry recreated, and I am not able to use the table either in QGIS or Geoserver at all as the table has "No geometry"

I've searched for ages trying to find a viable solution to regenerate these WKB-type geometries on the clipped polygons, as this is LARGE datasets and many different clips. Anyone have any idea how i can get this done?


texts: Avoiding using same road twice in pgrouting?

texts:
I am currently developing an application to create tours for motorbikes. At the current stage of development there is a start/endpoint and several points which outline the tour. Between these points the pgr_dijkstra function is used (with a customized cost function). The problem is best described by the picture attached.

Is there a way to avoid an edge getting traveled twice? I have tried to use the reverse cost column with a high cost but as the graph is not directed I get odd results. At least the pgrouting documentation states that the graph has to be directed to use the reverse cost if I understood that correctly.

sql := 'SELECT DISTINCT ON (seq)
        seq, id, the_geom, id1 AS node, id2 AS edge, cost, km
            FROM pgr_dijkstra(''
                 SELECT id,
                 source::integer,
                 target::integer,
                 my_cost::double precision AS cost,
                     FROM network
                     WHERE network.the_geom && bbox
            '',' || source || ',' || target || ', false, false),
            network WHERE id2 = network.id';

This is the statement to create a part of the route. It is nested inside a loop which passes new source and target in every iteration as long as there are some. This code is executed and returned record-wise.


texts: Standalone python script - Qgis - QgsVectorLayer points to wrong feature class in gpkg

texts:
GPKG file has multiple featureclasses in it. I create a layer to a specific featureclass using QgsVectorLayer but it always points to the first FC in the gpkg file.

This is a stand alone python script, not the console under the GUI.

output_gpkg = r'W:\srm\wml\Workarea\mamcgirr\Q_universal_overlap_tool\work\qscript_outputs\one_status_common_datasets_debug_version_dqu_shape_2.gpkg |layername = Legal_OGMA'

fc_input_layer = QgsVectorLayer(output_gpkg, 'test', 'ogr')

if not fc_input_layer.isValid():
    print("layer invalid.")
else:
    print("layer valid.")

for field in fc_input_layer.fields():
    print(field.name(), field.typeName())


texts:
I work with point layer in this size in the map:

When i create a Legend, the symbol size is to big:

I try to avoid from solve this issue by:

decrease the symbol size in the layer properties.
convert the Legend to graphics

How can i decrease the symbol size in the Legend and still remain with big symbol size in the map?


texts: Different colors on the same polygon in 3D view of QGIS

texts:
I have a .shp file with buildings represented as polygons.

My goal is that, in a 3D representation, I am able to have different colors in a single polygon (building) depending on its height. For example, from 0 to 10 meters, buildings are yellow, then from 10 to 20 m they are blue, and so on.
Is there any way of doing this in QGIS or using PyQGIS?
I have tried to add copies of the buildings to the shapefile with different heights and the number of copies depends on the height of the building, and then pick a different color for each building copy, but it does not work, still only one color appears.

EDIT: This is an exemplification of what I'm looking for:

EDIT2: Using the Properties > 3D View > Rule-based, and then View > New 3D Map View, I get a 2D map :( :


texts:
I am a total GIS beginner and have been trying to make some simple maps of the nation of São Tomé and Príncipe using QGIS and publicly available data. The readily available data I find comes from OSM.  Some features, such as waterways, coastlines, and cities/towns/villages I have been able to find as individual shape files and add to my QGIS projects as layers that can be turned on and off.  But other features I can only find as part of one "flattened" (computationally-speaking) file that renders like this:

This data was obtained by selecting 'OpenStreetMaps' from the 'XYZ Tiles' option in the QGIS browser pane, but I've found the same data various places on the internet.  There are lots of great features visible, but they all come as a single layer, and I can't find any way to select/de-select/modify them individually. For example, the district boundaries in purple and the park boundary in green.

However, my boss sent me some .png files given to him by some long-ago employee and wants me to create similar maps:

Here we have the same district and park boundaries, but whoever created this map was able to manipulate those features individually.  I've been looking high and low but can't find individual shape files for those features on the internet.

Does anyone know where they might have come from?

I also notice that the roads are broken out here as 'Primary roads' and 'Other roads', whereas the OSM shape file I've found only has one class: 'roads'.

I feel like (even with my rudimentary skills) I can re-create these maps as requested, if only I can find the underlying data!


texts: QGIS 2.10.1 GRASS Processing algorithms not working

texts:
I just did a fresh install of QGIS 2.10.1 from the qgis repos (Ubuntu 14.04). I also have grass70 installed (through ubuntugis-unstable).

I can see the GRASS7 commands in the processing toolbox panel, but when I try to use them I got the Missing dependency. This algorithm cannot be run :-( message.

In the plugins log I noticed that QGIS could not load libgrassplugin6.so, but I'm guessing that is not the issue since I'm trying to run grass7 algorithms.

Also, on the processing log I got this output when trying to run a grass 7 command:

GRASS GIS 7 execution commands
        g.proj -c proj4="+proj=utm +zone=30 +ellps=intl +towgs84=-87,-98,-121,0,0,0,0 +units=m +no_defs"
        v.in.ogr min_area=0.0001 snap=-1 input="/usr/share/qgis/python/plugins/processing/tests/data" layer=points output=tmp1444494282425 --overwrite -o
        g.region -a n=4458983.8488 s=4458921.97814 e=270855.745301 w=270778.60198 res=100
        v.voronoi input=tmp1444494282425 output=output22e863becb734ae4a06720411090abf0 --overwrite
        v.out.ogr -s -e input=output22e863becb734ae4a06720411090abf0 type=auto output="/tmp/processing/b5eb2443ef0840239dc22209a29a71db" format=ESRI_Shapefile olayer=output

But when I try to run a GRASS command (without version, I think that means grass 6, doesn't it?), I got no output on the log. I think that must mean something.

I started QGIS from a console and got this error when I try to execute a GRASS 7 algorithm:

ERROR 4: Unable to open /usr/share/qgis/python/plugins/processing/tests/data/points.shp or /usr/share/qgis/python/plugins/processing/tests/data/points.SHP.

At first I thought it was a permissions issue (the file exists), so i chmod it to 777 but the error persisted.

Browsing around the processing plugin code, it seems that every time you run a grass algorithm the plugin runs a test using that file.


texts: Using directional shading on boundaries for cartographic effect in ArcGIS Desktop?

texts:
The attached screenshot of a National Geographic map shows some beautiful directional shading of boundaries.  I'm trying to reproduce this sort of shading for a burn perimeter map, as I want to highlight features within the burn perimeter and not overshadow them with a transparent polygon or crosshatch.  Is there a way to accomplish this sort of boundary shading using ArcGIS 10?

Link to larger map titled: "1946 Political Subdivisions Of India By National Geographic.jpg"


texts: Replacing Street Directional Abbreviations (N, S, E, W) with spelled out value (North...) using ArcMap ModelBuilder

texts:
In converting street centerlines to Next Gen 911 schema compliance, how do I iterate a !St_PreDir!.replace("E", "East") to include replacing N with North, S with South, and W with West?

I am trying to automate in the ArcMap ModelBuilder and I have no Python experience.


texts:
I have a raster in UTM projection. The proj4 string is: +init=epsg:25833 +proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
I need to reproject the raster to DHDN / 3-degree Gauss-Kruger zone 4 (EPSG: 5678). So I did:
library('raster')

dtm <- read.table('333105624_dgm1.xyz', header = F, dec = '.')
dtm.raster <- rasterFromXYZ(dtm, res = c(1, 1), crs = crs('+init=epsg:25833'))
dtm.raster.rpj <- projectRaster(dtm.raster, crs = crs('+init=epsg:5678'), res = 1)

Unfortunatly, the reprojection doesn't align with the original:

The reprojected raster is shown above the original projected raster and as you can see it doesn't align with either the map or the original projected raster.
What do I miss here?


texts: What is this coordinate system from ESRI and how do I convert it to regular coordinates for GeoJSON Mercator?

texts:
I downloaded the ArcGIS health region boundaries file in ESRI Lambert Conformal Conic format here: https://www150.statcan.gc.ca/n1/pub/82-402-x/2018001/hrbf-flrs-eng.htm

Using mapshaper.org I exported the data to GeoJSON but the coordinates are not a system that I recognize, it looks like:

[-8650999.06360543,5474465.5029140925],[-8650958.754519965,5474511.452623501]

How do I convert these ESRI files to GeoJSON I can use with Leaflet maps?


texts: ArcGIS Python if statements say field not nullable?

texts:
I am creating a new field - let's call it A (SHORT INTEGER) in shp file using Modelbuilder and I need to fill it it with the values from two different fields. First step is to fill field A with the values from field B (that's a string field which is partially empty). The next step is to change all of the values (in A) which are 0, with values from another field with the result of equation C/3. Normally I would select the 0 values in the table and execute normal field calculation, but I'd like to have it as one statement in modelbuilder.
That's what I wrote so far in Python

def function(A,C):
  if (A == 0):
    A == (C/3)

function(!A!,!C!)

I get the error that the field A is not nullable, but while creating it, the option is set to 'nullable'.


texts:
I want to convert a json into topoJson using mapshaper.org.
The json file I am using has world wide names in it, hence is encoded in UTF-8.

However, when I upload the file on mapshaper and then click export - TopoJson, it saves it in an encoding which encodes some of the letters wrong.

Is there a way to tell mapshaper to encode in UTF-8?

Would be interesting to know how to export any file edited with mapshaper in an specific encoding.


texts: How to fix Error 14 -- Unable to open ~/.qgis//qgis.db?

texts:
I get this error:

Error 14 -- Unable to open ~/.qgis//qgis.db -- QGIS version 1.8.0-Lisboa -- Ubuntu 12.04

I had define a custom CRS and that error occurred when I tried to save that CRS under "settings"→"custom CRS" in QGIS 1.80.0-Lisboa on Ubuntu 12.04.

I notice that the qgis.db file is located at:

/usr/share/qgis/resources/qgis.db

I created a link but this did not work either. Can I tell QGIS to look for the qgis database at the above location rather that ~/.qgis//? Looks like permissions are read and write for all users.


texts: Landsat images from same path/row have different shapes

texts:
I'm starting to play with satellite imagery. In particular, I am using images of the Landsat 7 and Python. I open two TIF images of the same place but at different times to the numpy arrays and they have different shapes.

import numpy as np
import rasterio

with rasterio.open(image1_path) as src:
    B1_t1 = src.read().astype(read_astype)[0]
with rasterio.open(image2_path) as src:
    B1_t2 = src.read().astype(read_astype)[0]

B1_t1.shape
B1_t2.shape

#OutPut
#>>>(7811, 6951)
#>>>(7791, 6941)

These two images are from the same path/row.


texts: Creating Cloud free images out of a MOD09A1 MODIS image in GEE?

texts:
How can I mask clouds in a MODIS\MOD09A1 images in GEE?
I know I should do sth with the "Bitmask for QA", but how?

var myImgCln = ee.ImageCollection("MODIS/006/MOD09A1")
    .filterMetadata(???


texts:
I'm working with someone who exclusively uses ArcGIS Desktop. I run QGIS and no longer have access to a copy of ArcGIS Desktop.

Is there a way to develop a plugin in QGIS and then use it in ArcGIS Desktop?


texts:
Before I installed ArcGIS , I have installed python2.7 with many modules and then the ArcGIS installed another python2.7 in other folder.

Now I want to use the arcpy in own old installed python,because I dont want to install these modules again.

I have added Paths in the PYTHONPATH system environment variable included

D:\software\arcgis10.2\python2.7\ArcGIS10.2\Lib
D:\software\arcgis10.2\python2.7\ArcGIS10.2\Lib\site-packages

and add paths in the PATH included

D:\software\arcgis10.2\arcgis10.2\Desktop10.2\arcpy
D:\software\arcgis10.2\arcgis10.2\Desktop10.2\bin

However, when I import arcpy ,I still encountered errors ImportError: No module named arcpy


texts: Using JavaScript variable as date parameter in SQL view

texts:
I'm creating a map based on an SQL view in GeoServer, and I want the layer to change depending on the date parameter in the SQL view (%date%). It works when I set my variable to 'current_date', but if I want to go 1 week further with 'current_date + 7' it doesn't work. I don't get any errors when I inspect the page, the layer just won't appear on my web map.

It's like SQL doesn't recognize 'current_date + 7', but it works fine in pgAdmin. I have also tried it with "current_date + interval '7 days'" and that doesn't work neither.

The code that works:

var datum = "current_date";
var heatmap = new ol.layer.Image({
            source: new ol.source.ImageWMS({
                    url: 'http://localhost:8080/geoserver/wms?',
                    params: {
                            'LAYERS': 'heatmap:heatmap_met_slider',
                            'viewparams': 'date:' + datum,
                            'VERSION':'1.1.1'
                    },
                    serverType: 'geoserver'
            })
    });

The code that doesn't work:

var datum = "current_date + 7";
var heatmap = new ol.layer.Image({
            source: new ol.source.ImageWMS({
                    url: 'http://localhost:8080/geoserver/wms?',
                    params: {
                            'LAYERS': 'heatmap:heatmap_met_slider',
                            'viewparams': 'date:' + datum,
                            'VERSION':'1.1.1'
                    },
                    serverType: 'geoserver'
            })
    });

Any advice?


texts: Merging many map sheets together into one large sheet for later georeferencing?

texts:
I have hundreds of images, each one is a scan of a map sheet, and I want to georeference them. The maps cover a large area (a whole country), they just split the area into many smaller cells. Each map sheet aligns with (but doesn't overlap) another sheet. The sheets are predicably named, so I can "go north" of the current sheet etc. Each scanned sheet has a border (the map image wasn't pushed up against the edges of the scanner). Each sheet also has an internal paper border. Each scanned sheet was not aligned to the scanner, so each sheet will be skewed/tilted a little bit. I know I have to manually crop all of them.

I know how to georeference using QGis etc., but it's just the task of doing ~500 individual georeferencing that's daunting. Is there any way I can, once I have the cropped skewed images, 'merge' the sheets together to produce a small number of larger images?

The top left corner of the cropped image A, should match up with the bottom left corner of image B, etc., so is there a way I can automatically 'glue' image A and B together so I only have to georectify the final image?

I'm not a GIS professional, so apologies if this is a hard to understand question. I'm using Linux, so has to work there.


texts: Distance from point stored in geometry with SRID 4326

texts:
I have stored 2.5 million data points (addresses) using the geometry data type in column latlong with SRID 4326. The latlon - column has a GIST index defined. I would select all addresses within a distance given in metres.

I understand that SRID 4326 measures distance in degrees, so it doesn't make sense in my case. However the query

select *
from adresse
where st_distance(latlong::geography, ST_SetSRID(ST_Point(16.520, 47.846), 4326)::geography) < 200

is incredibly slow, I guess it is unable to use the index.

How can I speed up this query?


texts:
I have a number of data points in the form of (latitude,longitude,temperature). I am required to render the temperature on a map.

The final result should be similar to the the map here

I used google heat maps for this, and was able to render the data points successfully. The problem however is that the readings are far apart, so each reading appear as a disconnected circle. I could increase the radius of the rendered data points to solve this.

My question is: Is google heat map used to render the temperature correctly? or should I use another api?

Any advice is highly appreciated.


texts:
I have 3D shapefile which displays fine in QGIS 2D (and 3D in ArcScene), but I am having trouble changing to 3D view in QGIS. Shapefile is an imported polygon set (in a local CRS) > I get the message 'layer extent invalid'.
How do I fix this for my local CRS?


texts: ArcScene Extruding between constant TIN's to calculate bed volume within oil window?

texts:
I don't seem to be able to extrude between two TINs created from constant rasters.

Will the extrusion tool not work between TINs which are copied and pasted and altered in position by me after their creation?


texts: Align Center not aligning to page center in ArcGIS Pro after including margins from printer

texts:
After using the Include margins from printer option of ArcGIS Pro to display printer margins on a layout, the Align Center option stops aligning to the center of the page and aligns to the left printer margin instead.
To see this try the following using ArcGIS Pro 2.6.2:

Start ArcGIS Pro
From its Start Page use Start without a template
Insert a New Layout - I chose A4 from the ISO - Portrait category
Use the Straight text tool to drag a rectangle in a similar position to that in the picture below

Click on the page, away from the text, so that the text element handles can be seen
On the Text | Format tab of the ribbon, use the Align split button from the Arrange button group to first Align to Page, and then Align Center.  You should see the text move to the horizontal center of the page, which is to be expected.
Use the Undo button on the Quick Access Toolbar to move the text back to its previous position
Use the Layout tab of the ribbon to click on the Select Page From Printer button in the Page Setup button group (see picture below).  I am doing this because I want to display the printable area on my page.

Tick the Include margins from printer checkbox, and then OK.  You should see dashed lines appear just inside the page border to show the printable area of the page
Use the Text | Format tab of the ribbon again, to click Align Center again.  This time you should see the text move to be aligned with the left print margin (see picture below), which is NOT expected.

I believe this to be a bug, which I will report to Esri's local support.
While awaiting a fix for it, is there a workaround?


texts: Creating Postgis Polygon using ST_ConcaveHull from Geometry Points

texts:
I have a postgis 2.2 table with 20 columns of type geometry(Point,4326)

I'd like to generate a polygon which covers the outer boundary of the points - it seems like ST_ConcaveHull is a good option, but I can't see how to do it without first converting my points back to text (which seems to be missing the point).

Is st_concavehull the right option, and how do I go about constructing the query?


texts: Seeking literature on public transit network connectivity metrics?

texts:
Are you aware of any literature for measuring the connectivity of a public transit network, taking both route geometry and the schedule into account?

In other words, I'd like to design a process whereby one could feed in data about a transit network and its schedule, and the process would pop out some numbers that characterize the connectivity of the network - roughly, how convenient the connections are, how much time one has to wait to make a transfer from one route to another, or how likely it is that one can make a transfer from one route to another. Processes that require additional information, such as origin-destination demand or on-time performance, would be be interesting as well.

Such a process would be useful in quantifying the potential effects of a service change, in detecting and avoiding inadvertently breaking people's connectivity, and generally, in designing service to be optimally useful to passengers.


texts: What counts as a service request in ArcGIS for Server 10.3.1?

texts:
We've just updated our services to two new servers (Windows 2012) running ArcGIS for Server 10.3.1.

I'm trying to setup some statistics to monitor the usage on our new services, however I am having trouble finding any information in regards to defining what a "request" is.

I just finished publishing the last service, the environment has only been live for about 4 hours and it says there are already some folders with nearly 900 service requests.
If this number is accurate, great! I'm glad people are already using our services; though hard to believe...

So my main question here:
What signifies a service request?
Every time an image service is brought into arcmap? (this is hard to believe 900 people hit the same service within the last hour and a half)
Every time they zoom in or pan? (this is more believable if someone is panning around a lot and zooming on a single service)

Please help me better understand what these statistics are truly reporting!


texts: Inverse Distance Weighted (IDW) Interpolation with Radial Basis Function (RBF) in R

texts:
I want to perform an IDW interpolation with R. I am aware that I can do that with gstat::idw. But I want to change the parameter 'idp' with a radial basis function (RBF):

   where

 is the euclidean distance between of the i-th reference point and the j-th grid cell.

Do you know how I can implement this into the function of gstat::idw. Or do you know a script for IDW in R that I could edit?


texts: How to compute tabulate intersection in QGIS?

texts: I have a base layer of census data in census tracts (polygons). I'm looking to intersect those tracts with a buffer layer but I only want to include those census tracts whose area is over 50% within the buffer. The spatial query plugin doesn't offer that kind of precision and I can't find an equivalent of tabulate intersection. Is there any way to do this in QGIS?


texts: Converting EPSG:2263 to WGS84 using Python Pyproj

texts:
I'm trying to convert Projected bounds to WGS84 using Pyproj but getting strange results.

The Projected Bounds are : 909126.0155, 110626.2880, 1610215.3590, 424498.0529

The corresponding WGS84 Bounds (the target) are: -74.2700, 40.4700, -71.7500, 41.3100

The data is from: http://spatialreference.org/ref/epsg/nad83-new-york-long-island-ftus/

>>> from pyproj import Proj, transform
>>> inProj  = Proj("+init=EPSG:2236")
>>> outProj = Proj("+init=EPSG:3857")
>>> x1,y1 =  (110626.2880, 909126.0155)
>>> print(transform(inProj,outProj,x1,y1))
(-9122788.926222347, 3833518.013375292)

The outProj is based on my understanding of what Google Maps uses.


texts: Unable to perform centroids or point on surface operations because of "invalid geometry"

texts:
I am dealing with census blocks, and some of them have strange shapes, where they share a single point with two other census blocks.  Although check validity returned no issues, I am not able to create centroids with the census block features because of this geometry, because these operations return the "invalid geometry" error for those features. Note that the four small features that are (mostly) surrounded by the selected feature all have centroids, but the selected feature does not. The problematic point seems to be the one on the bottom left shared by the tear-drop shaped block, the long rectangle-shaped block, and the selected block.

Here is a link to the wkt for reference.

I have tried simplifying geometry, which returns the same "invalid geometry" error. I also tried v.clean.break with GRASS, but I am not able to create centroids with that output, and the view-only output produced by GRASS displayed "Inconsistent values of FID and field of same name" as an error when I attempted to save the output feature class. I'm very stumped. I probably should ignore this, but such geometry is not uncommon in US city boundaries, so I would like to make sure that I can create centroids with in QGIS.


texts: Adding new static page in CartoView / GeoNode?

texts:
I have installed the CartoView with the Windows stable build installer.

Then, I try to add a static page to the newly installed CartoView / GeoNode by modifying the urls.py file:

urlpatterns = [
# Static pages
url(r'^test/$', TemplateView.as_view(template_name='test.html'), name='test')
]

The test page entry doesn't work and localhost/test/ throws an error with no test page entry shown in the error log.

Has anyone successfully added a static page in GeoNode / CartoView?

Is there a compile action that I have been missing in forcing the urls.pyc update?


texts: Listing attributes of arcpy.Describe object in python?

texts:
I am trying to create a json structure that will describe a multitude of different arcpy data types (FeatureClass, FeatureDataset, Workspace, etc). I would like the script to be dynamic and give all possible information for the item that is passed into arcpy.Describe(), and since attributes of a Describe object change based on the input, I need a way to dynamically assign those attributes.

Here's what I've tried:

import arcpy
desc = arcpy.Describe("path\to\thing\being\inspected")

Attempt #1: (based on https://stackoverflow.com/questions/2675028/list-attributes-of-an-object)

desc.__dict__
# returns
AttributeError: DescribeData: Method __dict__ does not exist

Attempt #2: (based on https://stackoverflow.com/questions/2675028/list-attributes-of-an-object#comment19846153_2675148)

vars(desc)
# returns
TypeError: vars() argument must have __dict__ attribute

Attempt #3: (based on https://stackoverflow.com/a/687252/9052438)

dir(desc)
# returns
[]

Attempt #4: (based on https://stackoverflow.com/a/2675542/9052438)

desc.__dir__()
# returns
AttributeError: DescribeData: Method __dir__ does not exist

This all tells me that arcpy.Describe() isn't technically a Python object, but given that, is there still a way to do this without having to hard-code a list of attributes to check based on the dataType?


texts: Creating TINs in ArcGIS from triangular elements?

texts:
I am working with Finite Element Method (FEM) models with triangular elements, and I'd like to view the results as TIN surfaces in ArcGIS (and in ArcScene). In my external source, I have nodes (with X,Y,Z coords) and elements (forming a triangle from three nodes), which are the fundamentals of a TIN. How can I bring these together in ArcGIS to convert my triangular elements into a new TIN surface? This should ideally be a 1:1 conversion, but the documentation isn't helping me.

Here is an illustrated example of the data I have:

Nodes: n, X, Y, Z
1, 51.06, 84.98, -7.46
2, 30.45, 80.67, -7.47
3, 44.83, 72.70, 19.46
4, 64.99, 69.91, 61.48
5, 30.33, 63.40, 71.31
6, 46.95, 61.96, 48.44
7, 62.04, 51.97, 82.51
Elements: n1, n2, n3
6, 7, 4
4, 3, 6
3, 5, 6
5, 3, 2
1, 3, 4
2, 3, 1

I'm currently using ArcGIS 9.3 with the 3D Analyst Extension (among others), and would like some hints for a solution targeted for either Python or VBA.

Alternative approach: The raw data -> Esri TIN conversion doesn't look very easy, and I can't make heads or tails of the API for 3D Analyst, so I'll hold the Python/VBA solution off.

I have the same data represented in a Shapefile where all features are triangular PolygonZM shapes and the vertices have the height-field. (If you ask, this Shapefile was made using OGR via Python from the original dataset structured similarly as above). Viewing this file in ArcScene is really really slow, as I have ~100k triangular polygons in the Shapefile. This PolygonZM Shapefile should ideally convert triangle-by-triangle to a TIN dataset. I have the data, so interpretation or interpolation is not wanted!

So, how can I convert this PolygonZM Shapefile (pictured below) into a TIN? In ArcScene, in the "3D Analyst" menu, there is a "Create TIN From Features" tool, but I don't understand the language "soft clip" etc. The dialog provides no help.

LandXML

It turns out the best format is LandXML, thanks to the related answers below. For the example above, the file looks like this (note: the coordinate order is Y X Z):

<?xml version="1.0" ?>
<LandXML version="1.2" xmlns="http://www.landxml.org/schema/LandXML-1.2" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://www.landxml.org/schema/LandXML-1.2 http://www.landxml.org/schema/LandXML-1.2/LandXML-1.2.xsd">
  <Surfaces>
    <Surface name="My TIN">
      <Definition surfType="TIN">
        <Pnts>

84.98 51.06 -7.46

80.67 30.45 -7.47

72.70 44.83 19.46

69.91 64.99 61.48

63.40 30.33 71.31

61.96 46.95 48.44

51.97 62.04 82.51

        </Pnts>
        <Faces>
          <F id="1">6 7 4</F>
          <F id="2">4 3 6</F>
          <F id="3">3 5 6</F>
          <F id="4">5 3 2</F>
          <F id="5">1 3 4</F>
          <F id="6">2 3 1</F>
        </Faces>
      </Definition>
    </Surface>
  </Surfaces>
</LandXML>

Then use the LandXMLToTin_3d tool to process the result to an Esri TIN file:


texts: Project a set of coordinates - points from Excel to a shapefile with unknown coordinate system (map-match)

texts:
I am trying to project a set of coordinates (points) from an Excel spreadsheet, on top of a shapefile (map) with an unknown coordinate system.

The Excel file was created by first creating a map with some points in Google Earth and then exported it into a .kml file and then created the Excel.

The shapefile was given to me ready for the project I am doing. Basically I am trying to map-match the two files and then join-relate them to get some additional information for my project.

Here is what I did step by step:

Open ArcMap (V.10.1) and right clicked 'Layers' and defined a Geographic Coordinate System 'WGS 1984'
Imported the "vista_nicosia_links.shp" and received a warning for Unknown Spatial Reference

Exported the file into a shapefile using 'the same coordinate system as the Data frame' option.
Now in the Properties -> Source tab I can see the Geographic Coordinate System to be GCS_WGS_1984 and the 'Extend' info have 7 digits on the left of the comma.

I closed and re-opened ArcMap and performed actions of Step 1
Then I added X,Y data from the Excel file and then right-click -> display X,Y data and specified the X and Y field. In the description of Input Coordinate System it was by default GCS_WGS_1984.
I exported the data into a .shp file using 'the same coordinate system as the Data frame' option.
I re-opened ArcMap and repeated actions in step 1 and imported the two shapefiles I created in steps 3 and 7
The two shapefiles are displayed in entirely different positions the one far away from the other. I noticed that in right click -> Properties -> Source tab-> 'Extent' area of the Export_Output I see that there are 7 digits on the left of the comma (step 3), where as in Export_Output_2 I see only 2 digits on the left of the comma

Then I tried to project both shapefiles (one at a time) from "ArcToolbox -> Data Management Tools -> Projections and Transformations -> Define Projection" into Projected Coordinate system 'WGS 1984 World Mercator' but nothing seems to change or happening. They are both re-projected at the same (initial) positions.


texts:
I'm executing a long process WPS of around 180sec, and when it start after 60sec
the WPS (response) is abandoned, and the process continue to run until the end (checked on the log file).
I tried the solution proposed by Pywps wiki of reducing the timeout of
the Apache server (http://wiki.rsg.pml.ac.uk/pywps/Async_issue) without success !

I found the response which seems to be the only issue :
http://lists.wald.intevation.org/pipermail/pywps-devel/2013-April/001598.html
But I didn't find how can I insert my "call process" by self.cmd("process.sh...).


texts: Is there a way to set up data validation at the data entry level or restrict the entry to a choice set in CartoDB?

texts:
I’m wondering if there is a way to set up data validation at the data entry level or restrict the entry to a choice set in CartoDB.

Scenario: I want to display categorical data in a Visualization and continually update the database/Table. So… I add a point to my dataset/Table. Then, I click on it and edit the data. I now want to enter values for the various attributes for that point.

The problem: I can do it, but I don’t have any validation feature to restrict what I enter (other than if the attribute is set as “number”, “boolean", or “date”). I want to limit what I can enter to a choice set (e.g., like in with drop-down menu). If I could have the world, so to speak, I would love to have a survey form type of data entry feature – complete with “multiple choice”, “select all that apply”, formatting checks, and maybe even logic.

Question 1: Is this possible with CartoDB?

Question 2: Is there a way I could set this up, e.g., with embedding a survey form in the map pop-up window? I want the data to go straight into the CartoDB dataset/table that drives my Visualizations. I suppose I may be able to do this directly using a survey and have CartoDB reference or sync with that dataset/Table… but I want to be able to manually locate/pin the point on a map (to get Lat/Lon), and I don’t know of a survey platform that is capable of doing this.


texts: to be able to use mysql-connector-python inside QGIS python script

texts:
I am creating python script in QGIS python console to be able to connect to mysql db and manipulate data (retrieve from db, generate shapefile and show it as a layer).
I tried to install mysql connector from OSGeo4w shell:
pip install mysql-connector-python but it fails with error:

C:>pip install mysql-connector-python
Collecting mysql-connector-python
C:\OSGEO4~1\apps\Python27\lib\site-packages\pip_vendor\requests\packages\urllib
3\util\ssl_.py:318: SNIMissingWarning: An HTTPS request has been made, but the S
NI (Subject Name Indication) extension to TLS is not available on this platform.
 This may cause the server to present an incorrect TLS certificate, which can ca
use validation failures. You can upgrade to a newer version of Python to solve t
his. For more information, see https://urllib3.readthedocs.io/en/latest/security
.html#snimissingwarning.
  SNIMissingWarning
C:\OSGEO4~1\apps\Python27\lib\site-packages\pip_vendor\requests\packages\urllib
3\util\ssl_.py:122: InsecurePlatformWarning: A true SSLContext object is not ava
ilable. This prevents urllib3 from configuring SSL appropriately and may cause c
ertain SSL connections to fail. You can upgrade to a newer version of Python to
solve this. For more information, see https://urllib3.readthedocs.io/en/latest/s
ecurity.html#insecureplatformwarning.
  InsecurePlatformWarning
  Could not find a version that satisfies the requirement mysql-connector-python
 (from versions: )
No matching distribution found for mysql-connector-python

What is the best way to solve this?


texts: How to extract elevation data from binary .ADF elevation format at given coordinates?

texts:
I have a high resolution digital elevation model available. The ground resolution is 1 meter and the file format is adf which is some binary ESRI grid raster format.

I have a script the outputs x/y- or lat/lon-coordinates and as a next step I have to extract the elevation of the given coordinate from the binary adf file. But that's a huge one. The file is like 1.5 GB in size and I haven't found a way yet how to access the data via a headless command line or script.

How to access adf files via script? Is there any programmes available or are there any known libraries I could use?

Or should I rather convert the adf to some ASCII text files? But that will cost a lot of disk space and performance.

I've been reading about elevation databases. How to create one from my elevation model? Is it worth a try regarding performance?

Update: To break this down, in short - How to access adf files or elevation data most efficient using a script with coordinates as parameters?


texts: Leaflet Opacity Control Using Slider For Bootstrap

texts:
I am attempting to use the Slider For Bootstrap to update the opacity of a layer within Leaflet.

The slider also needs to be in a Bootstrap Popup but this part is working.

As per this JSFiddle you can see I have the slider working but I cannot get it to update/ set the opacity of the background layer, how can this be achieved?


texts:
I installed the new version of QGIS, QGIS 3.14 and I am having problems in running GDAL algorithms through Processing.
For instance "clip raster by mask layer" always gives me:
Loading resulting layers
The following layers were not correctly generated.
• C:/Users/../AppData/Local/Temp/processing_CBxlhD/70b11f1d802a45ffafcfad0804ef2610/OUTPUT.tif
You can check the 'Log Messages Panel' in QGIS main window to find more information about the execution of the algorithm.

Investigating more the LOG, I have this error:
ERROR 1: Failed to identify source layer from datasource.

How to solve it?


texts: Feature Class to Feature Class Keeps Layer Name In Fields 10.2?

texts:
I have a model that makes a join, adds two fields, and calculates the fields, then uses Feature Class to Feature Class to create a new output. The new output is named "Permit." Once the output is created the fields are labeled "Permit_" + the field name (from source data). I want to exclude the layer name ("Permits_") from the fields, and only keep the field name from the source data. Any thoughts on how to do this? (I have to use Modelbuilder.)

I've traced back where the "Permit_" has been added to my field names. Prior to creating the join, I have a Make Feature Layer, called "Permits Layer". After this there is the join, the fields are then called "Permit.FieldName", this is carried to the final output.

Image of Make Feature Layer joined to source data:

Image of my final output:

Image of entire model:


texts: Using an attribute index to find matching attributes of two layers faster?

texts:
Similar to Indexing attribute field of shapefile in QGIS, I am wondering if such a thing like an attribute index exists for PyQGIS. Goal of its usage would be to iterate over two vector layers and find matching attribute values of a specified field in each layer. So it would work like a spatial index, just using attributes instead. So far I could only find, that I can create an Index using createAttributeIndex() as stated here and here. But absolutely no further information about its usage, the way it works or examples.
Basically the idea is to speed up code written like this:
vectorlayer_a = QgsProject.instance().mapLayersByName("layer_a")[0]
vectorlayer_b = QgsProject.instance().mapLayersByName("layer_b")[0]

for feat_a in vectorlayer_a.getFeatures():
    value_a = feat_a.attribute(1)
    for feat_b in vectorlayer_b.getFeatures():
        value_b = feat_b.attribute(1)
        if value_a == value_b:
            print('Hurray, finally found (another) one. Can I find all of them faster with an attribute index?')
            # Do stuff some stuff like...
            geom_a = feat_a.geometry()
            geom_b = feat_b.geometry()

Also, could attribute(1) have any datatype or would such a thing only work with numerical values, if this 'thing' exists at all?


texts: Administrative level GeoJSON Polygons for African Countries

texts: I'm very new to GIS and I have spent the last few days trying to find some kind of database with the aforementioned polygons with no luck. The purpose is simply to label a location at district level given its longitude and latitude. I tried reverse geocoding but very quickly exceeded my rate-limit because of the size of the data. I've also tried GADM but the data is not labeled and there seemed to be more polygons than districts. MongoDB works nicely with GeoJSON hence my preference for data in this format. What options do I have?


texts: Determining if PDF has geospatial information?

texts:
Asked to convert geoPDF to geoTiff. I am familiar with Esri's tool:

http://desktop.arcgis.com/en/arcmap/10.3/tools/conversion-toolbox/pdf-to-tiff.htm

However, I can't seem to convert these specific PDFs. So if within a PDFs properties "Tagged PDF" is unchecked, does that mean it is lacking geospatial tags?(see image below)

Do I need a TerraGo extension to accomplish this?


texts: How to scrape data from a google map (in flash)?

texts:
I've contacted the author of this NYC homicide map from the NY Times several times to see if he could send me a spreadsheet of information the data is based off of, but I have not had any response (for a very long time).

Is there anyway I could scrap the data from the google map directly? I downloaded the swf file and converted to xml, but even within there I can not tell where the information is I want. Here is the zipped up file of the converted xml.

I'm hopeing someone with more knowledge of online mapping can quickly tell me if this possible or not (or where I should be looking in the xml file).


texts: Authentication methods missing when connecting to ArcGIS FeatureSevice in QGIS 3.14 (Pi)

texts:
I am trying to connect to an ArcGIS-10.3 FeatureService in QGIS 3.14 (Pi) recently installed.
In the connection details window, I have the URL but in the "Authentication" part, in "configuration", only "no authentication" is available. There is no drop down options available. I need an OAuth2 authentication (the plugin is installed - I checked in the Authentication window in Settings - options - authentication, installed plugins).


texts: In to_crs of GeoPandas, init proj4 string works, but new style doesn't

texts:
I wanted to unify crs into WGS 84, so based on
"Manually specifying the CRS"
I tried gdf.to_crs('epsg:4326'), but error occurs.
I want to know why the old proj4 style gdf.to_crs({'init':'epsg:4326'}) works, why not the recommended gdf.to_crs('epsg:4326'). I used python 3.8.1 32bit, geopandas 0.7.0 and pyproj 2.5.0.

import geopandas

gdf=geopandas.read_file('note/map/bnd_emd_2012/BND_ADM_EMD_2012.shp')
gdf.crs

<Projected CRS: PROJCS["PCS_ITRF2000_TM",GEOGCS["GCS_ITRF_2000",DA ...>
Name: PCS_ITRF2000_TM
Axis Info [cartesian]:
- [east]: Easting (metre)
- [north]: Northing (metre)
Area of Use:
- undefined
Coordinate Operation:
- name: unnamed
- method: Transverse Mercator
Datum: International Terrestrial Reference Frame 2000
- Ellipsoid: GRS 1980
- Prime Meridian: Greenwich

gdf1=gdf.to_crs({'init':'epsg:4326'})

Warning (from warnings module):
  File "C:\Python\lib\site-packages\pyproj\crs\crs.py", line 55
    return _prepare_from_string(" ".join(pjargs))
FutureWarning: '+init=<authority>:<code>' syntax is deprecated. '<authority>:<code>' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6

gdf1.crs

<Geographic 2D CRS: +init=epsg:4326 +type=crs>
Name: WGS 84
Axis Info [ellipsoidal]:
- lon[east]: Longitude (degree)
- lat[north]: Latitude (degree)
Area of Use:
- name: World
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

gdf2=gdf.to_crs('epsg:4326')

Traceback (most recent call last):
  File "<pyshell#26>", line 1, in <module>
    gdf2=gdf.to_crs('epsg:4326')
  File "C:\Python\lib\site-packages\geopandas\geodataframe.py", line 564, in to_crs
    geom = df.geometry.to_crs(crs=crs, epsg=epsg)
  File "C:\Python\lib\site-packages\geopandas\geoseries.py", line 436, in to_crs
    transformer = Transformer.from_crs(self.crs, crs, always_xy=True)
  File "C:\Python\lib\site-packages\pyproj\transformer.py", line 319, in from_crs
    _Transformer.from_crs(
  File "pyproj\_transformer.pyx", line 311, in pyproj._transformer._Transformer.from_crs
pyproj.exceptions.ProjError: Error creating Transformer from CRS.: (Internal Proj Error: proj_create_operations: geodeticCRS not found)

gdf2.crs

Traceback (most recent call last):
  File "<pyshell#27>", line 1, in <module>
    gdf2.crs
NameError: name 'gdf2' is not defined


texts:
I'd like to create several images of voronois polygons in QGIS; each images having differents values joined to the voronois.
I tried to load a qml style with the values I want to keep (8 classes : 0>700, 700>1400, 1400>2100, ..., 4900>5600) for all my images, but it doesn't work when I pick the column of values I'd like to classify into the voronois: the maximum value of the picked column replaces then the one I'd like it to be the maximum (5600).
Would you know a way to apply that same classification to the differents values I have? For instance through fixed classes?


texts: Performing cluster analysis in R with Point Datasets

texts:
I am a R-Beginner, but I want to use it for some spatial statistics. I want to check the relation between the occurrence of ill persons and the existence of water in a area of interest. To give you a more specific use case, I've created some example datasets:

Table - Persons
Table - City
Table - Water

Firstly I want to plot them in a map including a basemap. This already works!

persons_dataframe = read.csv("persons.csv", header = TRUE)
city_dataframe = read.csv("city.csv", header = TRUE)
water_dataframe = read.csv("water.csv", header = TRUE)
# plot them on a map
# load used packages
library(RgoogleMaps)
library(ggplot2)
library(ggmap)
library(sp)

persons_ggplot2 <- persons_dataframe
city_ggplot2 <- city_dataframe
water_ggplot2 <- water_dataframe
gc <- geocode('new york, usa')
center <- as.numeric(gc)
G <- ggmap(get_googlemap(center = center, color = 'bw', scale = 1, zoom = 10, maptype = "terrain", frame=T), extent="device")
G1 <- G + geom_point(aes(x=POINT_X, y=POINT_Y ),data=city_dataframe, shape = 22, color="black", fill = "yellow", size = 4) + geom_point(aes(x=POINT_X, y=POINT_Y ),data=persons_dataframe, shape = 8, color="red", size=2.5) + geom_point(aes(x=POINT_X, y=POINT_Y ),data=water_ggplot2, color="blue", size=1)
plot(G1)

This has the following output:

No I want to create something like a cluster map. I've already looked for some useful packages and I've ended with the packages SpatialEpi and DCluser. Both provide some functions, which might be useful.

But at this point I don't know, how to proceed. Some open questions:

Is it right that I need polygons to cluster the dataset? I've read it in the SpatialEpi pdf but I only have these point datasets.
Additionally I don't know, if the fact that all points are located on only three locations might be a problem?!

Do you have a suggestions how to continue with my work?


texts: Selecting entities based on first character of field being a letter using QGIS?

texts:
I would like to select entities whose addresses are incomplete, i.e. the first character is a letter instead of being a number. Does QGIS have an "is numeric" -type function to accomplish this?

Not Is Numeric (left("ADDRESS",'1') )


texts: Can someone explain to me what the "mini dump" error mean in QGIS?

texts:
I am struggling to understand what the error "Crash dumped: mini dump written to C:" means in QGIS. I am working with layers stored on a network drive. Whenever I try to read in a few particular .asc files QGIS returns the above error and closes.

These layers were working fine until I did some processing of one of them in SAGA and saved out to a new file. Now I can't read any of the original or processed layers, or any new .asc files exported from SAGA for that matter. The original layers still work fine on all other computers in the department, so they are not corrupted, it is something to do with my version.


texts: Generating Individual Shapes from Points with values

texts:
I have this shapefile of points, categorized in 8 colors.

what i want to do, using python and shapely is, for each category create a shapefile (fill polygon), but when i try this, i get this result.
(The points around the polygon i made on photo editor to ilustrate)

its possible to get this individuals shapes?


texts: Making polygon shapefile from text file with arcgisscripting?


texts: Click events in wrong place after div is moved when using MapView in version 4.1

texts:
I am having a problem when the div that a MapView is rendered in moves.  After moving the div, mouse events do not translate properly to the correct place on the map/screen.  These events appear to be translated based on the div’s original position.  It looks like the MapView is storing the div’s original position when it is created.  I am wondering if there is any way to manually update this after the MapView is created.  I do not see any likely methods in the api and have not been able to find any other similar issues.

You can recreate this issue by pasting the following code into the ArcGIS Sandbox.  My recreation is based on ESRI’s “Get started with MapView - Create a 2D map”.  After pasting the code and pressing Run, zoom into the map by holding the shift key and drawing a box on the map with your mouse.  After confirming that everything is okay, tap the “Scoot View” button a few times and try it again to illustrate the issue.  You will see that the zoom box does not appear where your mouse is.

I am working on an application where I would like to embed a map on a dashboard made up of movable widgets.  I could tear-down and recreate the MapView every time the widget is moved but that seems extreme.

<!DOCTYPE html>

    <meta charset="utf-8">
    <meta name="viewport" content="initial-scale=1,maximum-scale=1,user-scalable=no">
    <title>Get started with MapView - Create a 2D map - 4.1</title>

    <link rel="stylesheet" href="https://js.arcgis.com/4.1/esri/css/main.css">

<button onclick="scootView()">Scoot View</button>


texts: Problem with panning and zooming in OpenLayers

texts:
I'm interested in detecting when a user moves or pans a map so that I can run a new query on a database based on their location. However, the following code causes the panning and zooming controls to go crazy:

map.events.register("moveend", map, function(){
    getValues(map.center);
    });

If I comment out these lines, there appears to be no problem. If they are included, after two pans or zooms, the map suddenly jumps to half-way around the world (it looks similar to 0 lat, 0 long). I'm unsure what I am doing incorrectly.

The other problem I have regarding the pan/zooming is that the movend does not really seem to wait until the user has finished moving or panning. Queries seem to be fired off even as the map is being moved and panned. Again, I'm not sure if this is a problem with how I am detecting this but it does not seem to behave as I thought it would.

EDIT: Including getValues function also:

getValues = function(c){

    c.transform(map.getProjectionObject(), new OpenLayers.Projection("EPSG:4326") )
    console.log(c.lat, c.lon)

    $.ajax({
            type: "POST",
            url: "php/query.php" ,
            data:{x: c.lat, y:c.lon},
            async: true,
            success: function(r){
                    console.log(r);
            },
        dataType:'json'
        });

    }


texts: Automate database conversion from NAD83 to Web Mercator

texts:
A NAD83 database is updated daily while a Web Mercator clone is only update when time is available. Looking to automate this process of

Copying Layers
Projecting them to Web Mercator
Placing in correct .SDE

What would be the best way to do this? FME, model builder, python?

Thanks


texts: Installing Esri FGDB API with GDAL (Python 3.5)

texts:
I've recently migrated to ArcGIS Pro 1.4 and have had no problems moving from Python 2.7 to 3.5 - except one.  I try to use gdal/ogr when I can in order to share scripts easily with those who do not have an Arc license available.  In the 2.7 version, I could easily obtain support for read/write of Esri file geodatabases by placing the FileGDBAPI.dll into

Python27\ArcGIS10.4\Lib\site-packages\osgeo\gdalplugins

However, with ArcGIS Pro (Python 3.5) when I installed the gdal (2.1.3) wheel (obtained from Christoph Gohlke's site at http://www.lfd.uci.edu/~gohlke/pythonlibs), there is no gdalplugins folder.  I only have the following structure:

C:\Users\[my username]\AppData\Local\Programs\ArcGIS\Pro\bin\Python\envs\arcgispro-py3\lib\site-packages\osgeo

Under this, I have:

\__pycache__
\data
\include
\lib
\scripts

There are only 2 .dll files in the entire structure; gdal201.dll and geos_c.dll.

Could someone tell me where the file geodatabase drivers reside, and if adding Esri's FileGDBAPI.dll to the same directory is all that's needed to get things working?

My limitations are: Windows 7 64-bit, I need to be able to deal with pre-compiled binaries, I do not have a compiler available to me or full administrative rights on the machine.


texts: Getting error 'exceptions.NameError' in Python script?

texts:
I'm getting <type 'exceptions.NameError'>: name 'B' is not defined when I run this script.

Why is my 'B' not defined and is it a syntax error?

For clarification, I am trying to get the code to loop through specific shapefiles in the folder "\\My Documents\ArcGIS\..." and copy them into "\\My Documents\NewFolder\..."  The A variable represents a list of the specific shapefiles I want to copy from a folder and into a new one; the issue I am having is defining B, where the individual variable name from A is automatically carried over/read in any part of the code.  I am trying to tell Arc that every time it sees "+STATE+" (which should be seen as B: main issue) it must insert the specific variable name listed in A, process it (in this case copy), and then loop back to the next specified variable, until it reaches the end of the list.

I would like to copy the 'FL_RechargeArea' through 'GA_AreaofContribution' shapefiles from the "\\My Documents\ArcGIS\..." folder into the "\\My Documents\NewFolder\..." and get it to save as '(specific-varialbe-name)2.shp'.

It's much less about copying and more about being able to loop through and process specific shapefile names in a folder, rather than looping through and processing perhaps the entire list within a folder.

# Import arcpy module
import arcpy

for A in B:

    # Local variables:
    A = ['FL_RechargeArea','FL_AquiferProtectionArea','FL_AreaofContribution','GA_RechargeArea','GA_AquiferProtectionArea','GA_AreaofContribution']
    B = "STATE"
    New_Shapefile = "\\\\My Documents\\ArcGIS\\"+STATE+".shp"
    Output_from_New_Shapefile = "\\\\My Documents\\NewFolder\\"+STATE+"2.shp"

    # Process: Copy Features
    arcpy.CopyFeatures_management(New_Shapefile, Output_from_New_Shapefile, "", "0", "0", "0")


texts: OSRM-extract Sorting edges by start seems to stop

texts:
I try to use planet-latest.osm.pbf (2018-05-23) from https://planet.openstreetmap.org/pbf/.

I extract the pbf file on Ubuntu 16.04 and everything seems OK until Sorting edges by start ...

On internet, I saw that this operation takes usually only few seconds. But my server has been running it for more than 5 hours.

Should I just wait or has something gone wrong ?


texts: What is the case for Building Information Modeling (BIM) in GIS?

texts:
Working for a consulting firm has introduced me to many hot topic technologies, and key words/lingo. One of which is BIM, or building information modelling. It is definitely a popular topic among architects, structural and civil engineers that praise its utility, as well as software developers and architects that are data-savvy. There is definite talk about developing GIS workflows to support BIM project deliverables.

What it the case for BIM in GIS and where do we currently stand? How is BIM slated to improve GIS workflow, quality of deliverables, etc. Since when did GIS analysis become concerned with the "inside of buildings". Traditionally speaking, autocad/microstation technicians were left to design the interiors, and GIS was used for those exterior applications (zoning, land use managements), macro in retrospect.


texts: Geoserver 2.4.3 / OpenLayers - setting zoom level dynamically

texts:
Expected user behavior is that when user selects the County from a drop-down, the map should zoom to its max extents. I am using SQL View (Postgis on 9.3) and since I do not know the bounds, I gave the Max Bounds of the FL State. Not able to create layer without giving bounds. Can the boundaries be inferred from the SQL? The SQL for the view is as shown below:

   select ST_XMin(b.geom) as xmin,ST_YMin(b.geom) as ymin, ST_XMax(b.geom) as xmax,
   ST_YMax(b.geom) as ymax, a.zcta5ce10, a.geom
   from us_zip a, tl_2013_us_county b
   where b.statefp='12'
   and lower(b.name)='%COUNTY_NAME%'
   and ST_contains(b.geom,a.geom) order by zcta5ce10

If it is not possible, I'm thinking of getting the max bounds for the selected county before invoking the OpenLayers map.


texts: ogr.geometry() and ogr.GetGeometryRef() differences

texts: Can you tell me the difference between ogr methods ogr.geometry() and ogr.GetGeometryRef()? Is there any difference in returned object in Python?


texts: Making table that displays overlaying shapefile information?

texts:
I have a parcel I'm working with that has multiple layers laid over it. I'd like to make a table that displays various attributes of those different layers. For example, within that shapefile, I'd like to have a table display the watershed, school district, utilities, etc. and if that parcel changes, I'd like the table information to change too. I've experimented with data driven pages, but I'm not sure that's quite what I'm looking for.

Any suggestions?


texts: Coordinate system doesn't contain values in the range -180 to 180 and -90 to 90

texts:
I'm creating an android version of a PC program (I've contacted the company who owns the PC program, so it's legal). The program is in the core a GPS, but is used to navigate pre-defined routes for heavyweight transports, because they are not allowed to pass some bridges. However, that isn't very important for my problem. The program reads the coordinates from text files, so my app needs to do that too, but those coordinates are not in the range -180 to 180 and -90 to 90. So I can't really understand them. For example, one coordinate has these values:

X: 6603939,57000491
Y: 575841,753705909

I need to convert these coordinates to WGS84 coordinates (lon and lat) in order to be able to use them for what I want.


texts: Shattered, rough WMTS image in QGIS 2.18.12

texts:
After updating QGIS i found that wmts layers are being rendered in way that makes them useless. It look's like tiles are downloaded in greater resolution and then scaled down what makes them rough and all captions unreadable.

I'm on windows 10.


texts: QGIS Atlas Scale Bar in Miles with Round Numbers

texts:
I am working in QGIS 2.16.2 on an atlas comprising approximately 60 locations in 7 states from Texas to Florida to Ohio. Each location has a different coverage area that scales based on the size of circle polygons in a hidden coverage layer.

I would like my scalebar to show distances in miles because it is more intuitive to the client than meters or km. I have accomplished this by setting the units to feet with a label unit multiplier of 5280. Segments are set to "right 4" with fit segment widths of "min 50 mm" and "max 80 mm".

The problem I am having is that the divisions in the scale bar divide to a nice, round number before applying the multiplier (5280), so I get divisions like 1.51515 miles. I currently have the project projection in 3857 (pseudo-Web Mercator). I am open to other projections, and considered 2163: US National Equal Atlas Area, but that is still based on meters and I have the same issue. I could try a state plane based on feet, but I am not working in one particular state and I still have the 5280 multiplier. I am not aware of any projections with miles as units.

I have tried setting the map measurement units to miles under the Project Properties --> General menu, to no avail. If I do a fixed width segment of 5280 units on the scale bar, I get a nice result but it is always limited to 4 miles, even if the map is 400+ miles across.

Is there a way to get nice, even numbered divisions (like 1, 1.5, 2) while using a label unit multiplier like 5280? Or is there another way to have the automatically generated scale bars show miles nicely?

What I have:
  or

What I want:


texts: How to create ArcPy script that checks for duplicates and warns if adding one?

texts:
I am new to GIS(Arcpy)/Python programming , but I have experience in other languages. I have a question regarding the possibility to create a tool/script that does the following

-- Check a specific attribute field and see if there are any Duplicates (Our problem is that our DOC_ID field is supposed to be a unique value but some interns are accidently pasting/writing the same DOC_ID from a previous edit and saving it.

*The checking for duplicates is not that important as I have create a seconday field to check the data within the DOC_ID field to find duplicates*

-- Is it possible for the tool / script to be continuous. What i mean for this is that if a user accidently inputs a DOC_ID that is already within the list ... some kind of error or warning message would pop up saying that the "DOC_ID has already been used".

  ^(This is the main feature that I am trying to create or solve.)

If this is not possible, are there any other paths that I should consider.
(Validation Rules , etcetc)

This is not a big problem currently. There are only around 100-200 duplicate entries so far. We are trying to avoid the problem the best we can early once we start to implement more and more entries.


texts: Error when save feature attribute with Openlayers

texts:
I try OpenLayers.Control.GetFeature to get feature from WFS layer :

  function editTool() {
        var edit = new OpenLayers.Control.GetFeature({
            protocol: OpenLayers.Protocol.WFS.fromWMSLayer(Vidagis.CurrentMap.map.layers[1]),
            click: true,
            box: true,
            multipleKey: "shiftKey",
            toggleKey: "ctrlKey"
        }, {
            displayClass: 'pan no-follow',
            allowDepress: true
        });
        edit.events.register("featureselected", this, function (e) {
            currentFeature = e.feature; // get the current handler feature
            changeHtmlForm(storeFeature(e));
            unlockUpdate(storeFeature(e).data);
            $(function () {
                $("#smalldialog").dialog({
                    position: { my: 'right', at: 'right' },
                    maxHeight: 400,
                    maxWidth: 300,

                    autoOpen: false,
                    show: {
                        effect: "blind",
                        duration: 1000
                    },
                    hide: {
                        effect: "explode",
                        duration: 1000
                    }
                });
                $("#smalldialog").dialog("open");
            });

        });
        return edit; // return edit Control
    }

it work fine and i can get feature attributes from featureselected event .now i assign values to currentFeature

  currentFeature = e.feature;

Then i write a function to edit currentFeature attributes

 // Edit feature attributes
    function setFeatureAttributes(e) {
         e.attributes['code_2004'] = 124;
         return e;
    }

finally , save function

   function UpdateFeature() {
        var feature=setFeatureAttribute(currentFeature);
        feature.state = OpenLayers.State.UPDATE;
        saveStrategy.save([feature]);
    }

Declare variable :

  var currentFeature;
    var saveStrategy = new OpenLayers.Strategy.Save();
// trigger save event
    saveStrategy.events.on({
        'success': function (event) {
            alert('Changes saved');
        },
        'fail': function (event) {
            alert('Error! Changes not saved');
        },
        scope: this
    });

But when i call UpdateFeature() browser report the error

  Uncaught TypeError: Cannot read property 'projection' of null

How can i fix it?


texts: Generating GeoJSON tiles from GeoDjango - basic overview?

texts:
I'm working with the glorious PolyMaps and I want to build a map similar to the population example.

I understand that this needs GeoJSON tiles. I have a GeoDjango database, with a bunch of County objects, each of which has a Polygon field - these are what I'd like to display on the map.

Can anyone give me a high-level overview of how to turn these County polygons into GeoJSON tiles?

I've seen this question recommending TileStache, but is TileStache the most sensible thing to use with GeoDjango? And if so, how would TileStache fit in?

Bonus points for explaining what a GeoJSON tile, um, is.

Thanks :)


texts: Supervised Classification using Remote Sensing data?

texts:
Is there a correlation between the time that supervised classification takes to execute and the training areas that we choose?

Perhaps if we take larger areas for classification will the classification take more time, or if we take multiple small urban areas instead of a big one will it take less time.

Because I am trying to classify a MODIS Image with knn and it takes about 30-40 seconds and I don't know if my training areas are too big or if I should change to multiple small training areas and I don't really know where to look for to get a straight answer.


texts: ArcGIS 10.2/10.2.1 for Desktop with Data Interoperablilty not seeing GPX files?

texts: I have ArcGIS 10.2 and 10.2.1 both with Data Interoperablilty extension. But there is a problem with GPX files. ArcGIS simply doesn't see it. In 10.1 there was no problem to drag n drop them but now it's impossible. In quick import tool, when I open FME reader gallery I can't find gpx on the list. What's going on? Have you the same problem?


texts:
Sorry for a very basic questions about rtree.

I would like to check whether the point is inside the polygon and the distance (point to the nearest point in polygon) as well. I was thinking about using No-SQL database, but both MongoDb and GeoCouch support Polygon records, but don't provide searches like:
- what polygons are around the point
(for example, I want to find near what polygons a user is currently located; and after that check what in the precise distance from the user to the nearest polygon point)

So, I want to figure out how to make these kind of spatial searches using python and related libraries only (shapely, geopy, rtree).
Here are the steps I am following:

I have shapefile with the targeted polygons (saved from QGis)
I want to make Rtree Index, to improve performance when working with Shapely and geopy

I am wondering, whether I can:

make Rtree Index for shapefile directly?

If not. I can use shapefile lib to collect polygon coordinates or bounding boxes:

>>> shapes[33].points
[[120.9937913, 24.7902669], [120.9939797, 24.790369], [120.994035, 24.7902849], [120.9938467, 24.7901828], [120.9937913, 24.7902669]]

>>> shapes[33].bbox
[120.9937913, 24.7901828, 120.994035, 24.790369]

Now I am a bit confused, and trying to understand how to create rtree index file for my polygons:

>>> idx = index.Index()  //creating index
>>> left, bottom, right, top = (24.7901828, 24.790369, 120.9937913, 120.994035) // can I put bbox coordinates here? I assume that coordinates are in the form [xmin, xmax, ymin, ymax]
>>> idx.insert(0, (left, bottom, right, top))

Would it be a right way for me to make Rtree Index?:

import shapefile
from rtree import index
sf = shapefile.Reader("polygons.shp")
idx = index.Index()
shapes = sf.shapes()
id_num = 0
for shape in shapes:
    sh = shape.bbox
    xmin = sh[1]
    ymin = sh[0]
    xmax = sh[3]
    ymax = sh[2]
    idx.insert(id_num, (xmin,ymin,xmax,ymax))
    id_num = id_num + 1

After I execute the script above, I can get the following result from idx.bounds:

>>> idx.bounds
[24.7796055, 120.9839184, 24.7986571, 121.0025152]

To what bounds these coordinates are related to? Is it the bbox for the all the polygons?

Is this a right way to check whether there are polygons near the Point?
Assuming that the point coordinates are 24.789077, 120.999792

>>> list(idx.nearest((24.789077, 24.789077, 120.999792, 120.999792)))
[1, 2, 13, 29, 61, 62, 59, 60, 27, 28, 55, 5, 26, 54, 78, 49, 50, 3, 9, 10, 22, 72, 43, 39, 42, 85, 40, 8, 17, 76, 77, 36, 73, 71, 33, 34, 69, 63, 66, 64]

I am getting these results and not quite sure what they mean.. What does "nearest" exactly mean here? I cannot measure the distance here, like 5 meters.. so what is the "nearest" output here and how it's calculated?


texts: Optimizing minimum distance calculations between multiple polygons in R

texts:
I've already asked this question on stackoverflow, but here it seems more relevant.
I have a list of polygons and would like to create a minimum distance matrix by calculating the minimum distance between each possible pair. I am currently using the code below, but this is running really slow.
Is there a way to optimize this?
list #dataframe with all my multipolygons
distance<-as.data.frame(matrix(ncol=nrow(list),nrow=nrow(list)))#create an empty dataframe which I can then fill with distances.

for (i in 1:nrow(list)){
  for (j in i+1:nrow(list)){

      distance[j,i] <- st_distance(x = list[i,], y= list[j,], by_element = TRUE)

  }
}


texts: What geoprocessing tool do I use in order to assign a projected coordinate system to a raster?

texts:
I want to resample a raster to a smaller cell size. The raster cell size is given in decimal degrees with a GCS of North America 1983. The original raster covers the continental US, and I have clipped it to a state level. I would like to be able to resample according to a 30 m cell size, with a projection for the state.

Should I use the Project raster tool? I only see information on the GCS, so I am not sure if there is a planar projected coordinate system.


texts: Error Label with Python: encoding declaration in Unicode string (, line 0)

texts:
I get this error using Python while I'm trying to label a shapefile. The strange thing is that for other shapefiles the code works.

There is a table (relates.dbf) that has the data to be used for labelling and it has the fields: Schacht (string), Auswahl (short integer), max_Übers (double) and Max_Wieder (double).

The problem shouldn't be there.

NB: The SQL for this labelling class is [TN_UEB] > 0 and I use ArcGIS 10.3.1. The code is a variation of http://support.esri.com/technical-article/000012230

Any clues? What does it mean "line 0"! Where's the problem?

def FindLabel ([Name], [TN_UEB]):
  import arcpy

  key1 = [Name]
  key2 = "Schacht"
  key3 = "Auswahl"
  L =  "<FNT size = '10'><BOL>" + [Name] + "\n" +  " Ta = " + [TN_UEB] + " a" + "</BOL></FNT>" + "\n"
  L2 = "max_Übers"
  L3 = "Max_Wieder"
  myDataTable = r"P:\Duesseldorf\50009\FOG\Bestand_pz\Vergleich_Ergebnisse\KOSTRA_Auswertung\Variante_1_D_neuneu\relates.dbf"

  cur = arcpy.da.SearchCursor(myDataTable, [key2, L2, L3, key3])
  for row in cur:
    if str(key1) == str(row[0]) and row[3]>0:
      L = L + str(round(float(row[1]),1)) + " m" + "<SUP>"+"3"+"</SUP>" + " (T=" + str(round(float(row[2]),1)) + "a)" + "\n"
  return L


texts: How can filter low cloud coverange in a region with Sentinel-2 L2A in GEE?

texts:
I'm trying to filter clouds in Sentinel-2 through GEE in a ROI but I get the same number of images clipping or not clipping the geometry.
How can I calculate the cloud coverage in a ROI and filter an ImageCollection depending on the cloud coverage?
This is the code I was working on:
var s2Sr = ee.ImageCollection('COPERNICUS/S2_SR');
var s2Clouds = ee.ImageCollection('COPERNICUS/S2_CLOUD_PROBABILITY');

var START_DATE = ee.Date('2018-03-30');
var END_DATE = ee.Date('2018-06-10');
var MAX_CLOUD_PROBABILITY = 20;
var region = AOI
Map.centerObject(region, 12);

function maskClouds(img) {
  var clouds = ee.Image(img.get('cloud_mask')).select('probability');
  var isNotCloud = clouds.lt(MAX_CLOUD_PROBABILITY);
  return img.updateMask(isNotCloud);
}

// The masks for the 10m bands sometimes do not exclude bad data at
// scene edges, so we apply masks from the 20m and 60m bands as well.
// Example asset that needs this operation:
// COPERNICUS/S2_CLOUD_PROBABILITY/20190301T000239_20190301T000238_T55GDP
function maskEdges(s2_img) {
  return s2_img.updateMask(
      s2_img.select('B8A').mask().updateMask(s2_img.select('B9').mask()));
}

// Filter input collections by desired data range and region.
var criteria = ee.Filter.and(
    ee.Filter.bounds(region), ee.Filter.date(START_DATE, END_DATE));
s2Sr = s2Sr.filter(criteria).map(maskEdges);
s2Clouds = s2Clouds.filter(criteria);

// Join S2 SR with cloud probability dataset to add cloud mask.
var s2SrWithCloudMask = ee.Join.saveFirst('cloud_mask').apply({
  primary: s2Sr,
  secondary: s2Clouds,
  condition:
      ee.Filter.equals({leftField: 'system:index', rightField: 'system:index'})
});

var s2CloudMasked = ee.ImageCollection(
  s2SrWithCloudMask).map(maskClouds).median();
var rgbVis = {min: 0, max: 3000, bands: ['B4', 'B3', 'B2']};

Map.addLayer(s2CloudMasked,
  rgbVis,
  'S2 SR masked at ' + MAX_CLOUD_PROBABILITY + '%', true);

print(s2SrWithCloudMask);


texts: How to select and change only certain (same)values in attributetable in QGIS

texts:
I have an attributetable. Spread all over it are numbers like -99999997. I want to replace them by zero (0). There a many columns in which there are such numbers.
I have used the expression (an example with one column P-Ant_ARU)

case when   "P_ANT_ARU"  is  -99999997  then '0' else    "P_ANT_ARU"  end

for every feature in the column but I wonder if there is a faster way to use in the field calculator to get rid of all the -99999997 in one or perhaps two steps, because there are a lot of columns  :-)

I have also tried toint(0) but that doesn't work because you have to choose a field, and I want a value to be changed. The same problem occurs with the plugin QuickMultiAttributeEdit.

So how can  I select all those numbers in the whole table and  make them zero.


texts:
During a project on the spider fauna of Greece, we have to face the design of distribution maps of the species.

A good tool to work with is the design of Leaflet webmaps to be used on the website instead of static maps.

As the number of species is not small (1120 species = 1120 maps) am looking to find a way to create under a batch mode, webmaps with the Leaflet plugin.

All data are stored in a csv or shapefile with one column contain species name and rest the source (literature code) and coordinates in WGS84

Any idea on how can this work?


texts: Should I create multiple tables or one table for all PostGIS data?

texts:
My goal is to create an application in which I have 2 "types" of data that I will be querying against. I will have points and polygons. Ultimately I will be querying to find all of the polygons that a specific point exists in. I have downloaded the worlds admin level data, once in 6 separate shapefiles (one for each admin level), and in another instance, all in 1 shapefile (combining each admin level). I am also going to have to get different shapefiles from different sources besides GADM such as OSM, or census.gov.

Is it ideal to combine all of the polygon data into one table regardless of admin level, municipality, natural feature, etc, or to create a table for each potential type (adm1, adm1, locality, etc)?

My hunch is that combining them all makes the most sense because I wouldn't have to query a variable number of tables if I am looking for all of the polygons that a point lays within.

At the same time, I worry of such a potentially large table. Again, my goal is to load massive amounts of data into this table.

Has anyone done something similar to this, and can maybe shed some light on the problem?


texts: If/Then statement for calculating value based on the value of another item for ArcGIS Online

texts:
In the project I am working on some of my data needs to be doubled. I have created three main attributes. The first is pedvol, then there is intersection, and finally fixed_pedvol.

I want to take all cases of items where intersection = y and take the pedvol and double it and put that number into fixed_pedvol.

And ideas on how to do this on ArcGIS Online? I have been looking into IF/THEN statements and IF/ELSE but I cannot seem to figure this out.


texts: Is there a way to obtain polyline of pbf file?

texts:
I need to find boundaries of pbf file. Ideally - polyline, if this not possible - can be simple(min,max lat,lon rectangle). I download pbf from geofabrik(http://download.geofabrik.de/russia/northwestern-fed-district.html).

I try to obtain boundaries with

01:14:35 robots@HQ import → osmium fileinfo -eg data.bbox northwestern-fed-district-latest.osm.pbf
[======================================================================] 100%
(10.2919599,51.2626864,105.9294632,84.694631)

But this bbox seems to incorrect(Or I am doing something wrong).

However. I see on geofabrik web-page for downloading areas, on top right corner there is a map with an area of the pbf I am downloading and that is why I am pretty sure, that I somehow can extract this information from pbf.


texts: Displaying images as part of multiattribute labels using QGIS

texts:
I would like to display images referenced in an attribute column as part of a multiattribute label, so that I can eventually print them with the map.  I do not want to display geotagged images, since this allows me not to move the photos manually and also not to show multiple attributes at once.

Is there a way to style labels (by expression or html) as to also show the corresponding images of the individual features?


texts:
I would like know how can I make a thematic map (soybean production) with density points.

In this case the points are randomly plotted in a municipality area.

Please, see an example in the image:


texts: Getting 2 point features with max/min ID from each group with same value in two fields

texts:
I'm working in QGIS for a Digital Terrain Model.

I've got a table with point features, arguably vertices extracted from multipart iso-curves because the features are provided with a "vertex_par" value along with an "ID" (not unique, it recalls the original curve multipart feature ID).

Intent: I need to get the extremes of each original iso-curve lying on the property border into a table to feed my BIM-modeling app.

GOAL: I need to isolate all the points that are start and endpoints of the original curves (or parts).

Sketched in plain language:

For each group with the same "ID" value

    For each subgroup with same "vertex_par"

        I need to get the points with minimum/maximum "vertex_ind" value

QUESTION: How do I set an expression to filter these points?


texts: Update Point ID as Start Node and End Node to Line through FME

texts:
I need your help to perform the below process in FME.

I have two layers: Pole Data (Point) and Cable (Line).

Cable has no ID.

I want to keep the Cable which is going through the Pole.
Then make the cable as Single segment between two poles.
Then I want to update the Start ID and End ID (Pole ID) to each segment (Cable layer).

Please share the design of workbench and what parameters I will use in each transformer.


texts: How to auto rotate labels for Polygons in MapServer?

texts:
I want to auto rotate labels of polygons in Mapfile, all options that exist in MapServer refers to linear features only and not polygons. Is there any algorithm I can implement in the Mapfile to auto rotate label based on polygon rotation?

I need something similar to below image (taken from QGIS).


texts: Styling map elements based on their properties

texts:
I am currently generating tiles based on OSM data using Mapnik and Python and am quite happy using the XML styling but I am at a loss as to how to vary the styling on an element by element basis.

I would like to render OSM map tiles where certain elements are styled according to their attributes, for example, a roads colour is based on its length, or an area is coloured based on its size.

Is this possible?


texts: Determining percentage of intersection between two shapes In Oracle Spatial

texts:
I have to query polygons that intersect with one another.  Our polygons occasionally contain intersections that aren't intended (meaning just a few vertices cross boundaries).  I want to filter out these anomalies.

My current algorithm is to take the intersection of the two shapes and compare the area of the intersection with the area of the original polygon.  If the intersection area is 1/16th the original, I'm considering it valid for our business rules.

I have found sdo_intersection to be very slow (~30 seconds for 50 records).  I'm performing this operation against a national landgrid so the process runs for several days.

I've already fine-tuned my algorithm to exclude polygons with the 'inside' relationship and just work on the ones that I know have overlapbdyintersect.

Any ideas on how to make this faster.  At the end of the day, I want to know if there is any significant overlap between two polygons (not just a few vertices crossing over).


texts: Counting number of polygons within polygons using R?

texts:
I have two huge shapefiles. One contains polygons of one kind of vegetation and the other one representing a grid. A situation as follows but, of course, much bigger:

I need:

the efficient way to open and handle both files using R, because my computer "suffers" opening both using readOGR or shapefile even if I have a 32 GB RAM i7;
I need to get a matrix-like (or dataframe) containing the number of polygons for each square within the grid, something like this (from the previous example):


texts: Display map scale (like 1:20.000) next to ScaleBar

texts:
In my PageLayoutControl I add a ScaleBar. Additionaly I like to show the map scale/ratio (e.g. "1:20.000") next to it, but didn't find a way to achieve this with the ScaleBar. Am I missing something or is there another way to do this?

Here is, how I add the ScaleBar:

  IEnvelope envelope = new EnvelopeClass();
  envelope.PutCoords(5, 0.1, 12.4, 1); // Specify the location and size of the scalebar
  IUID uid = new UIDClass();
  uid.Value = "esriCarto.ScaleLine";

  // Create a Surround. Set the geometry of the MapSurroundFrame to give it a location
  // Activate it and add it to the PageLayout's graphics container
  IGraphicsContainer graphicsContainer = pageLayout as IGraphicsContainer;
  IActiveView activeView = pageLayout as IActiveView;
  IFrameElement frameElement = graphicsContainer.FindFrame(map);
  IMapFrame mapFrame = frameElement as IMapFrame;
  IMapSurroundFrame mapSurroundFrame = mapFrame.CreateSurroundFrame(uid as UID, null);
  IElement element = mapSurroundFrame as IElement;
  element.Geometry = envelope;
  element.Activate(activeView.ScreenDisplay);
  graphicsContainer.AddElement(element, 0);
  IMapSurround mapSurround = mapSurroundFrame.MapSurround;

  IScaleBar markerScaleBar = (IScaleBar)mapSurround;
  markerScaleBar.LabelPosition = esriVertPosEnum.esriBelow;
  markerScaleBar.UseMapSettings();


texts:
I have two fields in an attribute table and I want to add a third one which will contain the minimum of these two fields.

For instance,

Field 1  Field 2 Field 3

100      50      50

20       50      20


texts:
I need to automatically optimize text field character length to their current largest attribute character count. For example the field "MyName" is 250 characters long, but has attributes that are a maximum 6 characters long, the field "MyName" character length needs to updated from 250 characters to 6 characters.

I can do this processes manually (with possible errors), but I need a faster more reliable method to give to ArcGIS users. Before I start to write an ArcPy tool, has anymore written something like this before?

I'm guessing the ArcPy method would loop through each text field (A) and then loop through each row and note the largest attribute length (B), then it would pass the field name (A) and character length (B) to the output.


texts:
I am making a map on Auto Desk Infrastructure Map Guide Studio. Now I created a small ASP.NET website and hosted in local machine, Now on click of different buttons in flexible layout toolbar, I invoke URL for that asp.net page and it opens fine. I do some data querying on that page by selecting different options and execute query. The data changes in table in database and the effect is shown on the layer on map. Ok the issue is that I want to refresh the Auto desk map after invoke URL window is closed. And I want to do it by configuration in Map studio 2013. Does anyone have any idea how to refresh map automatically after invoke URL? Or in other sense if the data in the data source of layer is changed.

Addtional Info: the map always refresh when I click the refresh button from top or if I drag it or zoom it..Otherwise new layer data stays same.


texts: x and y coordinates will only display [x],[y] when field calculator is used to make a url

texts:
The standard formula for a google maps url is:

https://www.google.com/maps/@-15.623037,18.388672,8z

https://www.google.com/maps/@(latitude),(longitude),(zoom level)z

I am trying to make a new column on a point layer that would link you directly to google maps based on the layers x and y coordinates from other columns.

....but I am unable to get the values to show up it always fills it up with:

https://www.google.com/maps/@[x],[y],8z

Here is a sample that I am plugging into the field calculator.