texts:
I would like to georeference these fantastic old public domains maps
https://legacy.lib.utexas.edu/maps/ams/italy_50k/

I started from this one, and it has these projection parameters.

Then I have tried to georeference it, using the commands below and this proj parameters +proj=lcc +lat_0=45.90 +lon_0=14 +lat_1=45.90 +x_0=800000 +y_0=601000 +ellps=bessel +units=m +no_defs +k_0=0.998992911

# delete existent images
rm ./output_gcp.tif
rm ./output_warped.tif
# download the image
curl -L "http://legacy.lib.utexas.edu/maps/ams/italy_50k/txu-pclmaps-oclc-6540719-milano-45-iii.jpg" >./output.jpg
# add control points
gdal_translate -of GTiff -gcp 532.292 2433.53 405000 551000 -gcp 2618.57 2316.82 424000 551000 -gcp 2388.9 467.338 423000 567000 -gcp 537.666 569.428 407000 567000 -gcp 933.404 1364.4 410000 560000 -gcp 2090.5 1299.92 420000 560000 -gcp 1462.72 517.946 415000 567000 -gcp 1576.18 2376.68 415000 551000 ./output.jpg ./output_gcp.tif
# apply the proj parameters
gdalwarp -r near -order 1 -co COMPRESS=PACKBITS  -t_srs "+proj=lcc +lat_0=45.90 +lon_0=14 +lat_1=45.90 +x_0=800000 +y_0=601000 +ellps=bessel +units=m +no_defs +k_0=0.998992911" ./output_gcp.tif ./output_warped.tif

But I have a not so good result, a displacement of about 150 meters.

I think that the problem here is the definition of the standard parallel/s.

With this kind of projection (as I have read here), "two parallels and a scale of 1 are equivalent to
one parallel somewhere between them and a scale less than 1. Since there is a
scale, which is slightly less than 1, I think that there is one parallel"; it should be the Origin I can read in the projection info (45°54'N), and I have used it.

But it seems not work. Probably I must set two standard parallels, but which should i use?

EDIT

I have tried to use geographical coordinates of the 4 corners (using EPSG:1660):

# download the image
curl -L "http://legacy.lib.utexas.edu/maps/ams/italy_50k/txu-pclmaps-oclc-6540719-milano-45-iii.jpg" >./output.jpg
# add control points
gdal_translate -of GTiff -gcp 418.832 2486.39 -3.5 45.3333 -gcp 2689.8 2501.38 -3.25 45.3333 -gcp 2686.8 349.879 -3.25 45.5 -gcp 425.704 337.633 -3.5 45.5 ./output.jpg ./output_mm.tif
# apply the proj parameters
gdalwarp -r near -order 1 -co COMPRESS=PACKBITS  -t_srs "+proj=longlat +ellps=intl +towgs84=-104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68 +pm=rome +no_defs" ./output_mm.tif ./output_mm_warped.tif

It's a better result, but I would like to find a right way using the Lambert Projction.

Moreover it's wrong assumption because "The Monte Mario datum based on the Hayford (International 1909) ellipsoid was adopted in 1940, but this map explicitly says that it is projected from the Bessel ellispsoid" (see this reply).


texts: How can I write a zipped shapefile with ogr2ogr and /vsizip/?

texts:
After reading a bit of GDAL Virtual File Systems docs, I am trying to write a zipped Shapefile using this gdal command involving the use of /vsizip/:

ogr2ogr -f "ESRI Shapefile" /vsizip/test.zip\test.shp WFS:"https://www.wondermap.it/cgi-bin/qgis_mapserv.fcgi?&map=/home/ubuntu/qgis/projects/Demo_sci_WMS/demo_sci.qgs&SERVICE=WFS&VERSION=1.1.0&REQUEST=GetFeature&typename=domini_sciabili&bbox=544138,5098446,564138,5108446"

, and I get this error:

ERROR 1: Read-write random access not supported for /vsizip
ERROR 1: Cannot create test.shx while another file is being written in the .zip
ERROR 1: Failed to create file /vsizip/test.zip/test.shx: No such file or directory
ERROR 1: Terminating translation prematurely after failed
translation of layer domini_sciabili (use -skipfailures to skip errors)

The code is only writing test.shp (0KB though) inside test.zip.

I get what is happening (the error is self-explanatory), but sincerely, I don't know if I can go farther than this without probably using the GDAL API. Am I wrong?


texts: Edit data using arcobjects by SOE Web Services on Arcgis Server 10.0

texts:
I trying to develop WebService which main purpose is enable editing data. My Web Service get as param information about records, tables etc. Code looks like below:

IFeatureWorkspace ifw = null;
        IWorkspaceEdit iwep = null;
        try {
            ifw = new IFeatureWorkspaceProxy(new IVersionProxy(currentUserVersion));
            iwep = new IWorkspaceEditProxy(ifw);
            iwep.startEditing(false);
            iwep.startEditOperation();
            performOperation(mergeData, ifw, iwep);
//this method make store on records
        } catch (AutomationException e) {
            e.printStackTrace();
        } catch (IOException e) {
            e.printStackTrace();
        } catch (Exception e) {
            e.printStackTrace();
        } finally {
            try {
                iwep.stopEditOperation();
                iwep.stopEditing(true);
 } catch (AutomationException e) {

            } catch (IOException e) {

            }
        }

Everything works fine for one single request, but above code is horrible slow for multirequest for different version. Because different people has own version which is used to editing. I think the most time take startEditOperation and startEditing. Everything starts as capability of MXD file. Anyone can give me a solution how to using soe web services give capablities for editing data which will have a quite good performance?


texts: Creating buffer based on occurrence frequency using ArcGIS Desktop?

texts:
I am using ArcMap 10.4.1.

I have several point shapefiles containing different occurrence data for two different species of flower (from multiple sources, hence the multiple shapefiles). I want to create a buffer from a central point that includes a quantity of at least 50 each of each species. I'm familiar with the general buffer tool, so for example, I can create a buffer of a certain size and clip to it, but what I'm hoping to do is create a buffer that contains at least 50 occurrences of each species regardless of the size, from a certain point (i.e. a city).

I'm hoping this can actually be done in ArcMap? I've been searching Stack Exchange and elsewhere, and can't find this anywhere.

I've attached a screenshot of what the occurrence distribution looks like. Attribute tables consist simply of long/lat values. One of the species I am studying is protected, so I can't get more detailed than that. The primary study site is a natural area just outside of a major city, so we chose the city as an arbitrary "central point" to create the buffer from.

As of now, I have eight sources providing occurrence data, so eight different point shapefiles, that I created - so they should all be the same. I have an "Advanced" license, and the following extensions: 3D Analyst, ArcScan, Geostatistical Analyst, Network Analyst, Publisher, Schematics, Spatial Analyst, Tracking Analyst.


texts: Download data from OpenStreetMap that is not in English gives wrongs characters

texts: I have download data from OSM from Israel area but whenever I open the data in QGIS I am getting wrong characters and it doesn´t show Hebrew, even though on OSM it does show it in Hebrew. I know there is option to play with the Unicode but I don´t see this option in QGIS. Is there any way to change the Unicode of layers downloaded from OSM before opening?


texts:
I am trying to create a form of data driven page that incorporates linear referencing  very similar to the "Linear referencing scenarios" example:

I have a map that I will display on each page that shows a long polygon of a construction section for a pipeline. However, I just want the map to show the base imagery layer and the polygon.

The linear referenced data needs to be displayed outside the map area in an abstract manner much the same as the example, almost like a table underneath the map. If one were to take a straight edge and draw a vertical line on the page through the table, one would be able to determine each of the variables in question. (I keep thinking of the London Underground example, considering how the diagram shows a straight horizontal line where in reality it twists and turns)

I have eleven such routes (variables) to display for each page which shows about 8 km of the pipeline. (in the example, it would be accidents, lanes, material, speed, and quality). Seven of my routes are line referenced and 4 are point referenced. The line data are such things as political districts, or construction segments, etc. and the point data is road and stream crossings, etc.

I feel certain that someone must have paved the way in this area before. I am not the first person mapping multiple pages and multiple variables.

With data driven pages, in order to display some variable data on the page, you have to select one of the data driven page attributes. So I am thinking this route data has to be incorporated into my data driven page centers table.

I understand how routes work when drawn on a single map, but am not sure how to make the leap to data driven pages with linear referencing and dynamic segmentation, or how to abstract it outside the map, yet retain the linear horizontal reference of the map. (I have 60 pages to map) How do I incorporate these routes and their linear referenced data into a data driven page?

I am running ArcMap 10.2.


texts:
I would like to be able to export a layout in QGIS such that the resulting image is conic in shape, reflecting the CRS used.
In other words, something similar to what is show below:


texts: Adding shape of defined dimensions using PyQGIS in QGIS 3.4

texts:
I am attempting to add a feature with defined geometry as a layer. Having added the feature I checked on the dimensions of its sides and they measure differently.

This is the code that I am using.

x = [0,0,-5000,-5000,0,0]
y = [0,-1000,-1000,1000,1000,0]

vl = QgsVectorLayer("MultiLineString?crs=epsg:7405", "temp", "memory")
pr = vl.dataProvider()
pr.addAttributes([QgsField("name", QVariant.String)])
vl.updateFields()

PointList = []
for i in range(len(x)):
    p = QgsPointXY(x[i],y[i])
    PointList.append(p)

geom_line=QgsGeometry.fromPolylineXY(PointList)
f = QgsFeature()
f.setGeometry(geom_line)
pr.addFeature(f)
vl.updateExtents()
QgsProject.instance().addMapLayer(vl)

The lines do not measure as 2000m and 5000m, but as 1997m and 4992m (although opposite sides do not measure the same).

I wonder if it is to do with the CRS or the measure tool? I think that I am missing something fundamental here.


texts: Trigger a Python script to download external data for ArcGIS web application

texts:
I would like to create a very simple web application with ArcGIS. The web application would ask the user to provide the coordinates of a point and as a result the web application would provide  all the POIs (points of interests like parks, restaurants, places of worship, train stations etc) within 2 km from the provided point.

I have a Python script which is able to retrieve data from Google Maps under those circumstances and I would like to use it. The question is:
- How can I trigger this Python script when a user inputs the coordinates to download the data and show it in the ArcGIS application?


texts: Downloading MODIS products with R package rts

texts:
I'm having trouble batch downloading MODIS products with the rts package in R.
I'm using the code below to download the MODIS products:

library(rts)
library(raster)
library(RCurl)

bands = "0,1,0,0,1,0,0,0,0,0,0,0"

ModisDownload(x='MOD11A2',h=c(24,25),v=c(5,6),version='006',dates=c('2011.01.01','2011.03.31'),mosaic=F,bands_subset = bands, proj=F)

When I run the code, I get this error message:

Error in .._MD_curlHandle <<- RCurl::getCurlHandle() : cannot change value of locked binding for '.._MD_curlHandle'

I'm not sure what this error message means.


texts:
I have a .las file and want to show the outline of the big rocks located in this area.

So far I have done the following:

Obtain a ground model from the .las file including the rocks.
Calculate the slopes in the ground model.
Only show slopes greater than 30 degrees.

Step 2 and 3. are fine. My problem is to obtain a reasonable
ground model in step 1. (Step 3 unfortunately removes some of the smaller rocks, but this is not very important).

