texts: expression for "SelectLayerByAttribute_management" with value of a cursor

texts:
I am currently struggling with the correct syntax for a SelectLayerByAttribute_management query, which is created iteratively by a cursor. I want to intersect each municipality from a polygon file containing the municipality boundaries with another polygon file (using arcpy.SelectLayerByLocation_management - it works so far). Currently my problem is that I cannot pass the cursor value containing the municipality name into the query.

arcpy.MakeFeatureLayer_management(county, "countyLyr")
arcpy.MakeFeatureLayer_management(LASlist, "LASlistLyr")
reset_layers("countyLyr", "LASlistLyr")

with arcpy.da.SearchCursor("countyLyr", ("NAME", "ICC", "BEZIRKSNUM"), """"ICC" = 'CH' AND "BEZIRKSNUM" = 'CH01040000'""") as cursor:
    for row in cursor:
        #select = (row[0])
        #expression = "\"NAME\" = 'Steinmaur'"
        county = row[0]
        print(county)
        expression = "\"NAME\" = str(county)"
        print(expression)
        arcpy.SelectLayerByAttribute_management("countyLyr",'NEW_SELECTION', expression)
        print(row[0])
        matchcount = int(arcpy.GetCount_management("countyLyr")[0])
        print(matchcount)
        if matchcount <> 0:
            arcpy.SelectLayerByLocation_management("LASlistLyr",'INTERSECT', "countyLyr")
            matchcountlas = int(arcpy.GetCount_management("LASlistLyr")[0])
            print(str(matchcountlas) + ' LAS-Kacheln')
            #do something with the result
        else:
            print('Keine Selektion / Keine Daten')

If I put the name of the county in the form expression = "\"NAME\" = 'Steinmaur'" into the query , it works perfectly (but only for this county). When I output the cursor value print(row[0]), I also get the correct result (all counties). But when I combine the query with the cursor value I got unfurtunately the error message

  ExecuteError: ERROR 000358: Invalid expression

So my hopefully easy question how I can combine the cursor value with my querry, that I won't get the error message?


texts: Read and process point vector data in Spyder 3.7 without using GeoPandas?

texts:
I read the file in my Spyder 3.7 environment (without Anaconda) using Pandas and I need to convert the CSV file to a shapefile. I will then make circular buffers around the points but I do not want to/cannot use GeoPandas for this.
Any suggestions on the best possible way to go about this? matplotlib, numpy, pandas, scipy and sympy are the only optional modules that I have.


texts: Cesium Terrain Builder : Could not open GDAL Dataset

texts:
I am trying to run Cesium terrain Builder on Ubuntu 14.04.
On running
ctb-tile dem.tif

I get the error:

Error: could not open GDAL datasetERROR 4: `dim.tif' does not exist in the file system,
and is not recognised as a supported dataset name.
Error: could not open GDAL datasetError: could not open GDAL dataset

I have also tried with the full directory name.

What seems to be the issue ?


texts: Using Calculate Field tool to calculate on null fields?

texts: I'm using ArcGIS 10.1.  I want to use the calculate field tool in ModelBuilder on a table where if field1 is null, it should take the value of field2.  If field1 has a value then it should stay as is.


texts: Can I put hyperlinks in Carto infowindows?

texts:
I'm mapping a bunch of nursing homes where there have been incidents of abuse. I'd like to put hyperinks to the incident reports in the infowindows. Those would be PDFs hosted on Google Drive.
Is this possible? And if so, how do the PDFs launch and display?


texts:
Is there a way to specify the coordinate type in Zonal geometry output in ArcGIS? I am trying to perform zonal geometry on a raster and in the resulting output, the XCENTROID and YCENTROID are in metres but I want them in decimal degrees.

The raster has a NAD_1983_Albers projection.


texts:
Did anyone have problems with Google Maps API? I have an HTML with OpenLayers and Geoserver. Since yesterday, the Google Maps Hybrid map is not working, even if some images are loaded in Firebug. Sometimes I receive a message stating Google Maps API v3.2 is deprecated.

 var map, info;

    function load() {
        map = new OpenLayers.Map({
    div: "map",allOverlays: true, projection: new OpenLayers.Projection("EPSG:900913"),
    displayProjection: new OpenLayers.Projection("EPSG:4326"),
    layers: [
       // new OpenLayers.Layer.Google("Google Satellite",{type: google.maps.MapTypeId.SATELLITE, numZoomLevels:50}),
        new OpenLayers.Layer.Google("Google Hybrid",{type: google.maps.MapTypeId.HYBRID, numZoomLevels:20}),
        cladiri=new OpenLayers.Layer.WMS(
                "Cladiri",
                "http://127.0.0.1:8080/geoserver/proprietati/wms",
                {layers: 'proprietati:<?php echo "".$identificator."_cladiri";?>',transparent: true, format: 'image/png'},
                {isBaseLayer: true, opacity: 0.9, singleTile: true})
    ],
    controls: [
        new OpenLayers.Control.Navigation({
            dragPanOptions: {
                enableKinetic: true
            }
        }),
        new OpenLayers.Control.PanZoom(),
        new OpenLayers.Control.Attribution()
    ],
    center: new OpenLayers.LonLat(26.1473586305531, 44.4431676053183).transform(new OpenLayers.Projection("EPSG:4326"), new OpenLayers.Projection("EPSG:900913")),
    zoom: 15
});


texts:
I've built out numerous network datasets before, so I'm confused as to why I'm getting this error.

I've built out a polyline file for the sidewalks in my community.

Warning: No route from location "Graphic Pick 1" to location "Graphic Pick 2".
Error: No solution found.

I can solve a route on a single line, but not on any route that requires taking multiple lines.

I checked all the options under 'location snap options', but that didn't solve it.

Figuring connectivity might be an issue, I merged the entire file into a single feature, but I still got the same error.

Does anyone have an idea of why it won't link up? Did I need to make sure the polygons only crossed each other on nodes so they intersected?

I don't remember running across this problem before.


texts: highlight feature with click in mapbox gl

texts:
I have a geojson layer of streets that are highlighted when moused-over.

My objective is now to highlight individual streets red with a click event.  Only one street should be able to be highlighted at a time, and should remain highlighted until a different street is clicked.

All of my attempts so far have failed. Any idea as to what needs to be added to the following code?

<!DOCTYPE html>

<meta charset='utf-8' />
<title>HTML markers from geoJSON url</title>
<meta name='viewport' content='initial-scale=1,maximum-scale=1,user-scalable=no' />

<link href='https://api.tiles.mapbox.com/mapbox-gl-js/v0.15.0/mapbox-gl.css' rel='stylesheet'/>


texts: Creating a polygon from all raster values >0 in QGIS?

texts:
I have a raster layer with two values, 0 and 1. I want to create a polygon shape from all the 1 values.
I expected the "raster to vector" tool to achieve this, but this created another grid as a shapefile (i.e. every square of the raster became its own polygon).
Am I using the correct tool and/or the correct settings?


texts: looping to process multiple files in the same folder using Python

texts:
I have a folder with filenames bio1_36.tif, bio1_37.tif, bio2_36_tif, bio2_37.tif, and so on. There are 19 bio variables and 36 and 37 at the end refer to the tile numbers. My goal is to mosaic the two tiles for each variable, clip to the size of a polygon layer called “Katanga” and convert the final output to asci format. I’ve been using the code below however this will only generate the correct output if the two tiles for the same bio variable are specified in the code or these two files are the only input in the workspace folder.

How can I do the above procedure in a recursive manner for each bio variable, i.e. setting up a loop?

This way I don’t have to create separate folders for each bio variable or specify them in the code.

I am new to Python.

Below is the code I’ve been using (note: "test" folder only contains one bio variable with the two tiles):

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

start_time = time.clock()
env.workspace = "G:/Eck_health_disease/spatial data/WorldClimate data/test"
arcpy.env.overwriteOutput = True
ImgList = arcpy.ListRasters()
print "Processing Mosaic"
print "Number of tiles to be Mosaicked:" + str(len(ImgList))
arcpy.MosaicToNewRaster_management(ImgList, env.workspace, "bio1_mos.tif",    pixel_type="32_BIT_FLOAT", cellsize="", number_of_bands="1", mosaic_method="", mosaic_colormap_mode="MATCH")
Katanga_poly = "G:/Eck_health_disease/spatial data/COD_adm_shp/DRC_Katanga.shp"
arcpy.Clip_management("bio1_mos.tif","21.7447528839 -13.4556760785 30.7780862172 -4.99734274513", "bio1_mos_Kat.tif", Katanga_poly, "#", "ClippingGeometry", "MAINTAIN_EXTENT")
arcpy.RasterToASCII_conversion("bio1_mos_Kat.tif", "bio1_mos_Kat.asc")
print "Task Completed!"
print time.clock() - start_time, "seconds"


texts: Cannot compute area of shapefile in field calculator of QGIS?

texts:
This shapefile represents the counties of NM. I am trying to compute the area of each polygon (yes, I know there is an area field, but this is a simpler example of a problem I'm having which doesn't have the areas).

I load the shapefile into qgis and set the projection to UTM 13 NAD 83.
I open the attribute table, enable editing, open the field calculator, select new field, real number, 20 digits, 4 digits precision, and set the expression to $area.

Of the 33 counties, all but 3 show an area of 0. The 3 with non-zero values have values that bear no relationship to the actual area.

I also tried reprojecting the shapefile to Conus Albers, open that in a project with the Conus Albers CRS, and get the same result. Obviously I'm doing something very wrong, but I'm not seeing it.

qgis 2.02 64 bit
Windows 7 Pro, 64 bit


texts: Determine whether a polygon is completely surrounded using QGIS

texts:
I have a shapefile dataset of polygons. I would like to find those polygons that are completely surrounded/enclosed. Mathematically, this means that the perimeter of Polygon A should equal the sum of the line segments of its neighboring polygons B,C,D,E,...n+1. If they are NOT completely enclosed, then the perimeter of Polygon A /= the sum of intersecting line segments of neighboring polygons B,C,D,E...n+1.
The polygons are irregularly shaped and not always square, so there aren't always the same # of polygon neighbors.
I found this post discussing how to do this in SQL, but I'd like to know how to write something similar in PyQGIS 3. (Determine if a polygon is not enclosed by other polygons).
I'm new at Python programming, but I imagine I basically need to:

import pyQGIS libraries
iface select active layer
calculate new perimeter field (data is in planar proj)
loop for every polygon in active layer:
*select polygons that intersect
*sum intersecting line lengths
*create new field of summed intersection lengths
*select where the two fields are equal

Is there a way to do this in QGIS? Would I need to convert to lines first?
Example of what I want to select for based on this approach:


texts: How to use a plugin in the QGIS 3.4 toolbox

texts: I have version 3.4 of QGIS and I want to generate a model, so I need SCP tools (unsupervised classification). I can not find a way to add it to the tool box in order to finish my working model.


texts: ArcGISDynamicMapServiceLayer change layer hierarchy

texts:
I am building an application off an existing map server using the javaScript API.  The server has about 10 layers that I am loading into my application as a ArcGISDynamicMapServiceLayer.

Once I have loaded the layers into my application, is it possible to re-order their hierarchy?  For example move the roads on top of the streams or move the parcels under the roads.

I have already found that I cannot change the transparency of just one layer, and can only adjust for all layers at once.  My guess is that it is not possible to reorder the layers from the browser.

If this is the case how should I suggest my clients reconfigure their map server?  Should each layer they want in the application be an individual mapService?  Any other suggestions?


texts: Importing AutoCAD 2017 file into QGIS/Google Earth

texts:
I have a map available in specific coordinates (2D) in AutoCAD 2017.

How do I import it into QGIS 3.8.3/Google Earth for data analysis?


texts: listing existing joins in table of contents feature

texts:
I am looking for a way to use python in ArcGIS desktop to check and print out the joins that exist between tables and feature classes in my map rather than right-clicking, going to joins and relates and seeing what is joined to the feature classes.

Any idea on a python script that would print the joins that exist?  The rest of the coding I think I could figure out.


texts: How do you access the Draw toolbar symbol(s) from ArcObjects?

texts:
I'm dropping a text element onto a map and I'd like to use whatever settings the Draw toolbar has (font, size, color, etc).

Is there a convenient way to access symbol interfaces?  Or do you have to create those from the individual settings?


texts:
I'm developing a GIS solution using Postgis database and geoserver as web map server. The solution to develop will analyzes the defects (as a polygon) and process parameters values (as a linestring) on the surface of steel sheets.

The number of the collected process parameters is about 100.

Considering the process parameters, which is the best way to store them considering that I have to apply some filter on row data before the analysis stage?

Their values as they are (number) building at run time the linestring for the layer applying functions ST_makeline(ST_makepoint(...)). In this way I can filter raw data before build the points and linestring.
Their values as a points over the length of the coil point(0, parVal) point(1, parVal) and then calling the ST_Makeline function. In this way I can filter points using ST_X in where clause before applying the ST_makeline() function.
Store directly linestring. In this way I can delete some points from the already build line string. (Is this possible?)

Of course there is difference in the execution time of the query among the three solution.


texts: GeoJSON EPSG::3857 to normal format (example normal format for Mapbox studio)

texts:
I'm really a novice in the field of the Geographic Information Systems (not to say that I don't know anything about it), but I retrieve via a WFS service some GeoJSON data.
The problem is when I want to import these data into Mapbox dataset, I get an error saying Input failed. Datasets don't support features outside longitude +/-180, latitude +/-90. (which I thought was normal before I realized that there are several "standards" if I understood correctly). In Mapbox tilesets I  had the following message old-style crs. Looking at the json file I see: crs:EPSG::3857.
So my question is how to convert these GeoJSON data into a more "conventional" format, especially to be able to import them into Mapbox data.
PS here a sample of my json :
{
    "type": "FeatureCollection",
    "features": [{
        "type": "Feature",
        "id": "id1",
        "geometry": {
            "type": "MultiPolygon",
            "coordinates": [
                [
                    [
                        [-178543.91187612, 5882448.58513225],
                        [-178577.52460631, 5882357.92829031],
                        [-178596.91742747, 5882324.55133574],
                        [-178626.66849743, 5882300.02122665],
                        [-178661.20199806, 5882290.4923541],
                        [-178543.91187612, 5882448.58513225]
                    ]
                ]
            ]
        },
        "geometry_name": "the_geom",
        "properties": {
            "limite": "50m *",
            "bbox": [-178661.20199806, 5882290.4923541, -178543.91187612, 5882448.58513225]
        }
    }, {
        "type": "Feature",
        "id": "id2",
        "geometry": {
            "type": "MultiPolygon",
            "coordinates": [
                [
                    [
                        [-347809.38958008, 6020388.58369802],
                        [-348032.50431796, 6020261.6016117],
                        [-347957.11952919, 6020410.60267021],
                        [-347921.2199014, 6020429.137701],
                        [-347881.08191941, 6020433.52972921],
                        [-347842.04256863, 6020423.19474383],
                        [-347809.38958008, 6020388.58369802]
                    ]
                ]
            ]
        },
        "geometry_name": "the_geom",
        "properties": {
            "limite": "50m *",
            "bbox": [-348032.50431796, 6020261.6016117, -347809.38958008, 6020433.52972921]
        }
    }],
    "totalFeatures": 65908,
    "numberMatched": 65908,
    "numberReturned": 5,
    "timeStamp": "2020-12-28T23:54:17.923Z",
    "crs": {
        "type": "name",
        "properties": {
            "name": "urn:ogc:def:crs:EPSG::3857"
        }
    },
    "bbox": [-350088.63932757, 5882290.4923541, -178543.91187612, 6021659.01574107]
}


texts:
I have clusters containing population counts in a mapbox-gl map. The cluster data is based on geohash aggregations from elasticsearch. I am having trouble sizing the diameter of my clusters. They should be based on the population count within the cluster as well as zoom level. What is the recommended algorithm for sizing clusters?

Here is what I have...

circleWidth = 10 * Math.log10(docCount + 1) * this.zoomLevel + "px";

where docCount is the number of people in the cluster and this.zoom level ranges from 0 to 18 where 0 displays the entire world and 18 shows just one or 2 city blocks.


texts: Temperature using Sentinel 3 in Google Earth Engine

texts:
I am trying to find information on how to obtain surface temperature layers - LST - from Sentinel 3 on Google Earth Engine but I am not getting anywhere. Is there any method to do it?
The Sentinel 3 bands are:

If this is not possible, which platform do you recommend for the LST calculation?


texts: No matter which of the basemaps I use in QGIS (Open Map, Google Maps, Google Hybrid, Google Terrain & Streets) when I export as image or pdf the basemap is so faint that all the features are unreadable, even though in looks perfect in Layout View. Just want to get WYSIWYG.


texts: How to use the linear referencing plugin from QGIS with postGIS data?

texts:
I need a guide how to use this plugin properly, I'm a newbie on PostGIS and I'm trying to use the linear referencing plugin of QGIS. I've created a database on postgreSQL with the postGIS extension and on QGIS with the Spit plugin I created a new connection with that DB. But, it generated some problems when it tries to connect, it says: fe_sendauth: no password supplied. But, I decided to save the username and the password when I made the connection and that problem it isn't showed.

So, in the next step I tried the plugin, with the configure tool, it show Connection OK, metadata OK. Then I go to Calibrate route, and as it's showed

It stops the calculation on the fourth step, and I don't know why. What do I need to have on each Attribute table to this works properly, or what do I do wrong?


texts:
Have iMac with OS 10.7.5. This is the first time I am installing QGIS. Everything was downloaded from kyngchaos.

This is what I did in that order:
Downloaded QGIS-2.12.0-1. Tried to install.
It asked for GDAL 1.11. Got GDAL_complete 1.11(EB). Installed.
Next QGIS asked for matplotlib. Got matplotlib-1.4.3. Installed.
Next QGIS asked for PROJ_framework 4.9. Got PROJ_Framework-4.9.2-2. Installed.
Next QGIS asked for SQLite3_Framework. Downloaded SQLite3_Framework-3.8.7.3-1. Installed.
At this point QGIS installed.

When I click on Q in applications, it says 'QGIS quit unexpectedly'.
ERROR:EXC_BAD_ACCESS (SIGSEGV)

Clearly something is not as it should be.
What do I need to do?


texts: Can't get dataProvider object outside of QGIS python Interpreter

texts:
Within PyCharm

import qgis.core as qgis
vector_lyr = qgis.QgsVectorLayer("Point", "distance nodes", "memory")
pr = vector_lyr.dataProvider()
print pr
prints None

Within Python Interpreter in QGIS

import qgis.core as qgis
vector_lyr = qgis.QgsVectorLayer("Point", "distance nodes", "memory")
pr = vector_lyr.dataProvider()
print pr
prints <qgis.core.QgsVectorDataProvider object at 0x000000000DAFF0D0>

I've tried this Getting dataProvider from vector layer outside QGIS and still no luck. Thanks.

[Edit] I found this: Access Sextante (processing) in standalone QGIS app? is this the case for dataProviders? Is there a workaround?

[EDIT2]

Tried this to no avail. I'm trying to create an in memory vector layer, write values to it then do something else after.

import qgis.core as qgis

qgis.QgsApplication.initQgis()

vector_lyr = qgis.QgsVectorLayer("Point", "distance nodes", "memory")
pr = vector_lyr.dataProvider()
print pr
prints None


texts:
I need to summarize a set of incident data by several different stations. I have a set of five stations. The stations code emergency response using 10 incident codes, 1 being least severe and 10 being most.

I need to summarize the data so that each station has its own incident count, e.g., Station 1 had 34 incident 1 calls, 35 incident 2 calls, etc. Station 2 had 34 incident 1 calls, etc.

Essentially, I need to have a table output that has 5 rows and 10 sub-rows per main row. Is this possible?


texts: How to select from geometry column in postgresql?

texts:
I've got a table which it has a geometry column filed by some point (ex. "POINT Z (116.309627 39.991285 148)").
now I want to select point's id from the table. I wrote follow sql Code but I force an error:

  ERROR: Operation on two GEOMETRIES with different SRIDs SQL state:
  XX000

select id from database
    where (point_lat_lon = 'POINT Z (116.309627 39.991285 148)');


texts: QGIS is not working if the Python script is running with infinite loop

texts:
While executing the following code, I'm getting the error in the QGIS Software window. Err: "QGIS Desktop is not responding"
If the same code executed in the standalone Python application and added a few new points in the QGIS layer then the layer_details.dataProvider().featureCount() value is not changing. It is not considering the new points.
import time
while (True):
  layer_details = QgsProject.instance().mapLayersByName(layer.name())[0]
  print (layer_details.dataProvider().featureCount())
  time.sleep(10)


texts: Carto - problem with data import in carto format

texts: I want to do some CARTO guides (https://carto.com/learn/guides) but I can't import the data. It is already a carto format and it doesn't work when I try to drag&drop the file. Any ideas?


texts:
So I'm pretty sure I'm missing something very simple, but --

I'm trying to create a new dataset inside CARTO, by using the "Connect Dataset" button and choosing "Box" (I want to connect to a CSV file inside my Box account). Basically the exact workflow described here in the "Connect Dataset" section.

Unfortunately the list of files I get to choose from in Box doesn't seem to include my file? The interface says it's displaying 200 files -- a round number like that makes me suspect there's some kind of limit to how many files can be read from Box & it's somehow arbitrarily picking which ones it's making available to me? I have way more than 200 files in my Box account. The "Search" function doesn't seem to work at all, as it turns up no results for anything, no matter what I search for.

I went through this process perfectly once before, but now the list simply doesn't include the file I need to import/sync to -- though it's an accepted format. How can I make the CARTO import mechanism via Box see the CSV file I want?


texts: Is there a g.projector overlay with all world's administrative regions?

texts:
I'm searching for an overlay for g.projector including all (or most) of the world's administrative regions. (apparently, valid formats are CNO, CNOB and SHP)

With "administrative region" I mean the administrative level just below the nation level, such as US States, UK Regions, etc.

I've found some nation-specific ones at statsilk, but they aren't many, and they are a single file for each nation, so I wouldn't even know how to import it in g.projector…

How can I find such a file, assuming it exists?

Consider that I'm pretty new at this, so it's possible I'm missing something obvious. I've searched on google, following many links, browsing many sites, and noticed there are many different maps… only I couldn't find that one I'm looking for.


texts: High-level method for offline future road elevation data

texts:
I am new to GIS, so please do not hesitate to call me out in providing limited information.

I have a .tiff Terrain (elevation) data file for the east coast of the US and a .osm.pbf OSM file for North America. My objective in Python is as follows:

Using a gps coordinate, get the closest road a vehicle could be on (using the .osm.pbf OSM data?)
Get the gps coordinates for this road for the next and previous 5 miles
Finally, use these future road gps coordinates to get the estimated elevation data for the road (using the .tiff Terrain data?).

Requirements: Only offline use. (Ex. Can't use Google snap to road APIs)

I've been playing with the .tiff in QGIS and using the profile tool to draw polylines and seeing the elevation over the line.  I've also been reading up on programming with it using Python.

I'm looking for some advice on completing #1 and #2 and if it is doable with QGIS (is there QGIS compatible road data?) or will I have to interface with a pbf analyzer like "imposm.parser"?

There is a similar question, but it deals with different file types and is geared more towards step 3: Extract Elevation of Points from DEM using Python QGIS or GRASS


texts:
I'm facing a big issue...i moved from the ubuntugis version of grass to the svn one.

I managed to compile the svn (releasebranche6.4) version nicely but i can't  download the gdal-grass plugin for the latest version of gdal (the ubuntugis version) 1.8 because remotesensing.org/pub/gdal or download.osgeo.org doesn't work at least in my country (Sénégal).

I use ubuntu 11.04 natty, grass 6.4.2svn and gdal 1.8.0-2

So what can i do to install this plugin ?

i already scan many website during 3 days without success....any clue ?

PS : sorry for my basic english


texts: ArcGIS 10: Create an event linear layer from a point layer based on directionality

texts:
I have a point layer that reflects the speed limits and a line layer of the roads. The location of the speed sign indicates in which direction the speed limit applies.

How can I make linear event table on top of the road layer that reflects the speeds? So for each segment, return two speed attributes, one for each direction.


texts: Is there a world projection that looks like a 3D globe?

texts:
I would like to try and show data such that it appears as if it was on google earth. For example, in this case, the view of the earth is circular-like and the data shown curves around. Are there any map projections that enable me to do this?

(I do understand that geographic projections convert information from a 3D space into a 2D space. What I want to do, just for illustration purposes is to make the 2D space look like 3D!)


texts: How to replace accented strings in attribute tables using Python or VBA ArcGIS 10.2

texts:
I need to replace accented words, into the Field Calculator, to the same non-accented ones.

For example: María  ---> Maria

Need some help with this VB Code:

Dim data as integer
Dim contador
Dim words string

word= [POLYGON_Name]

Dim v_matriz(5)
Dim n_matriz(5)

v_matriz(0) = "Á"
v_matriz(1) = "É"
v_matriz(2) = "Í"
v_matriz(3) = "Ó"
v_matriz(4) = "Ú"

n_matriz(0) = "A"
n_matriz(1) = "E"
n_matriz(2) = "I"
n_matriz(3) = "O"
n_matriz(4) = "U"

for each i in v_matriz
    data=InStr(word,v_matriz(contador))
    if dato>0 then
        word=Replace(word,v_matriz(contador),n_matriz(contador))
    end if
    contador=contador + 1
next

But I'm having problems, I'm new in this scripting style, and maybe I'm doing something wrong.

I could use Python also, if you want.

SOLVED:
This code in Pre-Logic Script Code:

import unicodedata
def strip_accents(s):
   return ''.join(c for c in unicodedata.normalize('NFD', s)
                  if unicodedata.category(c) != 'Mn')

This in VALUE=

strip_accents(***your accented column***)

Thanks for all!


texts: Error with popup overlay in OpenLayers 3?

texts:
In OpenLayers (OL) 3.15.0 I declare a popup contains like :

popup = $(eval('

My contain

'));

And after, I declare my overly with this popup element :

var overlay_popup = new ol.Overlay({
                        element: popup
                    });

On OL 3.7.0 I've any problem all was okay.

But now with this new version of OL (3.15.0) Firebug say :

  TypeError: Argument 1 of Node.appendChild does not implement interface
  Node.

What can I do for working ?


texts: How to add a dynamic number of pictures into a Leaflet popup?

texts:
I have a huge dataset of playgrounds I want to represent in a Leaflet map.

Each playground has a different number of photos (between 1 and 13) connected as an attribute value.

The value contains the relative path to each photo, the structure could be changed to e.g. a comma-separated list or anything else that could help. The popup content as an HTML-table should place one photo in each row.

How do I manage to adapt the showing table rows inside the popup dynamically, depending on the number of photos?


texts: Is there a Select Feature by Attribute tool for Feature Classes when performing selection ModelBuilder?

texts: I want to perform a Select by Attribute selection on a Feature Class in ModelBuilder, but I know the Select by Attribute tool does not work in this situation.  I need to perform this function as I want to avoid creating an additional dataset.  This task is the second to last one, the last one being a Delete Rows.


texts: Seeking Android app to save offline map and POIs?

texts:
I don't know much about GIS and would like to find an application for travelling with the following features:

support Points of Interests (sights, restaurants, etc.); Clicking on any POI should display a pop-up with a bit of information (eg. address, history of a building, etc.)
download maps for offline use, so that no data connection is required
should display the user's current location using GPS. Apparently, a data connection is not required to use GPS, although it might speed things up when enabling GPS

Are there Android applications that can do this?

That would make it easier to travel rather than carrying infos on paper.

Android Google Maps doesn't support layers for offline maps, so adding one's locations through the web Google Maps won't work.

Users recommended MapsWithMe which is supposed to support "bookmarks" in the Pro version. If someone uses the Pro version, can you confirm that locations are downloaded and saved in offline maps?

Also, is it possible to build a list on a PC through eg. Google Maps/MapsEngine, and tell MapsWithMe Pro to download it?

E-mails to the author (info@, support@, and android@) returned no reply.


texts: Using QGIS, automatically create n polygons fully covering an area without overlap with an approximately equal number of points in each polygon

texts:
I'm trying to avoid/simplify the following manual process intended to create a polygons layer with exactly n polygons so that each polygon encloses the same or close to the same number of points in an existing points layer:

Select the approximately correct number of points and create a convex hull n times to match the required number of polygons.
Use v.clean\rmdupl to create separate polygons for overlapping areas.
Use Eliminate Sliver Polygon to get touching but not overlapping coverage with n polygons.
Use Points in Polygon to determine if the distribution is close to equal.
Manually move polygon boundaries and rerun the counts as many times as needed to arrive at the solution.

The polygons do not need to be regularly shaped.


texts: Trying to export stacked text to CAD using ArcGIS 10.2

texts:
I've been experiencing some unwanted results when trying to export stacked text to CAD using ArcGIS 10.2 for desktop. With the Export to CAD tool I can easily export polygons. However, I'm also trying to export the 'stacked text' labels inside those polygons.

With the Feature to Point tool I then export a 'CADtype' field that was created with the 'Add Field' routine and populated with values from other fields using the 'Field Calculator' tool.
My results are the correct text values exported to the centroid of each polygon but no longer stacked in columns, and a "?" at the end of each line. There are no error or warning messages generated.

I've tried both VBScript and Python expressions for the Field Calculator syntax but I'm getting the same results. I can't use the ArcGIS for AutoCAD plugin because the CAD on my machine is not current enough and I don't have FME.


texts: Looping through 2 folders to perform extract by mask in ArcPy?

texts:
I am trying to perform the extract by mask function on rasters stored in workspace1 using the rasters stored in works. Pairs of rasters have the same name.

outdir="E:\output"
workspace1="E:\workspace1"
arcpy.env.workspace=workspace1
inras_list=arcpy.ListRasters("*",".tif")

workspace2="E:\workspace2"
for inras in inras_list:
    rasx=arcpy.sa.Raster(inras)
    rasy=arcpy.sa.Raster(os.path.join(workspace2, inras))
    inras_extent=rasx.extent
    print inras_extent
    maskras_extent=rasy.extent
    print maskras_extent
    out_extract= arcpy.sa.ExtractByMask(rasx,rasy)
    outname = os.path.join(outdir, os.path.basename(rasx).split(".")[0] +
    "_opporwithX.tif")
    out_extract.save(outname)

The code works with no errors, but I see no outputs.

What is wrong with the code above and how do I fix it?


texts:
I have OSM data (2D) data from Geofabrik now i need 3D building data for my study. I found one web server where i got OSM 3D buildings
https://osmbuildings.org. Now i need to either crawl this data and convert it into json file which i can use it in QGIS. Anyone having idea how i can get this kind of 3D building data which i can fetch it in QGIS.


texts:
I am working on the data extract widget for the flex 2.5 viewer. It has installed properly and works fine but the output is not what I expected. The tool will not divide or split polygons. For example when I draw the AOI the return will omit any polygons that have been split or divided. I was expecting to get all data in the AOI. Is this the way the tool works? Or is there something I can do to fix this? My .xml file for the widget is below:

-

-->
-
http://s2-cor-pangaea:8399/arcgis/rest/services/demo1/GPServer/Extract%20Data%20Task
draw


texts:
I am very new to QGIS.

I have read some similar questions on here but none seem to get me to the correct result. I have a bit of experience with Global Mapper but very little with other GIS software or code so this is new territory for me.

I'm looking to create a hexbin map of point data (from a csv), where each hexbin represents the most common attribute (eg NAME, or a unique integer for each name) of the points within it. So within one hex there may be 500 1s, 200 2s, 3 26s, and 1 63. The output I want would be just 1.

Something like this:

To this:

Except in QGIS.


texts: Create a plugin to connect PostgreSQL db and run a query using QGIS Plugin

texts:
I have prepared a code to connect PostgreSQL database and run a query on the database using PyQGIS and it is running successfully on python console of QGIS. I want the same to happen using QGIS Python Plugin.
How can I write that in a Python Plugin File or inside which method shall I write the same.


texts:
I am using the ModelBuilder to apply some analysis on some feature classes(shape files). I have to apply the same model for different feature classes and of course the output naming would be according to the input data naming.I would like to know whether it is possible to define a loop in ModelBuilder to apply the model over a group of shapefiles and differentiate it based on the input data name.

I have heard that there is a possibility to use ArcPy to define loop but as I am not a pro in ArcPy , would rather to do it in ModelBuilder.


texts: Auto-complete function does not work In the QGIS python editor

texts:
I use QGIS 2.14.3-Essen. The auto-complete function works without any problems in the python console, but it does not work in the QGIS python editor. I also tried it on another PC running QGIS 2.14.3 - essen-based, both of which worked well. I even tried to re-install QGIS, but still was not able to solve this problem.

How do I solve this problem?


texts:
I'm trying to intersect two polyline feature classes and generate a point output using geometry objects. Its intersecting almost 2 million features with 250,000+ features, so running it with arcpy.Intersect_analysis has been running for...probably a week at this point with no progress. So I'm looking for a faster way. I'm getting AttributeError: 'list object has no attribute 'intersect', this occurs on line 13 at intersect = ln_geom.intersect(fn_geom, 1). All the documentation I've found says that is how to do it.

import os
import arcpy

Fishnet = r"C:\Users\AO_Fishnet_Testing"
Workspace = r"C:\Users\Testing.gdb"
Lines = r"C:\Users\Testing.gdb\AO_Testing_Ln"
Output = r"C:\Users\Testing.gdb\Testing_Points"

ln_geom = arcpy.CopyFeatures_management(Lines, arcpy.Geometry())
print "Lines Copied"
fn_geom = arcpy.CopyFeatures_management(Fishnet, arcpy.Geometry())
print "Fishnet Copied"
intersect = ln_geom.intersect(fn_geom, 1)
print "Intersected points"
arcpy.CopyFeatures_management(intersect, Output)


texts: Project an unknown coordinate system on a map

texts:
I am currently working with a data set, which uses an unknown coordinate system from the south island of New Zealand, using x and y coordinates. I have these in an Excel sheet, but no further information on the coordinate system is known.

I have used Excel calculations and GPS approximations from Google Maps to translate the unknown coordinate system to the NZTM one, but the results still show a deviation of about 5km for some points. As I need to have very precise positions, I am thinking about using ArcMap to project the data points on a New Zealand map, adjust the points so it fits to the map and then use some kind of translation to get from the data points to e.g. GPS coordinates.

I've tried to illustrate that with the pictures bellow.

Picture one: The blue dots represent the x,y coordinates of the unknown coordinate system. Is it possible to insert those data points into ArcMap like shown on the picture?

Picture two: The second step would be to adjust the data points s.t. it fits the New Zealand map. Is such an operation possible with ArcMap?

The last thing would be to use exchange the initial x,y coordinates with e.g. GPS coordinates.


texts:
My question is similar to this one Get Raster Values from a Polygon Overlay in Opensource GIS Solutions but I think I need another step.

I have a polygon layer of ecoregions which I brought into R with "readOGR". It has many attributes such as ecoregion and biome (biomes encompass multiple ecoregion polygons). I have a raster layer (I used "raster" on a tif) that is 0-2. I also have continuous raster layers, that I'll use later so I'm looking for a generalized solution.

I would like to be able to do a variety of evaluations like total area and proportion of each biome >0 or when I use the continuous raster, things like mean or range. Ideally I'd like to end up with some sort of attribute table that I can work with that combines the polygon attributes with the raster values. I'm not sure if it's most accurate & efficient to convert both to polygons, both to rasters, or to do summaries as is.

I am trying to work in R so that I can script the process; sample code is particularly helpful as I'm a novice.
I appreciate any thoughts.
Thanks a lot


texts: QGIS, Print Composer's legend regresses to old style classes

texts:
I'm using QGIS 2.16 to create some rather simple vector maps.
However, ever since I made some modifications to the style of one layer (e.g. I added one more class) the Print Composer regresses and gives me the old style classes as soon as I deselect the Auto update function. I already tried reinstalling the software (I had run 2.15 before) and saving/loading the new layer style, but to no avail.
Any advice on this?


texts:
According to Wikipedia,

  The predicates Contains and Within have subtle aspects to their definition which are contrary to intuition. For example,[10] a line L which is completely contained in the boundary of a polygon P is not considered to be contained in P. This quirk can be expressed as "Polygons do not contain their boundary". This issue is caused by the final clause of the Contains definition above: "at least one point of the interior of B lies in the interior of A". For this case, the predicate Covers has more intuitive semantics (see definition), avoiding boundary considerations.

If Polygons do not Contain their boundaries, do linestrings ever "contain" points?


texts: Can anyone point me to good development environment in Linux for client side search with Flex or JavaScript?

texts:
I want to create a client in MapFish, and/or opengeo, using OpenLayers, maybe Cartaro and would like to play in a development environment using Flex or JavaScript.

I need to take my data, make it searchable from table columns and features, and be able to search strings for geolocation content, any good tuts, links, or Linux development environments that you could recommend?


texts: Struggling to use Field Calculator - PyQGIS

texts:
I'm using QGIS 3.4 (Madeira) and I am struggling to use python to calculate attributes of a new field in the attribute table.

I have already joined a csv and a shapefile. I use the code below to create a new column in the attribute table called 'PatCNT' with type as 'Int'

join_layer = QgsVectorLayer('/Users/ep9k/Desktop/SandraMonson/cb_2017_us_zcta510_500k/cb_2017_us_zcta510_500k.shp', 'US Zip Codes', 'ogr')

caps = join_layer.dataProvider().capabilities()
if caps & QgsVectorDataProvider.AddAttributes:
    join_layer.dataProvider().addAttributes([QgsField('PatCNT', QVariant.Double)])

Now, I have another column in my attribute table called Patient_Data_PatientCount (which is what I joined from the csv file).
I want to copy that column as an integer. When imported from a csv, it's type is 'String'

I've made some attempts based on other Stack Overflow threads like:

-Is it possible to programmatically add calculated fields?

-How to use QGIS field calculator in python?

-Cannot calculate new values for an empty attribute using existing attributes in QGIS with python
However these are all for QGIS 2.18 and there are some differences (pendingFields doesn't exist anymore)

I've also been looking at the documentation and in the PyQGIS developers cookbook with no success

-https://qgis.org/pyqgis/master/core/Expression/QgsExpression.html

-https://docs.qgis.org/testing/en/docs/pyqgis_developer_cookbook/vector.html


texts: disable M value and Z value by using arcpy

texts:
I want to remove the M and Z values of the shapefiles, which can be done by using Feature Class to Geodatabase (multiple) and choosing Disabled for M and Z values. The following picture shows the details of the step.
I want to use arcpy to realize this step so I used Copy As Python Snippet and obtained the below code.
arcpy.FeatureClassToGeodatabase_conversion(Input_Features="LinearRef_Travis_LS;Merged_Final", Output_Geodatabase="C:/Users/TL/Desktop/1031/New File Geodatabase.gdb")

The code doesn't show the part of changing the environment setting. I wonder if there is way to remove the M and Z values of shapefiles by using arcpy?
I am using ArcGIS Desktop 10.3.1 Advanced License.


texts:
I'm curious about how -thin_with_grid is actually performing thinning tasks within lidR.
Which point is selected? One point within the grid is randomly sampled? Is it possible to select for example the lowest?
I know there is a full list of filter parameters with lidR::readLAS(filter = "-help") but I think this part of documentation is relatively poor (e.g. clear info about input parameters or how to combine different filters)


texts: how can i convert a shapefile to lat and lon boundaries?

texts: i have a shapefile of country subdivisions  and i would like to extract an array of bounding lat and lon for each division.. is it possible to do?


texts: How to understand degrees in real estate

texts:
We are wanting to understand where the easement is on a 10 acre property we purchased. We don't really understand about tags for posting.
This is the description:

  Ingress-egress easement for road purpose: Beginning 12.5 ft on each
  side of the following described center line: Beginning at the center
  of section 4, township 23, north range 32 west, McDonald County, MO.
  Thence north 68 degrees, 39 minutes, 25 seconds east 7l.96 feet.
  Thence south 75 degrees, 51 minutes, 17 seconds east, 55. 41 feet,
  thence south 81 degrees, 40 minutes, 50 seconds east, 162.27 feet to
  the point of termination.


texts: SQL Query not selecting all points that it should?

texts:
I am working on a map of 60+ monitoring wells and I plotted them through an excel data sheet and a few tools. In that data sheet, I have Benzene concentrations associated with almost each point. I need to have the ones that are above SHS to be green and the ones that are below SHS to be red. I have tried setting up SQL Query a few times and also tried selecting by attribute, but every time it only selects one or two points when it should be selecting 12 points (7, 17, 23, 25, 37, 58, 140, 180, 230,  270, 280, 1100), I've tried a few different queries including "Benzene" > '5', which selected 2 points (7 and 58), "Benzene" >= '7', which selected one point (7). I tried selecting by attribute and only 7 was selected every time. Also, when I've tried doing the opposite in label editor ( "Benzene" <= '5', it says it is a valid expression, but all the numbers above 5 are still showing up green).

What am I doing wrong?

My benzene column is a String.


texts: Load geographical data administrative boundaries for Europe by levels

texts:
How can I load all administrative boundaries by levels for Europe in coordinates?

I need .osm dates, but not all data inside Europe, just administrative boundaries. I import .osm for Europe in PostGIS base and by filter 'administrative' read data, but I am not sure is it right query.

Can I read those data from natural-earth and how?


texts: DotSpatial problems with reprojecting shapefiles

texts:
I have a shapefile with this projection:

PROJCS["RD_New",GEOGCS["GCS_Amersfoort",DATUM["D_Amersfoort",SPHEROID["Bessel_1841",6377397.155,299.1528128]],PRIMEM["Greenwich",0],UNIT["Degree",0.0174532925199432955]],PROJECTION["Double_Stereographic"],PARAMETER["False_Easting",155000],PARAMETER["False_Northing",463000],PARAMETER["Central_Meridian",5.38763888888889],PARAMETER["Scale_Factor",0.9999079],PARAMETER["Latitude_Of_Origin",52.15616055555555],UNIT["Meter",1]]

Now I want to reproject the shapefile to lat/lng coordinates.

I tried:

Shapefile tmpShape = Shapefile.OpenFile(file);
tmpShape.Reproject(KnownCoordinateSystems.Geographic.World.WGS1984);

The lat/lng area of the shapefile should be something like 52,.../4,... .
But after reprojecting the file with DotSpatial, I got something like -37,.../23,...

What did I do wrong?


texts: How do I create an IFeature object from an IPolyline or ILine object

texts:
I am having some problems with some code i am trying to write. The application i am creating allows the user to click on the map and a point will be createD, when they click again to create another point a line is created between the two points.

The code i have got at the moment is

for (int i = 0; i < pointCollection.PointCount - 1; i++)
{
    //create the line from the points
    ILine pLine = new Line();
    pLine.PutCoords(pointCollection.Point[i], pointCollection.Point[i + 1]);

    IPolyline newPolyline = new PolylineClass();
    newPolyline.FromPoint = pLine.FromPoint;
    newPolyline.ToPoint = pLine.ToPoint;

    //create a feature ****STUCK HERE****
    IFeature newline = new PolylineClass() as IFeature; ;
    newline.Shape = newPolyline;

    geometriesToUnion.AddGeometry(newline.Shape as IGeometry);
}
resultPolyline.ConstructUnion(geometriesToUnion as IEnumGeometry);

Has anyone got any idea if this is possible or a another way of doing it.

Thanks


texts:
I am try to upgrade Addin for ArcGIS finish to install new upgrade but Arcgis run the old one I am finish to remove and install still I see version 1 not version 2

How can I fix this are any temp folder I can delete ?


texts: Openlayers 4 - get all child layers in a layer group

texts:
How can I get all child layers of a given layer group in Openlayers?

I go as far as establishing that a given layer is in fact a group, but then I get stumped:

// Handle visibility control
$('i').on('click', function() {
    var layername = $(this).closest('li').data('layerid');
    var layer = findBy(map.getLayerGroup(), 'name', layername);
    var flagGroup = layer instanceof ol.layer.Group;
    if (flagGroup) {
       ...
    }
});


texts: ArcGIS Workflow Manager - How to setup a query to display current step?

texts: I'd like to setup a query that shows the current step as one of the fields in the results window.  For example, we have a query that shows all jobs in the system.  The challenge for me is since "step" is not one of the fields in the main jobs table, it has to be retrieved thru an intermediary table(s).  I'm not sure how to create a query that would include this field in the display. Using v10.2.  Another way to ask my question is how do I display a field that is not in the main jobs table, in the results view.


texts: Specifying a variable in ArcMap output if this is not present in the input

texts:
I have been using ArcMap 10.6 (working from within a .gdb, with vector files) to buffer selected variables from a number of files (classes of water, buildings, roads and a constructed other for the rest) at different distances (from 100 to 500m) from a defined sampling point in a city (there will be 60 of these points – so 60 points with 5 buffers each). Each buffered sampling point will consist of differing %’s of these variables making up the buffered area. The areas, defining different urban planning indicators, will be used to assess correlation against temperature at each sampling point to see what relationships can be seen at differing distances.

My issue is that I’m exporting this data to Excel and using pivot tables (which I’m familiar with and use quite a bit) to further process the data (min and max, percentages, volumes, etc) and refine into a summary table for graphing. I have the process sorted, however the problem I am having when linking to the pivot table is that, if there is no “water” in the buffered area, this will not be included in the buffered area produced by ArcMap (obviously). This random loss of one of the variables causes issues when summarizing data from the pivot tables, as the ability to do an easy sequential copy and paste link is not available (sometimes it’s there, sometimes not). I can do this manually, checking each link, however this is very time-consuming.

My question is essentially an ArcMap one. Using vector files, can I instruct ArcMap to somehow do the following If/Then - “if there is no water within the buffered area, provide a water class and give it a value of 0”, so this is recorded in the exported data, and each sample site then has 4 discrete classes (roads, buildings, water and other). For example, sample point 1, buffered at 100m, class = water will have an area_sqm = 0 if no water in the buffered area. To get to the final, merged buffers is a multipart process. I am looking at using a lot more variables in future so an easy fix for a missing variable would be great. The process works fine, it’s just this hiccup with the Excel automation, coming off the base ArcMap data. I did a bit of a search but could not find anything related to this issue. Maybe possible in Model Builder or using Python both of which I'm not overly familiar with?


texts: Reading geospatial properties of ArcGIS web services using C#?

texts: I need to display the Geospatial properties of ArcGIS web-services (like service name, service type, description, copyright, etc.) which are available in the "ArcGIS Services Directory", in my C# application.


texts: Unable to call GDAL functionalities from Python

texts:
I am trying to call GDAL functionalities from the following script I made in the past.

import os, fnmatch
from subprocess import call
call("dir", shell=True)

# define the variables
inFolder= os.path.join ('D:/', 'UMBE', 'PhD', 'Guglielmin', 'Permafrost',

'Alta_Valtellina', 'Clip_AltaValtellina', 'clip_sept2000')
outFolder= os.path.join ('D:/', 'UMBE', 'PhD', 'Guglielmin', 'Permafrost', 'Alta_Valtellina', 'permafrost_mapping', 'spectral/')
clip= os.path.join ('D:/', 'UMBE', 'PhD', 'Guglielmin', 'Permafrost', 'Alta_Valtellina', 'permafrost_mapping', 'RG_visible_parts.shp')

# path to a folder where to look for needed files
os.chdir (inFolder)

# define a list of raster within a specific path
def findRasters (path, filter):
    for root, dirs, files in os.walk(path, filter):
        for file in fnmatch.filter(files, filter):
            yield os.path.join (root, file)

# for loop-clip all the raster with a poligonal shp
for raster in findRasters (inFolder, '*_refl.tif'):
    (infilepath, infilename)= os.path.split (raster)
    print infilename
    outRaster= outFolder+ 'activeRG_'+ infilename[5:]
    print outRaster
    warp= 'gdalwarp -dstnodata 0 -q -cutline %s -crop_to_cutline -of GTiff %s %s' % (clip, raster, outRaster)
    call (warp) # this command is used to call external functionalities (in this case GDAL)

I am basically trying to batch clip some raster images with the same feature (a polgonal shapefile) using gdalwarp. I am receiving the following error:

Traceback (most recent call last):
  File "D:\UMBE\PhD\Guglielmin\Permafrost\Alta_Valtellina\permafrost_mapping\spectral\AltaValtellina_batchclip.py", line 26, in <module>
    call (warp) # this command is used to call external functionalities (in this case GDAL)
  File "C:\Python27\ArcGIS10.1\lib\subprocess.py", line 493, in call
    return Popen(*popenargs, **kwargs).wait()
  File "C:\Python27\ArcGIS10.1\lib\subprocess.py", line 679, in __init__
    errread, errwrite)
  File "C:\Python27\ArcGIS10.1\lib\subprocess.py", line 893, in _execute_child
    startupinfo)
WindowsError: [Error 2] Impossibile trovare il file specificato

I think this is due to the lack of GDAL variables in the Windows path. Moreover, I am using the default Python shell provided along with the ArcGIS installation.
I tried to workaround this (possible) problem by starting the Pythonwin shell in the QGIS folder ('~\QGIS Brighton\apps\Python27\Lib\site-packages\pythonwin\pythonwin.exe'), but it does not open, saying "The application can not locate win32ui.pyd(or Python)(126). Impossible to find the specified module", although the "win32ui.pyd" lies in the same folder.
I also tried to add some variables to the system path (e.g. 'C:\Program Files\QGIS Brighton\apps\Python27\Lib\site-packages\osgeo\'; 'C:\Program Files\QGIS Brighton\'), and creating and adding new variables (like 'GDAL_DATA') as explained here, but nothing...


texts: Computing number of points in a raster grid cell in R

texts:
I have a country shapefile and a dataframe with coordinates of bird sightings. My goal is to "rasterize" the shapefile and compute the number of points within each grid cell. Specifically, I want to end up with a grid-cell level dataframe containing lat/long of the centroid of each grid-cell, the district census code of that grid cell (which is contained in the shapefile attribute table), and the number of points falling in the grid cell. I would like to keep cells with no points in them and leave number of points as "NA". This is what I have so far. I get stuck after finding the district code of each cell. I guess the next step is to convert the raster to a polygon grid and then overlay the points on it? Is there an easier way?

# Load Packages
require(data.table)
require(rgdal)
require(rgeos)
require(sp)
require(raster)

# Load india 2011 districts
india.districts.2011 <- readOGR(paste(dist_shpf, "maps/india-district", sep=""),
                            'SDE_DATA_IN_F7DSTRBND_2011', stringsAsFactors = F)
# Read Bird Data
df.sample <- fread('ebd_IN_200001_201810_relOct-2018.txt')

# Initialize Grid
grid <- raster(extent(india.districts.2011))
res(grid) <- .25
proj4string(grid) <- proj4string(india.districts.2011)
india.grid <- rasterToPolygons(grid) #SpatialPolygonsDataFrame

# Get centroids
india.grid.centroids = gCentroid(india.grid ,byid=TRUE)
india.grid.centroids$c_code_2011 <-
over(india.grid.centroids,india.districts.2011)$district.code #SpatialPointsDataFrame


texts: How can I get information displayed at tooltips for each algorithm's name at Processing Toolbox for QGIS 3?

texts:
I use following code for searching adequate id name for each algorithm presents in Processing Toolbox (gdal, grass7, qgis, saga) for QGIS 3.

from os import listdir
from os.path import isfile, join

gdal_path_algs = processing.algs.gdal.__path__[0]
grass_path_algs = processing.algs.grass7.__path__[0] + '/description'
qgis_path_algs = processing.algs.qgis.__path__[0]
saga_path_algs = processing.algs.saga.__path__[0] + '/description'

gdal_algs = [ 'gdal:' + file[:-3].lower() for file in listdir(gdal_path_algs) if isfile(join(gdal_path_algs, file)) ]
grass_algs = [ 'grass7:' + file[:-4].lower() for file in listdir(grass_path_algs) if isfile(join(grass_path_algs, file)) ]
qgis_algs = [ 'qgis:' + file[:-3].lower() for file in listdir(qgis_path_algs) if isfile(join(qgis_path_algs, file)) ]
saga_algs = [ 'saga:' + file[:-4].lower() for file in listdir(saga_path_algs) if isfile(join(saga_path_algs, file)) ]

#print (gdal_algs)
#print (grass_algs)
#print (qgis_algs)
print (saga_algs)

However, when list algorithms is printed at Python Console, for instance saga, some of them has special characters ("-", "(", ")") or author's name of algorithm and id's name doesn't work for retrieving parameters list by using 'algorithmHelp' processing method.

Afterwards, I founded out that tooltips displayed when I put cursor over algorithm's name at Processing Toolbox (see following image) had correct information I need. How can I get information displayed at tooltips for each algorithm's name at Processing Toolbox for QGIS 3 as alternative approach of my above code?


texts: What is the difference between loading a background file or a data file onto a Trimble unit?

texts:
I am curious what the difference is between loading a background file and a data file onto a Trimble in Terrasync, as both seem to be able to display a .shp file.

Is it that a data file can then be edited while in the field but a background file can not?


texts: Updating attributes between layers with a spatial join

texts:
In QGIS 3.x how to copy the attributes from a point layer to a polygon layer where the point is within the polygon? In PostGIS I would do:
UPDATE poly_layer SET poly_layer.attribute = point_layer.attribute
FROM point_layer WHERE ST_Within( point_layer.geom, poly_layer.geom );

I realize that there is the "Join Attributes By Location" tool but it creates a new layer which is not convenient for my workflow. And I realize that if my two layers were both PostGIS layers I could simply execute my query in the DB Manager. But in this case I'm trying to copy attributes from a DXF point layer into an existing PostGIS polygon layer.


texts: What is a good resource for CA crop production data?

texts: To construct maps depicting crop yield throughout California, are there statewide datasets, or would I need to add data by county?


texts:
I'm trying to create hillshades by importing data used by OsmAnd but I can't figure out how to do so.

When adding the data as a Vector Layer, as I normally would do with sqlite data, it appears in the Layers Panel but doesn't show anything.

OsmAnd Hillshades are in sqlite format. Technical help from OsmAnd states that SQLite format used in OsmAnd is based on the "BigPlanet" SQLite as supported by MOBAC.

Looking at the sqlite file format shows that the table has columns x, y, z, s and an image blob. Is there any way to specify sqlite import format for QGIS?

PS: This conversation seems relevant but I don't understand it fully.


texts: Is there a commandline tool for intersecting Shapefiles?

texts:
in QGIS vectorlayers can be intersected. That is apparently done by ftools and seems to be integrated into the QGIS-program itself.

That is different from many other functions that call GDAL-tools.

I would like to intersect esri-shapefiles in a script and without having to load anything into QGIS. Is there a GDAL-tool or some other commandline-suite that will do this?

It would be nice to have a tool that handles the shapefile natively (GDAL-tools do that) without having to import it or convert it to some internal format.

An alternate solution would be a commandline-interface to ftools.

THX stn


texts: PyQGIS - polygon to lines with geometry length and select features with condition

texts:
How can I make a script to convert polygon layer to explode lines ? After that I need to get geometry length for conditioned features. I have a part of code ?

import processing
processing.runalg("qgis:polygonstolines","Polygon.shp","Line.shp")
processing.runalg('qgis:explodelines',"Line.shp","Explode_line.shp")
layer = qgis.utils.iface.activeLayer()
vlayer.selectAll() - "QGIS crashes...."

for f in features:
    geom = f.geometry()
    print "Length:", geom.length()

Here need select on a map all features in Explode_line where ID_attribute > 10


texts: Mapfish 3 merge datasource with multiple maps

texts:
I would like to combine the examples of merging datasources and having multiple maps.

So instead of having a table datasource (1st example) I have a map datasource like the second example.

The goal is to output a report with a map in the first page, and other maps with features in the following pages (from datasource).

How to accomplish that? (config.yaml)


texts:
Is it possible to make QGIS 3 portable with changed location for Qgis setup parameters folder.

For QGIS 2 it looked like this in qgis.bat file:

  qgis --optionspath "X:\qgis" --configpath "X:\qgisini"

For QGIS 3 it doesn't work. Setup parameters are default placed in:

  C:\Users\xxx\AppData\Roaming\QGIS


texts: Erasing part of line vector outside raster region using GRASS?

texts:
I need to delete all the vector 'part' that it's outside the yellow raster region.

How can I do it?


texts: Creating measured line segments along existing polylines from starting point?

texts:
I realize my problem may be solved using Linear Referencing, and if that's the case I'm looking for help figuring out the workflow involved there.

I'm working in ArcGIS 10.2 Desktop Basic, and this is what I started with. I had a bunch of points where there were animal sightings, and those point were mostly in between two other points called "camera locations" that are start & end points for these animal sighting points. These are separate point features. I was supposed to connect the animal sighting points into a generalized line from the start to the ending camera locations, to have a line. I have this for about a dozen different camera start & end locations, a connect-the-dots line in between them now. Now, I'm supposed to start at the "starting" camera point side of the connect-the-dots line, and measure distances from 0-400', 400-800', 800-1200', and 1200-1600' but if the line ends before 1600' just cut it off at the ending camera location point on the line. To be clear, these are like 12 different polylines in a feature (shapefile), and they each need to be broken up into the measures above.

I'm wondering if someone can take my problem and help me make sense of how to get those measured points along the lines I drew from connecting the previous dots. While I've been using GIS for years, somehow I've never used linear referencing and I'm not sure which of its tools I need to use and when.


texts: How to add the same field to all tables in CASE dataset with magik?

texts: I need to loop through all tables of my smallworld repository and  add UpdateUser field to all of them. The looping is easy but how can get a handle to case_objects collection?


texts:
I am trying to write a Python Script for ArcGIS 10.2.

#import system modules
import arcpy
import math
from arcpy import env

#Set environment options
env.workspace="C:/Users/el/ac/script1"

# Set local variables
inFeatures = "line2.shp"
fieldName = "angle1"
expression = "GetAzimuthPolyline(!Shape!)"
codeblock = """import math
def GetAzimuthPolyline(shape):
    radian = math.atan((shape.lastpoint.x - shape.firstpoint.x)/(shape.lastpoint.y - shape.firstpoint.y))
    degrees = radian * 180 / math.pi
    return degrees"""

# Execute AddField
arcpy.AddField_management(inFeatures, fieldName, "SHORT")
#print "Field has been added."

# Execute CalculateField
arcpy.CalculateField_management(inFeatures, fieldName, "GetAzimuthPolyline(!Shape!)", "PYTHON", codeblock)

Field was added, but the result is not calculated and is not writing into the field. But, if I write the code block in Calculate field in ArcMap, field calculates correctly. Why?

Code was amended. But if I run the script:

Traceback (most recent call last):
  File "C:\Users\el\ac\script1\scr1.py", line 35, in <module>
arcpy.CalculateField_management(inFeatures, fieldName, "GetAzimuthPolyline(!Shape!)", "PYTHON", codeblock)
  File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\management.py", line 3354, in CalculateField
    raise e
ExecuteError: ERROR 000539: Error running expression: GetAzimuthPolyline(GPVARIANTOBJECT0)
Traceback (most recent call last):
  File "<expression>", line 1, in <module>
  File "<string>", line 3, in GetAzimuthPolyline
AttributeError: 'str' object has no attribute 'x'


texts: Making fill color and outline color match in symbology using ArcMap

texts:
How can you make the outline color the same as the fill color in symbology?
I have saved fill colors that i would like to make the outline colors as well.
I am using ArcMap 10.4.


texts: How to use the replace string function in Field Calculator?

texts:
Within an attribute table I'm trying to replace one string with another using the 'replace' function .  As an example, there are a number of features whose field 'NAME' contain the string '(B)' and I'm wanting to replace this with the word 'County'.  I'm selecting the 'update existing field' checkbox and using the following expression

 replace( '(B)','(B)','County')

The end result is that the field 'NAME' for every feature is replaced by 'County' irrespective of whether or not the field 'NAME' originally contained the string '(B)'.  Any help would be appreciated.


texts: X, Y coordinates in thousands to millions

texts:
I have a shapefile of contour lines and the X and Y coordinates are in thousands not in millions. How can I transform them in millions?

Like:

7802,305 to 7802305,00.


texts: How many layers should be in a MapService - what are the tradeoffs?

texts:
When setting up MapServices for ArcGIS, I'm often given an mxd that has many layers in it and told to publish it.

Should I have more mapservices with fewer layers in each, or fewer mapservices with more layers in each?

What are the tradeoffs ?


texts: Slow performing ArcGIS Near tool returns -1?

texts:
We are running an ArcGIS 10.0 geoprocessing script written in Python that (as part of many steps above & below) runs the Near tool.

Both the input and output datasets (shapefiles & personal geodatabase, respectively) are on a network drive, in the same projection, and we are not using a search distance (which by default = infinity, I believe).  When we run the near command, sometime it successfully runs to completion and sometimes (usually) it only partially works.

When it partially works (but runs to completion), the results typically have anywhere from 1/4 to 1/2 of the features with a NEAR_DIST value = -1 (and it is always the top portion of the records, as sorted by the OID).  In these cases, it also takes forever to process the top records, sometimes going from full success runs in 20 seconds to partial success runs in 1.5 hours.

Has anyone else experienced this kind of issue?  Any suggestions on WHY the issue is occurring or HOW to fix the issue?


texts:
I'm having difficulties reading the Gridded height data (OS Meridian 2). I'm loading with OSGB 1936 reference system but have trouble finding the lat/lon data for each point. i.e.I open the attribute table to see three columns FEATURE_SE, FEATURE_CO, and HEIGHT. I know the data loads correctly, as it produces the right shape map but don't know how to query the height at a position as I can't figure out how the Feature_se and feature_co link to lat/lon, or OSGB references.
Thanks if more information is needed please ask and I'll try and find it.


texts: Adding bar charts to maps in ArcGIS Pro?

texts:
When it comes to Symbology, is there a way of showing bar charts on the maps or at least in a layout in ArcGIS Pro?

This was possible in ArcMap and so I am just wondering if this is still a feasible task in ArcGIS Pro.


texts: If there is no spatial autocorrelation, should the spatial cluster detection be performed?

texts:
I have a regional count data of disease incidence. First, I calculated empirical Bayes estimate of standardized morbidity rate. Then I performed global Moran's I to check for spatial autocorrelation of standardized morbidity rate. The result showed that spatial autocorrelation is not significant.
Should I still perform detection of the location of a spatial cluster?


texts: Is it possible to store data on Amazon S3 for use in OpenGeo Suite?

texts: We've got about 500GB of geospatial data we'd like to "web-enable" to various degrees and OpenGeo Suite seems like a good route to go down given our circumstances. We can sort out a server (probably VPS) for a good price, but storing all that data locally on the same server will bump up the price considerably. So is it possible to have the VPS running the OpenGeo Suite applications, while the data is stored elsewhere (e.g., on S3)?


texts: Creating a giant matrix for use in CARBayes without running out of memory

texts:
I am trying to do an areal unit analysis using the package CARBayes. As part of the analysis, I am using the below code. my issue comes when I try to create the neighbour matrix with nb2mat. My sp object has 170,000 odd polygons in it so it can't make the matrix with the memory I have.
library(spdep)
library(CARBayes)

W.nb <- poly2nb(sp)

W <- nb2mat(W.nb, style = "B", zero.policy = TRUE)

test <- S.CARbym(case ~ covariate1),
             family = "poisson",
             data = sp,
             W = W,
             burnin = 10000,
             n.sample = 30000,
             thin = 20)

I found the below code in another thread to make a bigmemory matrix but CARBayes won't recognise it as a matrix.
My question is, does anyone know a way to use bigmemory or spam /sparse matrix or something similar to create the matrix so that it can be used in the CARBayes package without throwing an error saying the W isn't a matrix.

my_listw2mat = function (listw)
{
  require(bigmemory)
  n <- length(listw$neighbours)
  if (n < 1)
    stop("non-positive number of entities")
  cardnb <- card(listw$neighbours)
  if (any(is.na(unlist(listw$weights))))
    stop("NAs in general weights list")
  #res <- matrix(0, nrow = n, ncol = n)
  res <- big.matrix(n, n, type='double', init=NULL)
  options(bigmemory.allow.dimnames=TRUE)

  for (i in 1:n) if (cardnb[i] > 0)
    res[i, listw$neighbours[[i]]] <- listw$weights[[i]]
  if (!is.null(attr(listw, "region.id")))
    row.names(res) <- attr(listw, "region.id")
  res
}

my_nb2mat = function (neighbours, glist = NULL, style = "W", zero.policy = NULL)
{
  if (is.null(zero.policy))
    zero.policy <- get("zeroPolicy", envir = .spdepOptions)
  stopifnot(is.logical(zero.policy))
  if (!inherits(neighbours, "nb"))
    stop("Not a neighbours list")
  listw <- nb2listw(neighbours, glist = glist, style = style,
                    zero.policy = zero.policy)
  res <- my_listw2mat(listw)
  attr(res, "call") <- match.call()
  res
}

W <- my_nb2mat(W.nb, style = "B", zero.policy = TRUE)

test <- S.CARbym(case ~ covariate1),
             family = "poisson",
             data = sp,
             W = W,
             burnin = 10000,
             n.sample = 30000,
             thin = 20)


texts:
Despite similar Q&As I can't seem to get my answer.

So I have this set of coordinate (124616.1535508, 486555.4631204) are mapped to the Bessel_1841 coordinate space. I'd like to transform these coordinates into a latlong's that can plotted in google maps.

I've tried several options none of which give me the desired location:

import utm
print(utm.to_latlon(124616.1535508, 486555.4631204, 32, 'U'))
(4.394248331543074, 5.61858563235919)

import pyproj
p = pyproj.Proj(init='epsg:27700')
lon, lat = p(124616.1535508, 486555.4631204, inverse=True)
print(lat, lon)
(54.2004456912487 -6.223076560046059)


texts: Creating separate lines for all pairs of points from a network

texts:
I have a network layer of roads and a point layer containing points (towns) on that network.
I need to create layers containing each pair of towns and its true shape route line. So for example from A-B, A-C, A-D, ... etc. So that every pair has its own shapefile with the line.

I tried the OD Matrix solution in Network Analyst in ArcMap 10.3 which created what I wanted as far as attributes go, but it draws only straight lines between each pair, it doesn't draw the line along the network. I need that.
Closest Facility solution provides true shapes along the network but it doesn't create what OD matrix does. I need one shapefile for every pair.

Any advice or direction?

I tried to implement what Paul in Answer 1 suggested but I couldn't complete the script. I still need a tool to address correct points in pairs so they could be solved. It could be done sql query maybe? Any further advice?

from collections import defaultdict
import itertools
import pprint

import arcpy
from arcpy import env
import os

# generate pairs
ddict = defaultdict(list)
    for start,end in itertools.combinations(range(10), 2):
    ddict[start].append(end)

# viz example 3
try:
#Check out the Network Analyst extension license
arcpy.CheckOutExtension("Network")

#Set environment settings
env.workspace = r"E:/Dokumenty/GISdata/PID"
env.overwriteOutput = True

#Set local variables
inNetworkDataset = r"E:/Dokumenty/GISdata/PID/silnice_krnap_ND.nd"
impedanceAttribute = "Length"
accumulateAttributeName = ["Length"]

for key,value in ddict.items():
    for v in value:
        inFacilities = r"PID/Obce_body_luzka" # melo by byt key
        inIncidents = r"PID/Obce_body_luzka" # melo by byt v
        outNALayerName = "Dvojice_bodu_{}_{}".format(key, v)
        outLayerFile = os.path.join(r"E:/dokumenty/gisdata/output2", outNALayerName + ".lyr")
        NAResultObject = arcpy.na.MakeClosestFacilityLayer(inNetworkDataset,outNALayerName,impedanceAttribute,"TRAVEL_TO","",1, accumulateAttributeName,"NO_UTURNS")

        #Get the layer object from the result object. The closest facility layer can
        #now be referenced using the layer object.
        outNALayer = NAResultObject.getOutput(0)

        #Get the names of all the sublayers within the closest facility layer.
        subLayerNames = arcpy.na.GetNAClassNames(outNALayer)
        #Stores the layer names that we will use later
        facilitiesLayerName = subLayerNames["Facilities"]
        incidentsLayerName = subLayerNames["Incidents"]

        # Import system modules
        import arcpy

        # Set data path
        intable = "E:\Dokumenty\GISdata\PID\Obce_body_luzka.dbf"

        # Get the fields from the input
        fields= arcpy.ListFields(intable)

        # Create a fieldinfo object
        fieldinfo = arcpy.FieldInfo()

        # Iterate through the fields and set them to fieldinfo
        for field in fields:
            fieldinfo.addField(field.name, field.name, "VISIBLE", "")
            arcpy.MakeTableView_management(intable, "crime_view", "", "", fieldinfo)

        arcpy.na.AddLocations(outNALayer, facilitiesLayerName, "Obce_body_luzka", "", "")
        arcpy.na.AddLocations(outNALayer, incidentsLayerName, "Obce_body_luzka", "", "")

        #Solve the closest facility layer
        arcpy.na.Solve(outNALayer)

        #Save the solved closest facility layer as a layer file on disk with
        #relative paths
        arcpy.management.SaveToLayerFile(outNALayer,outLayerFile,"RELATIVE")

print "Script completed successfully"

except Exception as e:
# If an error occurred, print line number and error message
import traceback, sys
tb = sys.exc_info()[2]
print "An error occured on line %i" % tb.tb_lineno
print str(e)

I worked with the idea of using the closest facility method, while thinking that "the closest" route is also "the only route" from A-B, etc.


texts: Calculating the average of a raster in raster calculator

texts:
I have a raster containing trout sightings spanning back to the 1800's and I want to calculate the average of all the data and then compare this to the  year with the most sightings.

Can anyone explain how to go about this please?

I can't seem to find any help about this elsewhere.


texts: pgrouting - network distance many-to-one shortest paths

texts:
I want to find the distance (km) between each origin point to the closest destination point (restaurants) using pgrouting.

I have 3 tables

my road network with nodes and edges
140 origin points - centroids
of neighbourhoods
3000 destinations - restaurants

The road network already has source, target, length (cost) and topology.

I have used tutorials that specify in the SQL query which nodes to calculate cost between such as this one:

SELECT seq, id1 AS source, id2 AS target, cost FROM pgr_kdijkstraCost(
'SELECT id, source, target, cost FROM edge_table WHERE cost >= 0',
10, array[4,12], false, false

Instead of explicitly referencing origin as "Node 10" and Destinations as "Nodes[4,12]" is there a way to reference my tables that list all 3,000 destinations?


texts: getFeatures() is inaccurate with memory layers

texts:
I was writing a plugin for QGIS when I wanted to copy all my data to memory layers and edit them instead of actual shapefiles to speed up the process. In a QgsMapTool-derived class, I am using a method call like this

def canvasPressEvent(self, e):
    mapCoords = self.toMapCoordinates(e.pos())
    features = self.__layer.getFeatures(QgsFeatureRequest(QgsRectangle(mapCoords, mapCoords)))
    #...

This worked just fine for the old, file-based layers, however, for the new in-memory ones, I'm getting horribly inaccurate results.

My code for copying the layer is pretty straightfoward:

with fiona.open(path) as f:
    geometry = f.schema["geometry"] # yeah, this is a bit hacky, sry

originalLayer = QgsVectorLayer(path, name + "_temp_layer", "ogr")

# build the creation string
props = geometry + "?"
props += "crs=WKT:" + originalLayer.crs().toWkt()
props += "&"
props += "index=yes"

# create layer and start filling it
layer = QgsVectorLayer(props, name, "memory")
#layer.setCrs(originalLayer.crs()) #this doesn't seem to change anything
layer.startEditing()

# copy fields
for field in originalLayer.fields():
    layer.addAttribute(QgsField(field.name(), field.type()))

# copy features
newFeatures = []
for feature in originalLayer.getFeatures():
    newFeature = QgsFeature()
    newFeature.setGeometry(feature.geometry())
    newFeature.setAttributes(feature.attributes())
    newFeatures.append(newFeature)
result = layer.dataProvider().addFeatures(newFeatures)

# commit
layer.commitChanges()
layer.updateExtents()

Now, this code seems to work just fine, the result looks just as it looked before, but when using the tool, I can click on a spot a few centimetres away from a feature and the tool detects it as if I clicked the feature itself. This also creates massive overlap between the features. These things worked just fine with file-based layers.

It feels like something is wrong with the CRS... Do I somehow have to apply the CRS to be taken into account during spatial querys?

I already tried manually creating spatial indices, but that yielded the same results.


texts: Point data extraction and printing as a TABLE

texts:
I am working with MODIS-NDVI image collection in Google Earth Engine. Below, you can find a script which is creating a chart of NDVI values over the given period for particular location. I would like to print it as a table or text in Console/Panel, but could not find code for it. Is there anyway to do it?

var point = ee.Geometry.Point([10.455856189782935, 50.440969544991205]);

var image = ee.ImageCollection('MODIS/006/MOD13A1').
      filterDate('2018-01-01','2018-12-31')
        .select('NDVI');

var div_ndvi = function(image) {
  var NDVI_0 = image.select("NDVI");
  var NDVI_1 = NDVI_0.divide(10000);
  return image.addBands(NDVI_1);
};

var withN = image.map(div_ndvi);

print(withN);

var chart = ui.Chart.image.seriesByRegion(withN, point, ee.Reducer.first(), 'NDVI_1', 500);

print(chart)


texts: Error with feature to raster conversion in ArcGIS for Desktop?

texts:
I am trying to convert an ESRI geodatabase into a raster using the 'Feature to raster' conversion tool in ArcMap 10.1. The geodatbase is available publicly from the USDA Forest Service website (South Sierra, zone 4), and I am using the WHRTYPE property to convert it in raster form. The output cell size is set to 30 (see attached picture).

The problem is that ArcGIS simply crashes during the operation, returning a serious application error popup. No explanations of what went wrong. A tif file is actually produced, although the corresponding 'WHRTYPE' field in its attribute table is empty.

I tested this on two different, quite recent machines and got the same problem. I checked the field for null values. I also disabled background geoprocessing, as suggested in other posts.

I tested the same process on another geodatabase from the same webpage (South Coast, zone 7) and it worked. So I am really confused at the moment. Can anyone help me out?


texts: Error: Attempt to create 2x0 dataset is illegal,sizes must be larger than zero

texts:
I have the following input file

gdalinfo SAVI_v2_T_0006_0006_savi.tif
Driver: GTiff/GeoTIFF
Files: SAVI_v2_T_0006_0006_savi.tif
Size is 4308, 208
Coordinate System is:
PROJCS["WGS 84 / UTM zone 32N",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",9],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","32632"]]
Origin = (742430.000000000000000,5870860.000000000000000)
Pixel Size = (15.000000000000000,-15.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  742430.000, 5870860.000) ( 12d36'26.78"E, 52d55'57.16"N)
Lower Left  (  742430.000, 5867740.000) ( 12d36'18.40"E, 52d54'16.39"N)
Upper Right (  807050.000, 5870860.000) ( 13d33'57.69"E, 52d53'58.27"N)
Lower Right (  807050.000, 5867740.000) ( 13d33'47.10"E, 52d52'17.62"N)
Center      (  774740.000, 5869300.000) ( 13d 5' 8.33"E, 52d54'10.83"N)
Band 1 Block=4308x1 Type=Float32, ColorInterp=Gray
  NoData Value=1.175494351e-38

And I would like to decrease its size by 1024x. I specify the target resolution for gdalwarp to be 30720 = 15 (pixel size in meters) * 1024 as follows:

gdalwarp -tr 30720 30720 SAVI_v2_T_0006_0006_savi.tif out.tif

and get

ERROR 1: Attempt to create 2x0 dataset is illegal,sizes must be larger than zero.

I suppose that formally the above command is correct but

why those 208 pixels are squashed into 0 pixels instead of 1 pixel?
how to deal with this situation?
is there any other reliable way to reduce image size by an arbitrary scale factor?


texts:
Recently got enamored with the idea of doing a time-based and an energy-based least cost path analysis. After doing some research, I see that GRASS is the best for executing time-based given the k-move functionality. I'm only accustomed to ArcGIS, so I was wondering if there are any GRASS tutorials around that might be helpful, specifically to r.walk, etc.

Some more specifics on what I'm doing: using ASTER gDEM v2 (30m res) in Turkey study area for archaeological project. Doing several points of origin/destination across a mntn range.


texts: Converting latitude, longitude to image coordinates

texts:
I have a 2D image and I know a given pixel in the image corresponds to a particular latitude longitude (so basically a reference point). The image I have is a camera image of an outdoor scene and I have the possibility to identify the (lat, long) of some of the points on that camera scene.
How could one transform a given latitude, longitude to a pixel coordinates given this information. I am guessing just one reference point is not enough and perhaps we need at last two to establish the spatial resolution. Is that sufficient information to transform latitude, longitude to image coordinates.


texts:
I'm not very familiar working with WMS/WFS. I found this list of datasets from Cyprus but I'm not able to open them in any software. Here is a link to the list of datasets:

http://freegisdata.org/place/50821/

For example I want to open Soil Threats;

http://freegisdata.org/record/urn:uuid:eccbd308-1580-45f4-ad18-87bf3cb1e44c-MS/

I tried to add WMS layer to QGIS but I'm not sure if I'm entering the right url, It doesn't work.


texts:
In OpenLayers 5, if loading of some tiles has failed (because of temporary network issue, for example), those tile are never attempted to load again.

Is it possible to force OL to retry loading of those tiles?

Workflow:

disable network
zoom in (to zoomlevel 2), all tiles fail
zoom in again
enable network
zoom out, you are back in zoomlevel 2 now. No tiles are displayed <- OL should reload those errorneous tiles

Goal is not to reload whole layer - OpenLayers should attempt to reload only failed tiles every time, when those are displayed on map.


texts: Issue with Calculate Field and text array

texts:
I've tried this expression in the Arcade playground and it works as I expect.
Number(Split($feature.PVI,"+")[1])
$feature.PVI is a text field that contains something that looks like "R+15".
What should be entered into the field is 15. When I try to use the same expression in the Calculate Field tool I get ERROR 002717: Invalid Arcade expression, Arcade Error: Index out of bounds.


texts: How to group 10k points into closest pairs?

texts:
I was first faced with a problem of grouping 1000 pairs into unique pairs to minimise the total distance - since this was a small number I solved it using linear optimisation:

Part 1 - Basic/Naive approach:

A bit more detail on it here but the gist is:

Create an objective function which wants to minimise all the entries in the lower triangle of a distance matrix (multiplied by the respective dummy). For example v1_10 means that point 1 and point 10 are connected. Such that each store is connected to only one (other store).

The result of the optimisation takes us from the first picture to the second. The structure of the code is thus:

shops = [[id, latitude, longitude],...]

The helper functions to create the optimisation:

def store_connect(shops):
    """
    Cross stores and return connection e.g. "v1-2" by ID
    """
    for i in range(len(shops)-1):
        for j in range(i+1,len(shops)):
            yield 'v'+str(shops[i][0])+'-'+str(shops[j][0])

def distances(shops):
    """
    Return distance in miles for crosses
    """
    for i in range(len(shops)-1):
        for j in range(i+1,len(shops)):
            yield haversine([shops[i][1],shops[i][2]],
                            [shops[j][1],shops[j][2]])

def all_constraints(shops):
    """
    Return constraint row
    """
    for a in shops:
        cons = []
        for b in shops:
            if a[0]<b[0]:
                cons.append("v"+str(a[0])+"-"+str(b[0]))
            elif a[0]>b[0]:
                cons.append("v"+str(b[0])+"-"+str(a[0]))
        yield cons

And then the optimisation itself:

prob = LpProblem("Minimise distance",LpMinimize)

# Variables
x = []
xnames = []
for one in store_connect(shops):
    xnames.append(one)
    x.append(LpVariable(one,0,1,LpInteger))
print("Created %d variables" % len(x))

# Objective
prob += pulp.lpSum([i[1]*x[i[0]] for i in enumerate(distances(shops))])
print("Created objective")

# Constraints
counter = 0
for constraint_rw in all_constraints(shops):
    prob += pulp.lpSum(x[xnames.index(one_cn)] for one_cn in constraint_rw) == 1
    counter += 1
    if counter % 10 == 0:
        print("Added %d contraints out of %d" % (counter,len(shops)))

Which is straight-forward to solve directly or upload the .lp file to CPLEX:

#Using PULP - if not too complicated ... otherwise CPLEX:
prob.solve()
# The status of the solution is printed to the screen
print("Status:", LpStatus[prob.status])

# Optimum variable values
sol = []
for v in prob.variables():
    if v.varValue:
        sol.append(v.name)
print("Total minimised distance: ", value(prob.objective))

shp_names = [row[0] for row in shops]
edges = []
for x in sol:
    v_list = (x)[1:].split("_") #remove the first letter 'v' and split
    conn = [it for sub in [shops[shp_names.index(float(v_list[0]))][1:3],
                           shops[shp_names.index(float(v_list[1]))][1:3]] for it in sub]
    edges.append(conn)

# Print diagram
plot_points(shops,edges)

Part 2 - Better approach (for lots of points):

If I expand the problem to 10,000 points I get roughly 50,000,000 variables. I was thinking of the following ways of simplifying the problem:

1) Use clustering to group the 10,000 points into something like 100 groups and then solve 100 times within those groups (however I'm not sure how to check whether the solution would be optimal).

2) Project the points; put them into a quad-tree and then instead of allowing connections between one point and any other point -> only allow connections to the closest 100 points (e.g.): however, I am not sure how to properly set this (the cut-off). I thought perhaps if a point is chosen that is the cut-off point then the algorithm will restart with a larger cut-off.

3) Run something like a minimum spanning tree or concorde tsp algorithm and then back-track the solution from that?


texts:
How can I do matrix processing in ArcGIS or ERDAS Imagine?
l have two different raster images for change detection. Matrix join same code filed information.

In this video link already done matrix in ERDAS and then process continue in ArcGIS.
https://www.youtube.com/watch?v=qg3xE_YAgWA&t=519s


texts: Using the GeoWebCache expireCacheList parameter

texts:
I've just started looking at the cache expiration parameters for GeoWebCache and came across something that really confused me. In the documentation it says in its explanation for the <minZoom> parameter;

  First zoom level. The list should always start with 0 and be
  monotonically increasing.

For a start, I don't understand what monotonically increasing really means. With a bit of saerching it looks like it means that it should go up in regular steps. It then gives this as an example

<expireCacheList>
  <expirationRule minZoom="0"  expiration="14400" />
  <expirationRule minZoom="10" expiration="7200" />
</expireCacheList>

So lets say I wanted to cache everything between 0 and 9 for a month, 10 and 15 for a week, and 16 and above for a day. Could I do it like this.

    <expireCacheList>
      <expirationRule minZoom="0"  expiration="[insert num seconds in month here]" />
      <expirationRule minZoom="10" expiration="[insert num seconds in week here]" />
      <expirationRule minZoom="15" expiration="[insert num seconds in day here]" />
    </expireCacheList>

Or would it have to be like this

    <expireCacheList>
      <expirationRule minZoom="0"  expiration="[insert num seconds in month here]" />
      <expirationRule minZoom="5"  expiration="[insert num seconds in month here]" />
      <expirationRule minZoom="10" expiration="[insert num seconds in week here]" />
      <expirationRule minZoom="15" expiration="[insert num seconds in day here]" />
    </expireCacheList>

Does anyone also happen to know off hand if this requires a reload of the GeoServer configuration or whether it applies as soon as you save the XML file. The documentation is a bit flaky and I'd like to understand it a bit better before going any further. I don't want to be responsible for bringing down our GeoServer (again...) so any help or pointers is much appreciated.


texts: How can I remove default zoom in/out [+/-] buttons from openlayers to keep the navigation bar on?

texts:
I want to remove the default zoom +/- buttons from openlayers.

Here is my code

    <script type="text/javascript">
    var map, ia_wms;
    function init(){
map = new OpenLayers.Map('map', {zoomDuration: 1,projection: 'EPSG:3857'
***/*controls: [
    new OpenLayers.Control.Navigation(),
    new OpenLayers.Control.ArgParser(),
    new OpenLayers.Control.Attribution()
    ]*/***
    });
    var sr_all = new OpenLayers.Layer.WMS("EOL","http://192.168.0.233/cgi-bin/mapserv?map=/var/www/mstest/plant.map",{'layers':"map_outline",'transparent':true,'format':'image/png'}
           ,{
            isBaseLayer: false,opacity: 100
            },
            {singleTile: true, ratio: 1}
            );
    var gmap = new OpenLayers.Layer.Google(
    "Google Streets", // the default
    { numZoomLevels: 25,visibility: false}
    ,{
        isBaseLayer: true
    }
    );
    map.addLayers([sr_all,gmap]);
    map.addControl(new OpenLayers.Control.PanZoomBar());

    map.zoomToExtent(new OpenLayers.Bounds(69.728426,22.311238,69.774227,22.356185));
    map.setCenter(new OpenLayers.LonLat(69.754227,22.332185).transform(
    new OpenLayers.Projection('EPSG:4326', 'EPSG:3857'),
    map.getProjectionObject()
    ), 14);
    }
</script>

when I added,
        controls: [
        new OpenLayers.Control.Navigation(),
        new OpenLayers.Control.ArgParser(),
        new OpenLayers.Control.Attribution()
        ]
to the map = new OpenLayers.Map, to remove default zoom +/- buttons, the map is not loading.
Please see the image also.

Both the navigation buttons exist. :(

I want to remove the small- 2 button- one.
I changed the code like this:

        map = new OpenLayers.Map('map', {zoomDuration: 1,projection: 'EPSG:3857'
    controls: []
    });

and added

        map.addControl(new OpenLayers.Control.PanZoomBar());
    map.addControl(new OpenLayers.Control.KeyboardDefaults());

But, the map is not loading. Do not know any syntax error..


texts: Is there a means of recording all user actions in QGIS?

texts:
Is there a means of recording all user actions in QGIS, in particular date/time of data access?

Most of our data is in Postgresql which is fine, I can run logging from there and QGIS is recorded there as the client application. What I'm stuck on is how to track access to flat files such as rasters.


texts: Reduce multiple images per day to a single image in time series collection (Earth Engine)

texts:
I have a code which coverts the MODIS Aqua dataset into monthly means. The code is here.

However, I'm want to do the same for a new data set, the NOAA/CDR/SST_WHOI/V2. However, rather than single daily values, this dataset has values 8 times a day. Is there a way of reducing this so it just collects one image per date, which I can then use the below code on?

var startyear = 2010;
var endyear = 2011;

var startmonth = 1;
var endmonth = 12;

var startdate = ee.Date.fromYMD(startyear, startmonth, 1);
var enddate = ee.Date.fromYMD(endyear , endmonth, 1);

var years = ee.List.sequence(startyear, endyear);

var months = ee.List.sequence(1, 12);
var monthlySST =  ee.ImageCollection.fromImages(
years.map(function (y) {
return months.map(function(m) {
  var w = UK_SST.filter(ee.Filter.calendarRange(y, y, 'year'))
                .filter(ee.Filter.calendarRange(m, m, 'month'))
                .mean();
  return w.set('year', y)
          .set('month', m)
          .set('system:time_start', ee.Date.fromYMD(y, m, 1));

});
}).flatten()
);


texts: Choosing dpi to export images from ArcMap?

texts:
I need to export images from ArcMap to embed them in a pdf via LaTeX. Some images should fit a whole page, while others should cover an area equal to approximately 1/4 of the paper page.

If I want to export them with an .eps format, what is the dpi I should give them?


texts: Reducing size of location cursor on ArcGIS Collector app?

texts:
In the Collector app the blue location cursor covers my point features when I get close to one I am navigating to. The point feature then disappears behind the blue dot and I can't tell if I'm centered on the blue dot or near the edge. I need to be precise (w/in 1m), as this location is used for repeat measurements and I would like to avoid installing permanent plot markers if possible.

As a work around I have simply increased the size of my point features to 15pt so they are slightly larger than the blue cursor. This results in a bunch of large dots on the map, which just looks messy.

I would like to reduce the size of the blue location cursor in Collector so I could then make my pnt features small again, thus resulting in a less cluttered map.

Does anyone know if this is possible and if so, how? I have tried a brief internet search but haven't found anything that addresses what I want to do.