The .las file is downloaded from hoydedata.no. It is available here. In the linked map there are two big rocks showing up in the DSM Hillshade layer. One rock next to the road, and one next to a small lake. The rock next to the lake does not show up in the DTM Slope degrees layer.

Further information:
The points in the .las files are classified as 1, 2, 7. It is somewhat arbitrary if big rocks are assigned class 1 or 2.

A rock is big if it is greater than 2 meters in at least 2 directions (and usually in all 3 directions). However, a method which can outline smaller rocks as well would be great.

The rocks are typically located in forests, swampland or mountains. There are usually no buildings around. There might be cliffs around, but these I want to show as well. Hence, an alternative could be to identify and remove the vegetation.

Are there any tools to remove vegetation points from .las files?

I have tried various tools for step 1:
pdal works when I increase the slope parameter slightly, but I worry that there are some trees showing up:

pdal ground -i /data/in.las -o /data/out.las --slope 0.7 --denoise

lidR has an impressive tree segmentation algorithm which color the rocks correctly, see 2. However, I was unable to extract the rocks or remove the vegetation. It looks like the rocks are classified as trees.

lastools works reasonably well with the following parametres (suggested by the data providers):

lasground.exe -i *.laz -ultra_fine -ignore_class 7 -step 5  -spike 1 -offset 2.5  -odir ground -olaz -cores 5
las2dem.exe -i ground\*.laz -slope -keep_class 2 -step 0.25 -opng -odir dtmslope

Unfortunately I do not have access to a licensed version of lastools, and I only have limited access to a Windows environment for running lastools.

Preferably I would like to use open source tools.

I am completely new to GIS and LiDAR.

How do I tweak the parameters or should I use different tooling for step 1?

For step 2 I currently use lidR and the following code:

require(lidR)
las <- readLAS("out.las", select="xyzc")
las <- lasfilter(las, Classification == 2)
gnd <- grid_terrain(las, res = 1, tin())
sl <- terrain(gnd, opt="slope")
png("tmp.png", width = 1920, height=1080)
plot(sl)
dev.off()

This produces the picture below from the output of pdal. In the picture you can see the outline of two rocks at the bottom center next to the road and the small lake.

.


texts: Buffers not showing on map but viewable in Table of Contents of ArcMap?

texts:
I am completing Assignment 9-2 in GIS Tutorial ArcGIS 10.1 Book. I am suppose to set the map to a UTM for California (either 10N or 11N). I am then suppose to add a California Counties layer and Earthquakes layer (which has the same projection as my map and was added from my File Geodatabase). Then I set up a definition query in my Earthquakes layer to have only the earthquakes with an MMI of 7 AND after the year of 1985.

I am suppose to have a 20 mile buffer for each earthquake. I went to Geoprocessing -> Buffer -> Selected my Earthquakes layer -> Saved the output in my File Geodatabase -> Linear Unit as 20 miles -> Dissolve Unit: None. The buffer was added to my Table of Contents, but nothing showed up on my map. I made sure that the buffer had the NAD 1983 projection which was what the map asked for. I also clicked "Zoom to Layer," and only saw small dots with lines coming out of them, but no buffer on my Earthquakes. I also tried to go to export the layer into a UTM projection. However, even with a UTM projection, the buffer layer wouldn't show on the map.

I am not sure why I am getting this error.


texts: Why is "Replace With Model" not replacing my multipatch with a SketchUp model?

texts:
I am following the steps recommended by ESRI:

extrude polyline in ArcScene
Save as Layer File
convert Layer file to multipatch
convert multipatch to Collada
import collada file to Sketchup
edit model in SketchUp
start editing in ArcScene
select multipatch with Edit Placement Tool
Replace with Model - Navigate to saved skp (NOTHING HAPPENS)

This only seems to (sort of) work if I export from SketchUp to Collada but only fully closed faces (and only some faces) are visible in ArcScene. Some faces can be seen from the reverse side. It's basically a mess. Not what I was expecting.

Why can't I replace the SketchUp model itself? Is this procedure only meant for very simple models?

I initially wanted to build the model in SketchUp but it seems impossible to export anything out of ArcGIS without the SketchUp 6 Plugin for ArcGIS.  Exporting complex surfaces to COLLADA almost always fails and since TIN cannot be exported to SketchUp directly, it has to first be converted to Multipatch. This either fails if the surface is not very simple or yields unusable results. I then decided to do it the other way around (build the model in ArcScene and add SketchUp models to it) but I don't seem to be able to replace simple multipatches with models edited in Sketchup.

I was able to create large Sketchup models in the past (up to version 9.3.1) with complex TIN, polygons, points and lines exported directly from ArcGIS using the SketchUp 6 plugin. I am not really sure what the workflow should now be if the recommended ESRI workflow doesn't really work.


texts: Cut a multispectral TIFF into smaller parts

texts:
I have a multispectral image with 6 bands that I want to cut into smaller parts. It seems like a trivial task but I can't make it work.
Approach 1: using rasterio.windows. I choose how many tiles I want, calculate how many pixels the individual tiles have and then read the corresponding part of the original image and write it to a new file:
with rasterio.open('path/raster.tif') as file:
    src = file.read()

bands, height, width = src.shape

grid_size = 5
grid_width = int(width/grid_size)
grid_height = int(height/grid_size)

with rasterio.open('path/raster.tif') as file:
    profile = file.profile
    profile.update(
    width=grid_width,
    height=grid_height)

    for i in range(grid_size):
        for j in range(grid_size):
            window = int(width*i/grid_size), int(height*j/grid_size), grid_width, grid_height

            src = file.read(window=Window(window[0], window[1], window[2], window[3]))

            with rasterio.open('data/cut_{}_{}.tif'.format(i,j), 'w', **profile) as sink:
                sink.write(src)

Problem: This correctly cuts the raster into pieces (in this code it is 5*5 = 25 tiles) but I don't know how to correctly set their coordinates (the "transform" parameter in profile). In my code, all the tiles are placed in the upper left corner of the original TIFF which makes sense because I use the original TIFF's profile for writing the tiles.
Approach 2: using gdal2tiles
import gdal2tiles

gdal2tiles.generate_tiles('data/clipped_data/kuopio.tif', 'data/')

Problem: I get an error: ERROR 6: PNG driver doesn't support 6 bands . So it seems this is made only for standard 3 band images. I haven't found an option to change that.
Question: How do I correctly cut a multispectral raster into a number of smaller rasters? I am using Python so this answer is not helpful for me.


texts: Placename matching algorithms and services for geocoding

texts:
Are there any placename-specific matching algorithms, routines or services available?  By placename specific I mean not just a simple string-difference match, or soundex/metaphone, but something that understands geographic abbreviations and conventions.

For example, match "W. Mt Washington Pkwy" with "West Mount Washington Parkway"; or "Oxon" with "Oxfordshire"; or "Oxford-st" with "Oxford Street"; or "N.S.W." with "New South Wales"; or "1234 Main St NW" with "1234 Apartment 345, North West Main Street", or "St. John's Ferry" with "Saint John Ferry".

Most of the addresses I'm dealing with (to geocode, to plot on maps for historical GIS) are at the town level. In most cases there are further elements of the address available (country, state, province, country) to avoid too many false positives.  I do have data to match against, the challenge is to get the names to match, allowing for varied spellings and abbreviations.

There are several stages to match, assuming we're in the same state/province, including but not limited to:

simple string matches, allowing for differing word order
expand basic geographic abbreviations (W, St, Cty, and so on), based on the language(s) and conventions of the country, really just a bunch of look-up tables
expand local admin area abbreviations and codes (NY, NSW, Oxon, and so on)
look up placename changes and language alternatives

This is all basic stuff (the later three depend more on data than on an algorithm).  Although some parts are tricky, like working out when to expand St as Saint, or Street (or try both).

So I assume this widely needed problem has been solved many times before (by people much more expert than me!), so I shouldn't re-invent the wheel.

But, even limiting it to English language countries (UK, Canada and so on, as well as the US), I haven't been able to find any published work or software for this.  I must be looking in the wrong places.  Does anybody have any suggestions for existing work I can build on?


texts: Adding multiple datasets to ArcMap with the same coordinate reference system but different transformations

texts:
I want to be able to load multiple datasets to my mxd, using various CRS (coordinate reference systems).  Some of these datasets will use the same CRS, but with different transformations.  When I load the first dataset to the project for a given CRS, I am prompted to specify the transformation.  However, when loading subsequent datasets on the same CRS - the option to specify a different transformation is not available (it presumably uses the previous selection by default).

Is there any way to specify the CRS and transformation for every dataset loaded to my project (i.e. every time I load any data), and effectively disable this default setting?  I do not want to have to use the project tool from the toolbox to effectively duplicate all of my data each time I wish to do this.


texts: How to use the ISpatialFilter interface with ArcGIS 10.1?

texts:
From what I've heard till now:

1) ArcGIS from version 10.1 and forth will only support http connections: local connections will not be supported any longer.

2) The code

ESRI.ArcGIS.Carto.IMapServer iMapServer = (ESRI.ArcGIS.Carto.IMapServer)ctx.ServerObject;

or the code

ESRI.ArcGIS.ADF.Web.DataSources.ArcGISServer.MapResourceLocal mfunct = (ESRI.ArcGIS.ADF.Web.DataSources.ArcGISServer.MapResourceLocal)gisResource;
                ESRI.ArcGIS.Carto.IMapServer iMapServer = (ESRI.ArcGIS.Carto.IMapServer)mfunct.MapServer;

can't both work with an internet service connection because they require a local connection.

3) The ESRI.ArcGIS.ADF.Web.SpatialFilter object cannot be set to identify what kind of query I'd like to perform. I mean, with the ESRI.ArcGIS.Geodatabase.ISpatialFilter interface you can do something like

spatialFilter.SpatialRel = ESRI.ArcGIS.Geodatabase.esriSpatialRelEnum.esriSpatialRelIntersects;

to specify what you want as result. SpatialFilter hasn't this feature. I need the ISpatialFilter spatial relationship feature.

4) If you put it all together, it figures out that you need an IMapServer-implementing object to perform the ISpatialFilter query

ESRI.ArcGIS.Geodatabase.IFIDSet fidSet = iMapServer.QueryFeatureIDs(mapDescription.Name, layerID, queryFilter);

So if you read back points 1 and 2... what am I supposed to do?

I've also asked this on the official forums but still no answers.

Please help because I'm getting crazy on this.


texts: How to clip a MID by a shape using python and gdal/ogr

texts:
For a project i have to clip point data contain in a mid file with a shape.
Mid file have no header, it's like a text file. Values are separated with comma, and field_2 and field_3 are wgs84 coordinates.
I know that I maybe can do it with gdal/ogr but I'm really lost after reading many documentation.
Somebody can give a way ?


texts: How to avoid the CRS Dialog when creating a vector layer memory?

texts:
I am trying to create a vector memory in a plugin. Basing on some examples,I was able to create a memory vector layer and to add it in QGis with : QgsMapLayer.instance().addMapLayer(myLayer,True)

But when the layer is added , QGis always displays the crs dialog to specify the crs of the layers

:

I tried to specify myself the crs with :
layer.setCrs(QgsCoordinateReference(4326)) but it still shows me the crs dialog.

Do you have any ideas to how avoid this dialog ?


texts: Storing feature layer content inside ArcPy script tool?

texts:
I’m developing a script tool using ArcPy. I am wondering if there any way to store layer contents in my script in an encoded way and programmatically create this layer instead of referencing external shp/geodatabase feature class?

I need this to distribute only the .tbx file without any dependent .shp files.
The nature of layer – rectangular grid of polygons with attributes, inside a relatively complex boundary(Iowa state border), so it difficult to create this grid programmatically without at least border geometry.
Now I’m using code like:

foo="C:\foo\bar.shp"
arcpy.MakeFeatureLayer_management(foo, "bar")

to add my layer. Would it be possible to somehow encode this feature layer geometry and attributes, store it as variable/list (like WKT, GeoJSON, KML) and use something like that

jsonEncodedString=" {\"type\": \"FeatureCollection\",\"features\": [{\"type\": \"Feature\",\"properties\": {},\"geometry\": {\"type\": \"Polygon\",\"coordinates\": [[[87.890625,67.60922060496382],[78.046875,64.47279382008166],[82.265625,63.54855223203644],[88.9453125,64.16810689799152],[90.3515625,65.80277639340238],[87.890625,67.60922060496382]]]  }},{\"type\": \"Feature\",\"properties\": {},\"geometry\": {\"type\": \"Polygon\",\"coordinates\":[[63.6328125,63.860035895395306],[55.54687499999999,64.16810689799152],[54.4921875,63.23362741232569],[59.765625,60.58696734225869],[63.6328125,63.860035895395306]]]}}]} "
some_function(jsonEncodedString, "layer_name")

to create a feature layer with attributes instead of referencing an external shapefile/feature class?


texts: Selecting only attributes specified using ArcPy?

texts:
This is a working code made using exporting the script from ModelBuilder.

import arcpy

# Local variables:
Links_Final = "Links_Final"
condition = "\"RR1\"= 777 OR \"RR2\"=777 OR \"RR3\"=777 OR \"RR4\"=777"

# Process: Select Layer By Attribute
arcpy.SelectLayerByAttribute_management(Links_Final, "NEW_SELECTION", condition)

However, the I want the value 777 in the condition to be taken from the user. So I added this code.

value = int(arcpy.GetParameterAsText(0))
condition = "\"RR1\"= value OR \"RR2\"= value OR \"RR3\"= value OR \"RR4\"= value"

Now the code throws an error saying, Invalid expression
Failed to execute (SelectLayerByAttribute)

I have also followed the thread https://stackoverflow.com/questions/32483187/trying-to-create-a-simple-select-by-attribute-script-based-on-an-input but that does not work either.

Note: My data type on the field is long integer and the data type taken from the user (as specified while creating the script) was long (as I could not find int type there)


texts: Export raster file geodatabase to jp2 or jpeg?

texts:
Can I export or copy a raster from a file Geodatabase and create a output to jpeg or jp2000 formate? I was reading two different helpline on ArcGIS 10.1 and wasn't sure if I was supposed to export jp2000 from another file geodatabase to file geodatabase.

I tried to do that the other day but it failed. Don't know if it was the space was issued or maybe I wasn't doing the proper way.

Thanks.


texts: QGIS Azimuthal Orthographic Projection with Raster

texts: I want to use QGIS to make a world map that has two or three azimuthal orthographic maps.  I would also like to use Natural Earth's Cross-blended Hypsometric Tints.  Is this possible? The example was done using ArcGIS Pro.  I would like to use QGIS if possible.


texts: Drawing the right coordinates in my grid in QGIS

texts:
I'm trying to put the numbers of system of coordinate in my map using WGS84/18S. All my shapes are in that projection and I'm using units of the map around 1500 point in "x" and "y" but, I think those aren't the right coordinates.
When I use just WGS 84, the numbers are correct. Maybe with general projection it works, but with a specific one it doesn't work, I don't know.


texts: Getting wrong values trying to convert to a latitude longitude system

texts:
By using listgeo I got the following information:

PCS = 31468 (name unknown)

If I am not mistaken this translates to DHDN / 3-degree Gauss-Kruger zone 4, and by using this site I got a proj4 string:

+proj=tmerc +lat_0=0 +lon_0=12 +k=1 +x_0=4500000 +y_0=0 +ellps=bessel +datum=potsdam +units=m +no_defs

Trying to convert values drom this projection system to latlong using proj4, I am doing the following:

double x1, y1;
x1 = 4386596.4101;
y1 = 5237914.5325;

pj_merc = pj_init_plus("+proj=tmerc +lat_0=0 +lon_0=12 +k=1 +x_0=4500000 +y_0=0 +ellps=bessel +datum=potsdam +units=m +no_defs")

pj_latlong = pj_init_plus("+proj=latlong +datum=WGS84")

pj_transform(pj_merc, pj_latlong, 1, 1, &x1, &y1, NULL )

What I get back though is 0.226465 0.917617.

What am I doing wrong here?


texts: non-editable date / UUID is not saved in attribute table

texts:
I'm using QGis 2.14.12 Essen on Windows 7 Professional.

I want to set up a project that includes a feature class which has a non-editable° creation_date  field and a non-editable UUID field.  I have found a way to do it but it doesn't work the way I expect it to. I was able to recreate the problem on a colleague's computer.

Here's the problem:

I created a new .shp. In the layer settings, I've created a field "creation_date" (field type "date") which has the element setting° "date/time" and a UUID field (field type "Text/string", length 50) which is "UUID Generator"°.  If I tick the box "editable"° in each dialogue, everything works just fine when creating new features.

If I untick the box "editable", the pop-up upon creation of a feature displays the correct date and it's not editable. Same for the UUID.
Upon saving, neither of the entries show up in the attribute table. I can see them via the identify tool, though.

Can you reproduce this problem, is this a bug or am I making some mistake?

(I'm using the German version, so I'm not sure if I translated the words (marked with °) correctly. I hope you know which buttons/settings I'm talking about. Please help me out here to make the question more understandable.)

EDIT
If I open the .dbf in Excel, the editable columns has data, but not the non-editable.
Removing layer from Project and adding it again, does not Show the Information in the Attribute table and it stops the Information from being available in the identify tool.


texts: ST_MakePoint returns an error when using Python

texts:
I am trying to use ST_Point or ST_MakePoint to create points from a table of latitude and longitudes coordinates. However I am consistently getting an error even when trying the examples.

For example

pd.read_sql("SELECT ST_SetSRID(ST_Point(-71.1043443253471, 42.3150676015829),4326)", my_db_connection)

returns an error

DatabaseError: Execution failed on sql '
SELECT ST_MakePoint(-71.1043443253471, 42.3150676015829)
': function st_makepoint(numeric, numeric) does not exist
LINE 2: SELECT ST_MakePoint(-71.1043443253471, 42.3150676015829)
               ^
HINT:  No function matches the given name and argument types. You might need to add explicit type casts.

It seems that this has to do with how the query is being passed from Python (I've tried using pandas and psycopg2 but no luck) but I am not sure why this would matter when I am trying to run my query using real data (SELECT ST_Point(lon, lat) FROM my_table).

Casting doesn't seem to help either. Has anyone come across this before?


texts: How can I add 3 axis orientation information to my photos?

texts:
I have a survey project where I need to record not only the location, but the direction (bearing and angle of elevation/depression) that the camera was pointing at the time the picture was taken.  Accuracy is not critical.  +/- 10 degrees would do.

I would prefer not to have to buy a new camera.

I don't need another GPS. I need to use a handheld GPS to get to each photo spot anyway, and integrating the GPS info with the photos is an easy ExifTool script.

I don't want external cables. This is an all weather project in bushy terrain; cables catch on things, and require open port covers. Similarly I don't want 'big lumpy things' like FotoMapr. (Which also has GPS that I don't need.)

The ideal solution would be a device the size of a pair of stacked nickels that would attach to the camera's hot shoe. The device would have 3 axis so that I could get compass direction and elevation angle. It would record these along with a time stamp at the time the picture was taken. Calibration to turn raw data to true north, true elevation would happen at data merge time. The device would have a mini-USB port for calibration and data retrieval.

If it requires power, I'd like it to get it from the hot shoe, but battery is acceptable. The unit should have a cost under $100

I originally asked this over at the Photography exchange.  They suggested I try here.

Rebuttal to Answer # 1

(Comments only allow 500 characters and no paragraphs.)

It's a bad solution:

It does not record dip, but just uses a 3 axis magnetometer to correct for tilted camera angles.
It has cables, and requires open ports on the camera.
It sticks out.  It will catch on brush.
It's a GPS.  I don't need to pay for another GPS.

The 3 axis compass is a single chip under 1 cm square.  It retails on a 1 inch printed circuit board for $10 quantity 1.  It should be possible to make one that would in essence be a bump on the flash shoe.

If it can't get a time sync from the camera, it would also need to have a clock chip.  If it has a USB port, then it gets clock reset whenever attached to the computer.  The computer could also use the amount of drift to apply a linear correction to the time stamp.  So now our chip count is up to 4:  Memory, USB controller, 3 axis magetometer, clock.


texts: Using multiple attributes feature for symbology in ArcMap?

texts:
I am trying to make a map using US census data showing vacant housing to total housing on the same map.

I came across the "multiple attributes" feature under the symbology tab and wondered if this is the way to go and, if so, how to use it?

I have only been able to get one type of symbol, the gradient feature to appear.

I'm using ArcGIS Desktop 10.


texts: GeoJSON or TopoJSON outline of Milky Way

texts:
I'm looking for a GeoJSON or TopoJSON file with the Milky Way outline on the Sky sphere, for an astronomical almanac web application I'm building with the D3js JavaScript library.

I have searched in vain with google...

The best would be if it has a few layers roughly indicating the brightness levels (brighter in the middle of the disk, and towards the MW centre).

Alternatively, does anybody have an idea how I could create such a file?


texts: How to reference piped input in ogr2ogr SQL statement

texts:
I have a geojson data set that I am piping into an an ogr2ogr command using /vsistdin/. However, I am having trouble figuring out how to actually refer to this input in the SQL statement I'm running. The below simplified example has been erroring out when I run it (I'm doing some previous manipulation/cleanup on the data, not just using cat, that I have left out for brevity's sake). It's obvious what the error is, I just can't figure out what else to refer to the piped input as.

Example:

cat input.json | ogr2ogr -f geojson test.json /vsistdin/ -dialect sqlite -sql "SELECT *
FROM /vsistdin/"

Error:

ERROR 1: In ExecuteSQL(): sqlite3_prepare_v2(SELECT *
FROM '/vsistdin/'): no such table: /vsistdin/


texts: Converting tax map polygons from Shapefile to table of map number and corner coordinates

texts:
Fairfax County, Virginia provides a Shapefile here which shows the boundaries of the tax maps as polygons. Each polygon simply consists of four corners, and what I'd like to do is create a table showing the map number and the coordinates of the corners, such as the following:

Map number, bottom-left, bottom-right, top-right, top-left
77-1, -77 32, -76 32, -76 33, -77 33

(though they aren't quite square like this)

I don't have access to ArcGIS software, but do have QGIS and OpenJUMP, and can install any other software that is freely available. I believe this should be very easy but I'm just missing something. Thanks!


texts: QGIS: vector union produces a single feature only, same as the input layer

texts:
I have two layers in QGIS 3.16,

a MultiPolygon layer with a single feature, the area of Papua New Guinea
a MultiLineString layer of OpenStreetMap roads, clipped to the extent of the first layer

Using the Union tool, with the first layer as input and the second as overlay, I was hoping to obtain a new layer with the roads plus all non-road areas - but all I get as output is a single feature the size and shape of the first layer. Can anyone advise?
I tried using Dissolve to merge all the roads into one vector before doing the Union, but that didn't work either. A validity check throws up no problems.
I am following an online tutorial in least cost path analysis called "Evaluating Landscape Permeability in Quantum" - it tells me "A union must be completed on all of the vector layers that do not have full data coverage of the boundary file. The union will act to fill any NULL values in the dataset to be 0s instead; this will become important when using the Raster Calculator in later steps". Eventually, I want to use GRASS and a DEM to calculate a cost distance raster to the nearest roads.


texts: Varying Line Widths According to Attribute Values

texts:
I'm trying to draw a vector map of Europe on QGIS 1.7.3 with lines symbolizing the transmission of energy (power transmission lines) between the different countries on the map. I am trying to correlate the transmitted power in gigawatts (that I have entered as an attribute for each line) with the thickness of that line (instead of a graduated or categorized color ramp). Can you please tell me if there's a way to do that?

Thank you.


texts: QGIS Layout - Add image in attribute table

texts:
I am creating a layout of map using QGIS 3.10
I want to show attribute table with each point logo in table.
I have a field which has the complete path to each logo and this field is already being used in symbology to show each respective college logo.e.g. in below picture attribute table is shown and in pic_path column instead of actual path i want logo to be shown.


texts: SDO_GEOM materialized view: Why is static field slower than dynamic field?

texts:
EGDB: Oracle 18c; ArcGIS 10.7.1
External system: Oracle 12c

I have a materialized view in an EGDB that syncs data from an external system (on a schedule via a dblink).
The materialized view has an SDO_GEOM column (since using SDE.ST_GEOMETRY isn't possible due to an Oracle limitation here and here).
create materialized view my_gis_schema.wo_mv --EGDB
build immediate
refresh force
on demand
start with sysdate next sysdate + (5/(60*60*24))
as
select
    cast(workorderid as number(38,0)) as objectid,
    wonum,
    status,
    --other fields
    longitudex,
    latitudey,
    sdo_geometry(2001, 80085, sdo_point_type(longitudex, latitudey, null), null, null) as shape
from
    maximo.workorder@gistomax --dblink

create materialized view log on maximo.workorder with rowid;  --external system

Question:
When I add the MV to ArcMap, the layer works fine.
However, the performance is surprisingly poor. It takes 10 seconds to load 6,000 points in the map.
In contrast, if I ignore the static SDO_GEOM column, and instead construct a shape column on-the-fly from the non-spatial XY columns, the performance is much better (1-second):

I generated an SDO_GEOM column on-the-fly in a regular view (fast!)
And alternatively, I created an XY event layer from the non-spatial XY columns in the MV (also fast)

Why is the static SDO_GEOM column so much slower than columns that are calculated on-the-fly?
This seems very counterintuitive to me.


texts: Retrieving LiDAR metrics from relative vertical layers of point cloud?

texts:
I need to obtain the proportion of all returns within 5 equally spaced height bins between the cutoff (1 meter) and the maximum height value. I need to get this proportion of all returns for 1000 cylinders with diameter of 28 m (I have the center coordinates of each cylinder).

I am planning to use FUSION to do this task. It seems that PolyClipData will work to clip the cylinders if I generate a shapefile file buffering my points. Then, CloudMetrics (strata:[#,#,#,...) would compute the returns in various height strata. However, the example in the user's manual shows that strata will be fixed (from the default strata) i.e. 0.15, 1.37, 5, 10, 20, 30.

What I want is to have 5 equally spaced height bins between the cutoff and the maximum height value, so the height of the strata will vary according to the highest return. How can specify this function in the switch strata?


texts: Creating a point in QGIS 3.2.2 with coordinates?

texts: I am trying to plot a point using specific co-ordinates in QGIS 3.2.2 and have been told to use 'numerical digitize' or similar to do this. I can't find a plugin that QGIS 3.2.2 will accept. Can anyone help?


texts: What mathematical processes does a GPS receiver use to determine position?

texts: I'm writing a report on the mathematics of GPS, but I'm wondering what methods the GPS actually uses. What I've done is set up the four sphere equations and solve them using the Newton-Raphson method for four variables (using the Jacobian matrix of partial derivatives). But is this what a GPS receiver actually does, with the an initial guess and iterations? Because I saw that on several websites, but I've also read about other methods such as a least squares solution. Or is this something that depends on the receiver?


texts:
I have time data about transactions by month. I'd like a Torque visualisation with a month-by-month slider. The raw data looks like this:

month,lat,lng,transactions
January 2014,52.2,0.0,7
January 2014,52.7,0.0,10
February 2014,52.2,0.0,10
February 2014,52.7,0.0,20

I'm not sure how to do this. CartoDB's "date" columns want dates which specify a day (not unreasonably), so it becomes more like this:

01/01/2014,52.7,0.0,20
01/02/2014,52.7,0.0,20

Then when I use the Torque visualisation, I get a brief flash on these days, but the map is mostly black. I guess I could get around this by adding a row for every single day, but that would make my data file unmanageably large (it's already over 100MB).

The other way I could get around this is to specify the steps in the slider as always being the first of the month. However, this doesn't seem to be possible right now.

How can I show data by month?


texts: Inserting null value into geometry column using PostGIS

texts:
I'm trying to insert a null value into a geometry column, which I cannot do:

  ERROR: parse error - invalid geometry

So what I'm trying to accomplish is to erase a point from a map and save that information - the deletion of the point - to the database.

I am new to spatial data and all I can find from the internet is how to insert points, polygons and whatnot, but not how to remove that information.

The table where the geometry column resides has other information as well, so removing a row is not something I can do.


texts:
I am trying to run QGIS 3.6 on a Mac with the latest OS.

I have the following gdal error message:

>>>>>>> gdalwarp -of GTiff -tr 250.00115934006 -250.00168180475995 -tap
>>>>>>> -cutline "/Users/Laura/Documents/2019 Adjusted Maps/CA Base Map
>>>>>>> Dissolved/CACTYBASE7dis.shp" -crop_to_cutline
>> "/Users/Laura/Documents/2019
>>>>>>> Adjusted Maps/NLCD adjusted/NLCDAlbers.tif"
>> "/Users/Laura/Documents/2019
>>>>>>> Adjusted Maps/NLCD adjusted/NLCDClipCA.tif"
>>>>>>>
>>>>>>> GDAL command output:
>>>>>>>
>>>>>>> /bin/sh: gdalwarp: command not found
>>>>>>>
>>>>>>> Execution completed in 0.19 seconds
>>>>>>>
>>>>>>> Loading resulting layers
>>>>>>>
>>>>>>> The following layers were not correctly
>>>>>>> generated.

/Users/Laura/Documents/2019 Adjusted Maps/NLCD
>>>>>>> adjusted/NLCDClipCA.tif

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

I have applied all the steps found on Kyngchaos.com, but it still isn't working

I think there may be a problem with this pathway but I am unsure which part needs to be fixed. Any suggestions?

/Library/Frameworks/Python.framework/Versions/3.6/bin:/Library/
Frameworks/GDAL.framework/Versions/2.4/Programs:


texts: Change symbol types for point features on mapview in R

texts:
Using the R mapview package, I'd like to use different symbol types/shapes, linked to a column in the table. e.g.

library(mapview)
library(sf)

head(breweries)
set.seed(999)
breweries$coop <- factor(replicate(nrow(breweries), sample(c("yes", "no"), 1, replace = T)))
table(breweries$coop)

mapview(breweries, zcol = "coop")

e.g. a functionality like this:

plot(breweries[,"founded"], pch = c(1,2)[breweries$coop])

Are there different symbols available for point features? If not, how can I use custom symbols but maintain a functionality that allows graduated colours, as in the default circles.


texts: Georeferencing Control Points Aren't Saving

texts:
I have some CAD data that needs to be georeferenced. I create my first control point, and everything goes well, the drawing shifts to the correct control point.  As soon as I try to create a second control point, the first control point disappears.  Then it acts as if my second control point is now the first control point.

I've done this before, but I'm not sure what the problem is now.

Process:

I bring in my .dwg file.
Open the georeferencing toolbar, select the file I want to georeference, then I click the "Add Control Points" button.
I click where I want my first control point, then right click and select "Input DMS of lon and lat." (I prefer to create control point using this method).
Once I enter the lat and long, that creates my first control point.
I repeat this process to create my other control points, however instead of creating a second and third control point, it "starts over" at control point #1 and will NOT create a second and/or third point.


texts: Can I convert an IFeature to ILayer and add it to the map?

texts: I want to use FeatureClass.Search() to get a specific Feature and I want to convert this one to ILayer and set it to the Map, is there any way to do that?


texts: I can't find CRS EPSG 27700 in version 3.6.0

texts:
I have tried to use the filter element in the Project properties menu to find EPSG 27700 with no success.

Any suggestions please for converting to British National Grid


texts: Convert geometry set to simple features collection using `sf` in R

texts:
How can geometry set objects be converted to simple features collection objects  with the package sf in R?  I want to add attribute data to the geometry but am prevented from doing so because geometry set excludes including data.

Converting simple features collection to geometry set is done using st_geometry, but what is the equivalent in reverse? See sample code below:

#load sample data set
nc <- st_read(system.file("shape/nc.shp", package="sf"))
#Simple feature collection with 100 features and 14 fields
#geometry type:  MULTIPOLYGON
#etc

#convert to geometry collection
nc_geom<-st_geometry(nc)
#Geometry set for 100 features
#geometry type:  MULTIPOLYGON
#etc


texts: Trouble implementing Key authentication in geoserver

texts:
I'm trying to implement key authentication in geoserver 2.3.3. The steps I followed are:

Added the authkey plugin
Created the authkey filter
Secured the layers by assigning workspaces to different users
Added authkey filter to the filter chain

I'm passing the authkey(obtained from authkey.properties file) as a parameter along with the request while accessing the layers from the client side of my application. But still the authentication does'nt work and the layer's are inaccessible. Is there anything that I'm missing out or anything wrong with my procedure.Kindly help me with this.

Thanks.


texts: Identify what projection this Australian map uses?

texts:
Im looking to put this map into google maps as custom tiles, but i now learn what projection is and that this map does not use the same format as google's default.

You can create custom projections i read in the API.

I suspect this map uses a common projection, likely the same as this one:

I believe it will likely be this one:

e.g.
National Geoscience Datasets, Geoscience Australia
and
Conic Projection – Lambert Conformal Conic


texts: Calculation of Standard Vegetation Index from EVI for Image Collection in Google Earth Engine

texts:
I'm relatively new to scripting in the Google Earth Engine and have an issue with calculating the Standard Vegetation Index (SVI) for a MODIS image collection. I want to write a function, which uses ee.Reducer.mean() and ee.Reducer.stdDev() of an image to calculate the SVI and map it over the image collection. However, GEE gives me the following Error-Code:

Error in map(ID=XXX_XX_XX): Image.constant: Invalid Image.constant type.

Where "ID=XXXX_XX_XX" seems to be a random date.
Has anybody an idea how to solve that problem?

//Add ROI
var ROI = ee.Geometry.Polygon(
  [[-2.2535247523754443,5.881894109545312],
  [-0.31993100237544425,5.881894109545312],
  [-0.31993100237544425,9.669268383910472],
  [-2.2535247523754443,5.881894109545312]]);

//Add MODIS ImageCollection
var MODIS = ee.ImageCollection("MODIS/006/MOD13Q1")
  .select(['EVI'])
  .filterDate('2000-01-01', '2020-01-01')
  .filterBounds(ROI)
  .map(function(image){return image.clip(ROI)});

//Function to Map SVI over image collection
var CalculateSVI = MODIS.map(function(SVI){
  var SelectImage = SVI.select(['EVI'], ['SVI']);

//Calculate Mean of Image
  var Mean = ee.Number(SelectImage.reduceRegion({
  reducer: ee.Reducer.mean(),
  geometry: ROI,
  bestEffort: true,
  }));

//Calculate Stdv of Image
  var Stdv = ee.Number(SelectImage.reduceRegion({
  reducer: ee.Reducer.stdDev(),
  geometry: ROI,
  bestEffort: true,
  }));

//Calculate SVI
  var CalcSVI =  SelectImage.expression(
  ((SelectImage.subtract(Mean)).divide(Stdv))
  );

  return CalcSVI.addBands(SVI);
});
print(CalculateSVI);


texts:
I have set up GeoServer on a mac mini running OSX Mavericks. I have my Mac hooked up to my LAN and would like to start developing geographic web applications.

I have a bunch of data sitting in my PostgreSQL/PostGIS database. I have a handful of map services set up in GeoServer. I even have a few nice styles put together to make everything look nice. I put together a little introductory web application the get the geography party started on my LAN but for the life of me, I can't seem to consume the geography remotely.

I access the service with a little bit of Javascript:

var map;

function mapinitialize() {
    bluemap = new L.TileLayer.WMS(GEOSERVERBASE +     "/geoserver/phillystreetmap/wms",
{
    layers: "phillystreetmap:Philly",
    format: 'image/png',
    attribution: ""
});

map = new L.Map('map',
    {
        center: new L.LatLng(37.45, -77.5),
        zoom: 12,
        layers: [bluemap],
        zoomControl: true
    });
}

and it looks great running locally on the server. However, when I try to access the site from another computer on my LAN, everything but the geography shows up. The error I get for the tiles is:

GET http://localhost:8080/geoserver/phillystreetmap/wms?service=WMS&request=Get…x=-8619650.805662755,4500612.225431179,-8609866.86604225,4510396.165051683 net::ERR_CONNECTION_REFUSED

Any idea why the data would show up when I run the application locally on the server but not when I access the site remotely?


texts: Styling the legend in QGIS with line layer

texts:
Working on QGIS 3.4.13.

I have a vector line layer, there are rivers. I put a rule based symbology with two categories and different scale max / min for each category.

I want to put style for the name of each type of river. I configured my legend in vector layer properties (following those steps : https://docs.qgis.org/3.4/en/docs/user_manual/working_with_vector/vector_properties.html?highlight=l%C3%A9gende#legend-properties)

I also checked this :
QGIS styling a legend

So at this point I have this in QGIS

and that in QGIS composer

But I would like to have my symbol AND my legend styling apparent, something like this :

seen here : Placing the label style alongside the layer style within the legend - Arcmap 10

I also checked in QGIS 3.10.1, to see if there some more functions, apparently,no.

Any idea?


texts: Creating line between two points with same value using QGIS?

texts:
I have two separate point layers. The first point layer contains the locations of
100 sensors which are aligned in one specific direction (for example: Asten).
The second point layer contains the location-coordinates of the possible directions.

The attribut table of a sensor node looks like following:

The attribut table of a direction node looks like following:

What I am trying to do is to compare each point of point-layer-1 with each point of point-layer-2. If the "location"-attribut of the point-layer-2 node is equal to the "direction"-attribut of the point-layer-1 node, those two points should be connected to a vector. The result should look like this:

The result vectors should be stored in a new line-layer I suppose. Is there
a way to do this?


texts:
Snapping with measurement tool on all layers/single layer not working.

this.featureLayer = new FeatureLayer(this.Khasraurl,
     {
         mode: FeatureLayer.MODE_ONDEMAND,
         // infoTemplate: infoTemplate,
         outFields: ["*"]
     });

var symbol = new esri.symbol.SimpleFillSymbol(esri.symbol.SimpleFillSymbol.STYLE_SOLID, new esri.symbol.SimpleLineSymbol(esri.symbol.SimpleFillSymbol.STYLE_SOLID, new dojo.Color([255, 0, 0]), 2), new dojo.Color([0, 0, 255, 0.20]));

this.featureLayer.setRenderer(new SimpleRenderer(symbol));
this.map.addLayer(this.featureLayer);

var layerInfos = [{
    layer: this.featureLayer
}];

// alert(featureLayer);
//alert(JSON.stringify(featureLayer));
//dojo.keys.copyKey maps to CTRL on windows and Cmd on Mac., but has wrong code for Chrome on Mac

var snapManager = this.map.enableSnapping({

    alwaysSnap: false,
    snapKey: 'CTRL', // key string ==> keys['CTRL']
    tolerance: 15
});

snapManager.setLayerInfos(layerInfos);


texts: Can i calculate viewshed in MS-SQL server? Viewshed analysis is Raster operation. I have elevation value vector observation as point vectors. Are there any methods to calculate Viewshed (using algorithm or something like that) ?


texts: How to extract planting lines from a raster image?

texts:
I have a high quality drone image and I want to remove the soil and leave only the planting line.So far already able to classify the images leaving 1 for the plantations and 0 for the soil, it is a gray image of simple band, I made the separation using the NIR image and using r.reclass.

The difficulty is to crop the RGB image to display only the crops, I do not know how to do that.I tried to vectorize the reclassified image to trim the RBG but it did not work, it's still in the process of vectorizing, it's taking too long!I do not know if it is possible to remove RGB soil using the raster calculator.

Any suggestions are welcome. I currently use QGIS and SAGA.
Here is an image to illustrate.


texts:
Experimenting with the W3C Geolocation API.
Is there a way to get the method how the location was processed? I mean GPS, IP, Cell IDs etc. I didn't find any object properties doing that.


texts: QGIS 2.18 - Can't get click coordinates from mapCanvas

texts:
Trying to get click coordinates from map canvas but code is not working. Newbie for Python, pyqgis, pyqt

I already check these topics but still code is NOT working

1 - QGIS 3 Get coordinates from mouse click

2 - Programatically check for mouse click in PyQGIS?

3 - How to get co ordinates of points on mouse click in PyQGIS

Here is my code that I am trying

import os

from PyQt4 import QtGui, uic
from qgis.gui import QgsMapToolEmitPoint
from qgis.utils import iface
from PyQt4.QtCore import *
from PyQt4.QtGui import *

FORM_CLASS, _ = uic.loadUiType(os.path.join(
    os.path.dirname(__file__), 'click_coor_dialog_base.ui'))

class ClickCoorDialog(QtGui.QDialog, FORM_CLASS):
    def __init__(self, parent=None):
        """Constructor."""
        super(ClickCoorDialog, self).__init__(parent)
        self.setupUi(self)
        self.pushButton.clicked.connect(self.showCoor)
    def showCoor(self):
        # this QGIS tool emits as QgsPoint after each click on the map canvas
        pointTool = QgsMapToolEmitPoint(iface.mapCanvas())
        pointTool.canvasClicked.connect(self.display_point)
        iface.mapCanvas().setMapTool( self.pointTool )
    def display_point(self,point):
        print point[0]

No matter what changes i do in above code, the display_point function does NOT get called, when I clicked on qgis map canvas. But it supposed to be called

What mistake I am doing here?


texts: How to use relative path in Mac for open internal and external file?

texts:
I would like that my project would be portable, I mean that I can share it without lose all the path I have in my project. This seems not so easy and many times I have to reset all the path when I give the project, with all the folders with rasters and shapefile into to a colleague.

So I would like to know:

can I set by default only relative path to the folder in which my project is (for opening raster and shapefile that QGIS can handle) and this remain relative even if I share the whole project to a colleague with a different operating system.
For opening an external file with action in properties->Action I use Open file with this action [% "Ref" %] where Ref is the attribute where I write the file path. I supposed that writing the path with  ./  before the path (e.g. ./RASTER/test.jpg) could be use to create a relative path but this seems not working. Why? And this however would work only in Unix like system!

I'm using QGIS 1.8.


texts: QGIS 2.6.1 doesn't load python plugins on Ubuntu after upgrade

texts: I have experienced this scenario once before and now today, whenever I updated my QGIS installation on Ubuntu, it didn't reload python plugins under ~/.qgis2/python/plugins when started, only c++  plugins are available. I can't really recall what I did previously but when I restarted my machine it seemed to work out of the box. Any ideas?


texts: PyQGIS script works on one computer but not on another?

texts:
The script below works perfectly fine in my computer. But when running the script on another computer, I receive an error:

  Error: Algorithm not found

The computer has similar specifications to my own (Windows 7, 64-bit, no admin restrictions). I did the following steps for the computer to match the scripts requirements:

Installed OSGeo4W 64-bit, QGIS 2.6.1
Copied the processing folder from C:\OSGeo4W64\apps\qgis\python\plugins\processing\ to C:\Users\user_name\.qgis2\python\plugins\

Did I miss something obvious?

import os, sys, glob

from qgis.core import *
from qgis.gui import QgsMapCanvas
from PyQt4.QtGui import *
from os.path import expanduser
home = expanduser("~")

# Set a custom config path
QgsApplication( [], False, home + "\AppData\Local\Temp" )
QgsApplication.setPrefixPath("C:\\OSGeo4W64\\apps\\qgis", True)
QgsApplication.initQgis()
app = QApplication([])

# Folder path of the Results for shapefiles
path_dir = home + "\Desktop\Test\\"
path_res = path_dir + "Results\\"

# Prepare canvas framework from qgis_interface.py
sys.path.append( home + '/Desktop/Test//' )

# Get an iface object
canvas = QgsMapCanvas()
from qgis_interface import QgisInterface
iface = QgisInterface( canvas )

# Prepare processing framework
sys.path.append( home + '/.qgis2/python/plugins' )

# Initialize the Processing plugin passing an iface object
from processing.ProcessingPlugin import ProcessingPlugin
plugin = ProcessingPlugin(iface)
from processing.tools import *

Cellsize = 1000
layerPath = path_dir + "Layer.shp"
extent = QgsVectorLayer( layerPath, '', 'ogr' ).extent()
centerx = (extent.xMinimum() + extent.xMaximum()) / 2
centery = (extent.yMinimum() + extent.yMaximum()) / 2
width = extent.xMaximum() - extent.xMinimum()
height = extent.yMaximum() - extent.yMinimum()

def run():

    outputs_1=general.runalg("qgis:creategrid", Cellsize, Cellsize, width, height, centerx, centery, 1, 'EPSG:7405', path_res  + "/"+ fname)

run()
QgsApplication.exitQgis()
app.exit()

It appears that the script cannot find any processing algorithm (from qgis, GRASS or SAGA). It also does not recognise PyQt4.Core whereas it does on my machine. Could this be relevant?


texts: How to change the coordinate transformation function in Oracle?

texts:
I'm trying to reproject some data from one CRS to another (SRID:27493 to SRID:3763) using the SDO_CS.TRANSFORM function.

The function works and returns no errors, but the coordinate transformation function that Oracle is using must be inadequate, because the results are inaccurate.

So, I would like to alter or define a new coordinate transformation function that suits my data.

I've realized that in the SDO_COORD_REF_SYS table there is an PROJECTION_CONV_ID attribute that should lead somewhere, but the transformation function and its parameters seem to be dispersed through several tables.

Anyone has experience in this and / or can point some manuals / tutorials?

Problem 1.

One of the problems I had with the coordinate transformation between data in 'Datum 73 / Modified Portuguese Grid' (SRID:27493) and ETRS89/PT-TM06 (SRID:3763), was that the results I was getting were different of the ones I would get through PostgreSQL/PostGIS, which uses Proj4 libraries and gives identical results than the default coordinate transformation assume by EPSG.

I realized in time that Oracle was using deprecated parameters for the Tranverse Mercator projection from Datum73 - SRID:4274 (the geodetic coordinate system in which SRID:27493 is based) and WGS84 - SRID:4326, and that was influencing the projections of any coordinate system based on Datum73.

Problem 2.

The coordinate transformation method recomended by Direcção-Geral do Território (DGT) between SRID:27493 and SRID:3763, therefore the most accurate, is based on NTv2 transformation and the Grid Shift Matrices provided by DGT.

So I would like to had this new transformation method and its parameters to Oracle Spatial.

NOTE: the Oracle documentation mentioned in comments was no definite help on this


texts: I am using R for geostatistical analysis and spatiotemporal variogram calculations on a server. Although it has very high specifications but using variogramST takes a very long time (12 hours so far and still running) because I'm working on a relatively big data set. I have found some R packages that use CUDA. Is there any GPU support for geostatistics packages like gstat and spacetime?


texts: Changing highest and lowest value of image using ArcGIS for Desktop?

texts:
I want to change the highest and the lowest value of raster image. I want it to have the highest value of 345 and the lowest of 0. The highest value of original image now is 265 and the lowest is 0.02.

Is there any way to do this using ArcGIS?


texts: How to solve an "Error found when retrieving conflicts." error in ArcMap?

texts:
I uncovered this error when testing ArcMap 10 SP2 for a customer:

This occurs when the conflict is detected--ArcMap is attempting to show the conflicts window but fails.

I get this error on all workstations at the customer's location (all Windows XP).  But I don't get the error on my systems (I've tried XP and Windows 7).

I've tried various oracle/sql server combinations and I still receive the error on the customer's workstation.

It appears to be a compatibility issue and I have an open Esri incident which hasn't yielded any results, so I thought I'd tap the gis.se.

Has anyone experienced this error with this set(or similar set) of conditions?

Update From Comments:
I initially got this in an Oracle10/9.3.1 sde database, but also replicated the error with a sqlexpress/arcsde 10.  I only have this issue with the customer's ArcGIS Desktop workstations and each experience the same issue.

Data from crash:

00:04:14.578: First chance exception 0xC0000005 (Access Violation) occurred in "c:\program files\arcgis\desktop10.0\bin\VERSIONUI.DLL" at address 0x39458983.
00:04:14.578: LoadLibraryW("DbgHelp.dll") called from "c:\program files\arcgis\desktop10.0\bin\AFCORE.DLL" at address 0x00714A04.
00:04:14.594: Loaded "c:\program files\arcgis\desktop10.0\bin\DBGHELP.DLL" at address 0x0B770000.  Successfully hooked module.
00:04:14.594: DllMain(0x0B770000, DLL_PROCESS_ATTACH, 0x00000000) in "c:\program files\arcgis\desktop10.0\bin\DBGHELP.DLL" called.
00:04:14.594: DllMain(0x0B770000, DLL_PROCESS_ATTACH, 0x00000000) in "c:\program files\arcgis\desktop10.0\bin\DBGHELP.DLL" returned 1 (0x1).
00:04:14.594: LoadLibraryW("DbgHelp.dll") returned 0x0B770000.
00:04:14.610: GetProcAddress(0x0B770000 [c:\program files\arcgis\desktop10.0\bin\DBGHELP.DLL], "MiniDumpWriteDump") called from "c:\program files\arcgis\desktop10.0\bin\AFCORE.DLL" at address 0x00714A14 and returned 0x0B7CD4F0.
00:04:14.610: Second chance exception 0xC0000005 (Access Violation) occurred in "c:\program files\arcgis\desktop10.0\bin\VERSIONUI.DLL" at address 0x39458983.
00:04:14.657: Exited "c:\program files\arcgis\desktop10.0\bin\ARCMAP.EXE" (process 0xC0C) with code -1073741819 (0xC0000005).

New Information (8/4)

Client is getting this in ArcGIS 9.3.1 as well--they so seldom get
conflicts that this hasn't been an issue until now.
I'm nearly certain that this issue is caused by some other software that was pushed out to these machines.  We're still troubleshooting and I'll update this question with an answer once I know more.


texts: moveToLocation and open Marker-Popup automatically

texts:
For reference: I use L.Control.Search to search for my Marker.

It's this Plugin: https://github.com/stefanocudini/leaflet-search/blob/master/examples/multiple-layers.html

And this is the working search function: https://labs.easyblog.it/maps/leaflet-search/examples/multiple-layers.html

My problem:
Now, I want the popup of the marker I searched for to open automatically without me clicking. I butchered the code from a polygon layer, but can I make it work for only markers? I don't get any errors but don't get a popup either.

       var searchControl = new L.Control.Search({
            layer: myLayer,
            propertyName: 'name',
            collapsed: false,
            container: 'findbox',
            textPlaceholder: 'Search',
            marker: false,
            initial: false,
            buildTip: function (text, val) {
                var type = val.layer.feature.properties.iconcategory;
                return '' + text + '' + type + '';
            },

          moveToLocation: function (latlng, title, map, feature) {
                var zoom = 16;
                map.setView(latlng, zoom);
            },
        });

        searchControl.on('search:locationfound', function(e, feature) {

    if(e.layer._popup)
        e.layer.openPopup("test");

}).on('search:collapsed', function(e) {

    filterLayer.eachLayer(function(layer) {
        filterLayer.resetStyle(layer);
    });
});

        map.addControl(searchControl);

This is a glimpse of the code I used for the popups.

onEachFeature: function forEachFeature(feature, layer) {

            var popups = "Name " +
                feature.properties.name + ";

            if (feature.properties && feature.properties.popups) {
                popups += feature.properties.popups;
            } layer.bindPopup(popups);

        },


texts: PyQGIS composer add picture by layer path field

texts:
I want to create a atlas composer map using PyQGIS3.
How can I take dynamic path by field from atlas layer?

In PyQGIS3 I need some like this :

pict = QgsLayoutItemPicture(layout)
pict.setPicturePath(myAtlas['pict'])
pict.setFrameEnabled(True)

where myAtlas is my layer and pict is my field with paths.


texts: Executing Dissolve method in ArcGIS Engine?

texts:
I'm always having a problem in GP.execute() when trying to execute this method. Here is my code:

 Friend Function DissolveFeatures(ByVal pInFLayer As ESRI.ArcGIS.Carto.IFeatureLayer, ByVal DissolveField As String, ByVal OutName As String) As ESRI.ArcGIS.Carto.IFeatureLayer

            Dim GP As New ESRI.ArcGIS.Geoprocessor.Geoprocessor
            Dim Dissolve As ESRI.ArcGIS.DataManagementTools.Dissolve = New ESRI.ArcGIS.DataManagementTools.Dissolve
            Dim pFLayer As ESRI.ArcGIS.Carto.IFeatureLayer = New ESRI.ArcGIS.Carto.FeatureLayer

            Try
                GP.AddOutputsToMap = True
                GP.OverwriteOutput = True
                GP.SetEnvironmentValue("workspace", "Database Connections\IFE connexion.sde")

                Dissolve.in_features = pInFLayer
                Dissolve.dissolve_field = DissolveField
                Dissolve.out_feature_class = "Database Connections\IFE connexion.sde\" & OutName

                RunTool(GP, Dissolve, Nothing)

                pFLayer.FeatureClass = CType(GP.Open(OutName), IFeatureClass)
                pFLayer.Name = pFLayer.FeatureClass.AliasName

            Catch ex As Exception
                MsgBox(ex.Message & "Dissolve Features")
            End Try

            Return pFLayer

        End Function

Private Shared Sub RunTool(ByVal geoprocessor As ESRI.ArcGIS.Geoprocessor.Geoprocessor, ByVal process As IGPProcess, ByVal TC As ITrackCancel)

        ' Set the overwrite output option to true
        geoprocessor.OverwriteOutput = True

        Try
            geoprocessor.Execute(process, Nothing)
            ReturnMessages(geoprocessor)

        Catch err As Exception
            Console.WriteLine(err.Message)
            ReturnMessages(geoprocessor)
        End Try
    End Sub

    ' Function for returning the tool messages.
    Private Shared Sub ReturnMessages(ByVal gp As ESRI.ArcGIS.Geoprocessor.Geoprocessor)
        ' Print out the messages from tool executions
        Dim Count As Integer
        If gp.MessageCount > 0 Then
            For Count = 0 To gp.MessageCount - 1
                Console.WriteLine(gp.GetMessage(Count))
            Next
        End If
    End Sub

This returns an error:

  Error HRESULT E_FAIL has been returned from a call to a COM component


texts: Highlight Polygons on a map Arcgis API for Javascript

texts:
I am able to query a feature layer using the following code

          var query = new Query();
          var geometry = null;
          query.returnGeometry = true;
          query.outFields = ["*"];
          query.where = "CONST LIKE '%"+searchString+"%'";
         resultsLayer.queryFeatures(query).then(function(results){
              var resultsArray = results.features;
              console.log();
          });

However after collecting the results I would like to be highlighting the returned features on the mapView. The results however didn't return a 'graphics' attribute for me to add to the map layer which is why I am stuck.

How do I add features to the mapView using the returned features array?


texts: Using search and update cursors for random number generation in ArcPy?

texts:
I am having some trouble getting a search cursor to return a value for each record. The purpose of my code is to scan through a list of records in a field (FieldA) and based on those values populate a second field (FieldB) with a random number between 1 and the value in FieldA. FieldA contains numeric values which vary depending on their association. My code is as follows:

CursorA = arcpy.SearchCursor(FC)
for row in cursorA:
    rowcount =(row.getValue("FieldA"))
    row = cursorA.next()
    cursorB = arcpy.UpdateCursor(FC, "FieldB")
    for rowB in cursorB:
        rowB[0] = randint(1, rowcount)
        cursorB.updateRow(rowB)

When I run the code it appears as though the searchcursor locates the value of the last record in FieldA and uses that value as the highest value in the random sequence. I could be wrong on this however, and it could be related to improper coding on my part.


texts: Getting colour ramps and other funky styles in the ArcGIS REST API's DrawingInfo or elsewhere

texts:
I have some layers with "funky" styling, they contain colour ramps, or hashes across the line (mimicking train tracks), or dashes with different spacings. I'm reading all of these layers from an ArcGIS map server and using the drawingInfo to generate Mapbox styling objects.

However while I can get the renderer type (simple, Unique value, and Class breaks) and line/circle/fill/etc type, I can't see these funky additions anywhere in the returned drawingInfo objects from the rest API. I can see that color ramp objects definitely exist.

How am I able to programmatically discover these extra styling flourishes on my layers?


texts: Esri enterprise geodatabase migration to PostGIS/PostgreSQL using FME

texts:
I am trying to migrate my Esri enterprise geodatabase (10.2.2) on MSSQL server 2012 to PostGIS version 2.0/PostgreSQL version 9.4 using trial version of FME. I am getting the following errors;

bulk copy failed on table 'table name' using delimiter ':'. Error was Error column has Z dimension but geometry doesnot Context: Copy 'table name', line 1, column geom:''

Any suggestions on how I can resolve this issue?


texts: Convert CSV GPS Logs to GPX on Android (GPSBabel on Android)

texts:
I bought a GPS Logger - a small device that will, while turned on, log the GPS Coordinates to a SD card every second. It has no display and no other UI than 3 buttons and a couple of LEDs.
The Device will write its log into a CSV files with the coordinates + meta data.

GPSBabel was a great help to get this csv into a GPX XML file very quick and easy.

However, as i will go on a hiking trip and i'd like to see my track while hiking, i would need an app that i could run on my Android phone or tablet (can't carry my linux laptop ;)).

Is there any app that might have this as a secondary feature that is so far at the end of the feature list that Google can't find it?


texts: Create Georeferenced Raster For JPEG file with Mapnik Python

texts:
I am able to use Mapnik generate_image.py to generate image with a BBOX and projection.

Also, I have added a python script to generate a JGW file. But to completely georeference the generated image. I guess I need to add projection information to it.

I found out it can be done via .prj file. But not sure if it can be generated via Python script. Is there any way to do the same? Is only generating .jgw file the way to completely georeference the image?


texts: Re-project raster image from Mercator to Equirectangular

texts:
I have a raster image of the world map in Mercator projection, square 1:1 format. I want transform the image into an Equirectangular projection accurately.

I assumed from my research that GDAL would be the way to go, specifically gdalwarp. I used gdaltranslate to add gcp points to the image on the four corners, and then used gdalwarp to make the transformation.

It produced an image with the correct size – 2:1 ratio – but the projection has not been accurately transformed. It looks exactly the same as if I had taken it into Photoshop and simply stretched it to that shape. The proportional relationship between Mercator and Equirectangular is obviously more complicated than a simple, global stretch like this.

So my question is: is there a way to accurately transform between projections for raster images?

Edit / Additional Information

Here is the starting image: http://imgur.com/48wtsgn

I then applied to GCP points using the following gdal_translate command:

gdal_translate -gcp 0.0 0.0 -180.0 90.0 -gcp 1000.0 0.0 180.0 90.0 -gcp 1000.0 1000.0 180.0 -90.0 -gcp 0.0 1000.0 -180.0 -90.0 Mercator.tif MercatorGCP.tif

Then I used the following gdalwarp command to make the intended reprojection from Mercator (ESPG:3857) to Equirectangular (ESPG:32662):

gdalwarp -s_srs EPSG:3857 -t_srs EPSG:32662 MercatorGCP.tif Equirectangular.tif

Here is the resulting image: http://imgur.com/5LGk16o

Compare it to a map with proper Equirectangular projection, like this: http://imgur.com/3ncIW2R

You can see that the proportions have not been properly adjusted, the map not actually reprojected at all, simply squashed uniformly.

Am I doing something wrong, or is this simply not possible with gdal?
If not, is there another way to do this?

Thanks!


texts:
I work on a particular area, my data are geospatial types and are stored in a geodataframe which contains:

a structure column (house, building, town hall, store etc).
a neighbors_structure column which has been deduced from the graph representing the n closest structures for each structure.
a geometry column, of type POINTS in my case. All structures are located and are an ID that I do not put here for the sake of brevity. Also, the GPS coordinates given here are hypothetical.

For example, with n=2 we can have:
   structure            neighbors_structure                    geometry
0  house                      house - house     POINT (42.70275 9.94481)
1  building        store - building - house     POINT (42.70030 9.94783)
2  house           building - store - house     POINT (42.69700 9.97133)
3  house                   store - building     POINT (44.18208 10.27975)
4  house              house - house - house     POINT (43.17099 10.52153)
...

Using this geodataframe, I would like to answer the following types of questions:

If we have a structure of type _ there is a _% chance that its closest neighbor is a _

For example : if we have a structure of type house there is 30% chance that its nearest neighbor is a building.

If we have a structure of type _ there is a _% chance that in a window of _km2 we find at least one _
If we have a structure of type _ there is a _% chance that within a radius of at least _km we find another structure

Does a powerful tool to do this exist?
I'm not sure that Pysal for example allows to do exactly that with simple POINTS.


texts: Can't choose interpolation attributes of my shapefile

texts: I am attempting to interpolate data of a point shapefile. In the attribute table all attributes are shown correctly. When I am trying to interpolate them with Raster >> Interpolation no attributes to choose are shown.


texts: Web Tier Authentication in ArcGIS for Server?

texts:
Can anyone explain me the steps, how to view map service from LDAP ,Web Tier Authenticated ArcGIS Server 10.2.1 in .net Application.

I have installed web adapter in Web server with windows authentication enabled.

When I open my application in my login I can see the map but with other system I can't see the map.

I am using Forms authentication in .net application.

Did I miss any configuring?


texts: Build network based on shortest distance between points

texts:
Is there a way to build a network based on the shortest distances between each point of a point layer.

I have occurrence data from species (as points) and want to estimate whether there are one or more populations. As I know (more or less) the potential dispersal ability of the species I can estimate whether the occurrence points could be connected populations or whether they are isolated from each other.


texts: Create an animation of movement along a shortest path

texts: I am new using QGIS. In my project I need to find the shortest path between different layer points. I did this using the network analysis shortest path tool for single points and OD-matrix tool with QNEAT 3 plugin for multiple points. Now I want to show vehicles' movement through the shortest path found. I have already found how to extract the list of points through the shortest path route. How can I use them to represent a continuous movement of vehicles in the route?


texts: Making reclassified raster layers automatically display break values in their attribute tables?

texts:
I am using ArcGIS 10.2. I have many raster layers which I need to reclassify (using the quantile method), I Also need to know the break values of these reclassifications.

I know that the break values can be viewed in the reclassify dialog box, however, due to the sheer number of layers I am dealing with, I do not have time to copy these values from this dialog box, before copying them again into the attribute tables of the reclassified layers.

Is there some way that I can have the break values automatically displayed as a field of the attribute table  belonging to the reclassified output?


texts: Getting boundary of raster image as polygon in ArcGIS Desktop?

texts:
I have a rectangular raster image in ArcGIS Desktop.

How can I produce a polygon which is just the extent of the raster image?

That is, I want to have a polygon layer with a single quadrilateral which is the edge of the image.

I've tried simply doing a Raster->Polygon conversion, but this tries to use the data in the image to generate the outline - all I want is the outline of the edge of the image.

I'm using ArcGIS Desktop 10.


texts:
How to install postgis2.0 on Debian6.0 and with Postgressql9.1 ?
Didn't find any good repository ?
Do you have to go all the way and compile ?


texts: How to fix a gap in a Google Transit trip that spans the beginning/end of a loop?

texts:
In my GTFS data, I have bus routes that are on a loop and that is specified in my frequencies.txt. However, I noticed when querying trips on Google that if the trip crosses over the starting/end point of the loop it doesn't connect and it tells you to walk.

Query example

My shapes.txt file draws out the full loop for the bus routes so I would expect everything to work but it always seems to have a gap from the last stop in stop_times.txt to the first.

Here are my GTFS files. https://github.com/joshuahouston/theRIDE/tree/gh-pages/feed_files


texts: How to interpret spatial random permutations against the observed value of the statistic (PySAL)

texts:
I have a map for which I calculate the value of the Global Moran's I. With the use of PySAL, there is also an option to do a computational inference. The result of such analysis is shown below.

From the videos I've seen and some papers I've read, what I understood from this is that random spatial permutations work like this:

I take my original map, shuffle the arrangement of the polygons
For each iteration, calculate the Moran's statistic
Do N times
Fit a curve to see the distribution of the simulated statistics
compare with the observed value of the statistic

In the figure above, my observed statistic is far to the right of the reference distribution.

Does this mean that, by simulating randomness, none of the configurations were able to come up with a value similar to that of the observed statistic? If so, does this mean that my map presents non-random patterns? Or does it mean that it is not probable to even have such  a configuration/pattern?

Thanks for any info on this!


texts: Buses movement data visualization in QGIS

texts:
What is the easiest way to visualize the movement of city buses in QGIS, as it is shown at this website?  https://anitagraser.com/2020/05/24/movement-data-in-gis-30-synchronized-trajectory-animations-with-qgis-temporal-controller/

There is a precise timetable for buses movement from the official site of public transport service of my city. It contains the information for each bus to bus-stop arrival time. Also I'm going to create a line layer of roads and point layer with bus-stops (I can create a point layer with points, laying on the roads near the bus-stops (using snapping options)) to show precisely places, where the buses exactly stop.

How  can I create such an animation (as at Anita Graser's website) with not very much difficulty?


texts: How to prevent automatic new line when labeling point feature class?

texts:
I have a point feature class that I am labelling with the string "EXIT" for a majority of points, but some points will be longer, such as "EXIT into Pacific".

I wrote a labeling script that will create a new line after "EXIT", but not after "into" because I want to keep "into Pacific" on one line.

I want the label to show up like this:

     EXIT
 into Pacific

But instead, it automatically creates a new line and shows up like this:

     EXIT
     into
    Pacific

Is there a way to keep ArcMap from automatically creating a new line for the label? Or is there a good invisible character that will make the program think that "into Pacific" is one string with no white space?

Here's my script:

function FindLabel ([TYPE]){
var exitSlice
var intoSlice
intoSlice = ""
var label
var field
field = [TYPE]
exitSlice = field.slice(0, 4)
if (field.length > 4) {
    intoSlice = field.slice(5)
    label = exitSlice + "\n" + "<ITA><FNT name='Melior' size='8'>" + intoSlice + "</FNT></ITA>"
    }
else {
    label = exitSlice
    }
return label}

I'm pretty confident the script isn't causing the extra newline between "into" and "Pacific".


texts: Census Tracts entirely within another Census Tract

texts: I am working with a U.S. Census Tract layer. I am just wondering why there are tracts entirely within another tract when, in general, the territory is divided by partitioning the space into a number of simple rectangular regions, divided up by axis parallel (more or less) splits.


texts: Do ArcGIS clients that direct-connect to postgresql/postgis databases override configuration settings?

texts:
I went in to modify some settings on a Postgresql database via the PGAdmin config editor and some of the settings don't have their default values and they are higher because they've been overridden (e.g.  maintenance_work_mem ).  I know PGAdmin will do this. And that I can do it manually in PSQL.

Do the ArcGIS clients ever override postgresql configuration settings?  And if so, what triggers it? Is it tool-based? Or when the software starts up? Or is Esri entirely hands-off?


texts: st_combine and st_union produce empty attribute table for shapefile (how to edit before st_write?)

texts:
I am starting with R programming and I am having a hard time figuring out how to edit the attribute table of a shapefile after using st_combine and st_union.
My shapefile has multiple polygons and I want to create a shapefile with a unique outline and a well-defined attribute table like this one:
| A | B | C |
|---|---|---|
| X | Y | Z |

This is part of my code:
# Combining features into one
Shapefile <- st_read(Shapefile2Process)
st_combine(Shapefile)
Outline <- st_union(Shapefile)

# Simplification depending on TolInMeters.
Outlinetol <- st_simplify(Outline, preserveTopology = FALSE, dTolerance = TolInMeters)

# Save new shapefile with outline feature
st_write(Outlinetol, NewShapefile)

print(names(Outlinetol)) returns NULL, the shapefile is created.
I want to add columns with values to the empty attribute table but I don't know if this is possible with a sf object and before st_write, and how to do it?
I tried to convert this sf object to a "Spatial" object using as(Outlinetol, "Spatial"), create Outlinetol$A <- "X" and then writeOGR(Outlinetol,NewShapefile,driver="ESRI Shapefile") but I was not successful...


texts: Calculate real distances based on OSM/Google roads between many points ArcGIS

texts:
In ArcGIS I've two shapefiles with many points called "source" with hundreds of points, and another called "destinations". Based on the roads data of OSM/Google/other, I'd like to calculate the real distances between each point of the couple source-destination.

I've done some searchig, but I found only a solution for linear distances (tool Point distance or similar) o using network analysis but for one point at time.


texts:
I am wondering if there was an easy way to assign id's to existing points, but not in a manual way and not consistenly.

For example i have this file which contains 2 points in one position:

I need the first point (and second but in the same position) to be numerated 1 and 80, then 2 and 79 and so on.


texts: Updating a subset of raster values using open source Python?

texts:
I have a source raster data set (tif file) and a number of polygon features. For each one of those polygons I need to:

a) 'clip' the raster.
b) change values of those clipped raster cells (with the average value).
c) write new values back to the source raster.

I've managed to do a and b, which means I now have a numpy.ndarray variable where all cells have been populated with the average value. How can I write this back to the source raster?

I believe I'd be able to write a script which could do such a thing. I have arrays dimensions (width and height) and top left cell's position, so I guess it'd just be a matter of operating with these arrays. However, I was wondering whether there's any existing function that does that already.

My code:

import os, sys, datetime, time
import geopandas as gpd
import gdal
import rasterio
from rasterio.mask import mask
from fiona.crs import from_epsg
import numpy as np
import pycrs

alr_path = r'Z:\GRAU_Team_Admin\Alfonso Jimenez\Flood Modelling\Building Footprint\OUTPUT\NI_RIVER_ALR_1000.tif'

gdb_vml_polygs = r'\\lwukwvdi11\data\Data\UK\UK_Mapping\VML\data\Europa20190403\RSA_VML_BUILDINGS_NI.gdb'
fc_vml_polygs = r'VML_Buildings_NI'

alr = rasterio.open(alr_path)

vml_polygs = gpd.read_file(gdb_vml_polygs, driver='FileGDB', layer=fc_vml_polygs)

for index, row in vml_polygs.iterrows():
    #row[3] contains the geometry
    window, out_transform = mask(alr, row[3], all_touched=True, crop=True)
    if np.all([window < 0]):
        avg = -1
    else:
        avg = window[window != -1].mean().item()
        window_avg = (np.where(window!=-1, int(round(avg,0)), window))

In the script window_avg is a numpy array containing the new values. How can I write these values back into the source raster?


texts: How to draw line in leaflet using postGIS road data

texts:
I have few points from server and I want to draw a line with this point on road. But the leaflet drawing point to point line as shown in picture in orange color. But I want to draw the line according to the blue line. How it can be done?

Thanks


texts:
Does the a wms layer in openlayers have "maxfeatures" option?

I'm asking this because my wms layer contains allot of points ( > 100.000). When zoomed in I don't want the rest of the features (not sure if you call it features when using wms) to be rendered.

When I was using WFS, this was my code:

    wfs = new OpenLayers.Layer.Vector("WFS", {
        strategies: [new OpenLayers.Strategy.BBOX({
                        resFactor: 1,
                        ratio:1
                    })],
        protocol: new OpenLayers.Protocol.WFS({
                        maxFeatures:1000,
                        url:  "my url",
                        featureType: "my_featuretype",
                        featureNS: "my_feature ns",
                        version: "1.1.0"
                    })
});

How can I accomplish the same render functionality with a wms layer?

Edit

Solution: Hide the layer at a certain zoomlevel. This prevents a big amount of features to be loaded.

code:

    map.events.on({ "zoomend": function (e) {
        if (this.getZoom() > 7) {
            layer.setVisibility(true);
        }
        else {
            layer.setVisibility(false);
        }
    }
});


texts: Map slightly off-center with simple (cartesian) CRS and tilelayer in Leaflet

texts:
After several attempts, some tinkering and a lot of asking here, I got my fictional map running with tiles by setting up a simple custom CRS.
Alas, the map is slightly off center. I got the feelig that fitToBounds does not really work and if I add setView it does not change anything either.

My map size is 1024x1024 so the map should fit nicely at zoom level 2, but look for yourself.

The code looks like this:

 ////////////////////////////////////////////////////////////////////////////
    // Conversion from (x, y) raster image coordinates to equivalent of latLng
    // Taken from Leaflet tutorial "Non-geographical maps"
    // http://leafletjs.com/examples/crs-simple/crs-simple.html
    ////////////////////////////////////////////////////////////////////////////
    var yx = L.latLng;

    var xy = function(x, y) {
      if (L.Util.isArray(x)) { // When doing xy([x, y]);
        return yx(x[1], x[0]);
      }
      return yx(y, x); // When doing xy(x, y);
    };

    var minZoom = 2;
    var maxZoom = 12;
    var img = [
      16384, // original width of image `karta.jpg` => x / ~longitude (0 is left)
      16384 // original height of image => y / ~ reverse of latitude (0 is top)
    ];

    /**
     * See https://stackoverflow.com/questions/34638887/leaflet-custom-coordinates-on-image
     * http://leafletjs.com/reference-1.0.3.html#transformation
     * 256 px * 2^4 = 4096 px
     */

    L.CRS.MySimple = L.extend({}, L.CRS.Simple, {
      //                      coefficients: a      b    c     d
      transformation: new L.Transformation(1 / 48, 0, 1 / 48, 0) // Compute a and c coefficients so that  tile 0/0/0 is from [0, 0] to [img]
    });

    var bounds = [xy(0, 0), xy(img)];

    var map = L.map("map", {
      crs: L.CRS.MySimple, // http://leafletjs.com/reference-1.0.3.html#map-crs
      maxBounds: bounds, // http://leafletjs.com/reference-1.0.3.html#map-maxbounds
      minZoom: minZoom,
      maxZoom: maxZoom
    });

    map.fitBounds(bounds);  //this makes the map fit the bounds we defined

    L.tileLayer('tiles/{z}/{x}/{y}.png', {
      bounds: bounds, // http://leafletjs.com/reference-1.0.3.html#gridlayer-bounds
      noWrap: true,
      attribution: 'Your Map',
      minZoom: minZoom,
      maxZoom: maxZoom,
      maxNativeZoom: 6  //how deep the actual map tiles go
    }).addTo(map);

//Experimental Zoom Viewer

    var ZoomViewer = L.Control.extend({
        onAdd: function(){
            var gauge = L.DomUtil.create('div');
            gauge.style.width = '200px';
            gauge.style.background = 'rgba(255,255,255,0.5)';
            gauge.style.textAlign = 'left';
            map.on('zoomstart zoom zoomend', function(ev){
                gauge.innerHTML = 'Zoom level: ' + map.getZoom();
            })
            return gauge;
        }
    });
(new ZoomViewer).addTo(map);


texts:
Can anyone point me to instructions for setting up PyCharm (or any IDE) on OS X to run QGIS for building standalone scripts/applications?

I have seen many tutorials for Windows, but nothing out there (that I can find) for OS X.


texts: Accessing single value in field using ArcPy with ArcGIS Pro

texts:
I'm pretty new to Python scripting (especially in ArcGIS). I would like to access a single element from a field and modify it. I understand that I can do this with cursors by searching the field column for the value (arcpy.da.SearchCursor) and then replacing it (arcpy.da.UpdateCursor), but this feels inefficient.
Is there a way to directly access an element without looping continuously?


texts:
I would like some external users to be able to access WMS-services from my internally enabled GeoServer. By opening up my firewall for this specific IP-address it should get access to my WMS-services, right?  I'm sure I'm missing something crucial here since this is new to me.

Is there any other settings/permissions to be set in GeoServer for this specific layer/store/workspace to work out properly? The URL is working fine in internal applications.


texts:
I have a particular .asc file, with elevation data in it.
I am wondering if somebody knows, what is the unit for the "cellsize" variable in the header of the .asc file?

Here is how the header looks like:

ncols   236
nrows   121
xllcorner   2.022083333
yllcorner   41.32875
cellsize    0.000833333
NODATA_value    -32768

What unit is the "0.000833333"?
A meter, or an arc?

I found a couple of pages on internet which explain the content of the .asc files headers, but still none explain what is the length unit of the "cellsize"?

Does anyone know?


texts: Get error when QgsVectorLayer is not valid

texts:
In QGIS 2.18.16, with a custom plugin, I tried this code to create a new vector layer for WFS:

v_layer = QgsVectorLayer(uri.uri(), "my WFS layer", "WFS", False)
if not v_layer.isValid():
    print v_layer.error()            # 1
    print v_layer.error().message()  # 2
    print v_layer.error().summary()  # 3

print 1: is a qgis._core.QgsError object
print 2: empty string
print 3: qgis crash with segmentation fault

I need to get an error type or the WFS response to understand, for example, if the url is not correct or if the service is down.


texts: Items (scripts) in custom toolbar are not working

texts:
I created custom toolbar in my working mxd and added few menus to it:

Using customize dialog I added some system and my own tools (scripts) to it:

It worked great for last two months or so, providing very fast access to most frequently used tools. However yesterday after long 4 days weekend I found that the only items in menu still working are the ones pointing to system tools. Nothing happens when I click on the item pointing to my own tools.

I did not:

remove tools from COMMANDS using CUSTOMIZE dialog. Moreover if I click description within that dialog, it gives me correct tool description
physically moved scripts.

What went wrong? How to prevent it from happening again?


texts: why a vertical white band appears after saving/exporting from print composer?

texts: I have a problem with my map composer. After exporting images a white strip appears on my map. How can I fix that?


texts: Missing Dependency in SAGA. Orfeo Toolbox does not work

texts:
I am trying to perform segmentation from Orfeo Toolbox (QGIS last version 2.2), whenever I try to use the toolbox, I get the following message:

Oooops! The following output layers could not be open
Output labeled image: D:/ERDAS_examples/Sardegna
1954/segmnentation.tif
The above files could not be opened, which probably indicates that
they were not correctly produced by the executed algorithm
Checking the log information might help you see why those layers were
not created as expected

The same when I try to use SAGA toolbox:

Missing dependency.This algorithm cannot be run :-(
This algorithm requires SAGA to be run.Unfortunately, it seems that
SAGA is not installed in your system, or it is not correctly
configured to be used from QGI,
Click here to know more about how to install and configure SAGA to be
used with QGIS

How do I solve this?