texts: How do I convert an unsigned raster to a signed raster?

texts:
I'm using ArcGIS 10.2.  I need to convert two rasters from 8bit unsigned to an 8bit signed raster so I can do a subtraction and get both positive and negative values.

I've tried to use the copy raster tool that allows you to set the raster pixel type, but it always comes out the other side as an unsigned raster.

I tried converting my GRID to a .tiff, then using the copy command, but the Copy Raster tool will only output a GRID.  After converting and copying I still got an unsigned raster output.


texts: QGIS how to compute a trendline of two point lines?

texts:
Does anyone know how to compute a trendline for two lines of points? I have two lines of points of a 10 km long urban road. 1st line is one driving direction and the 2nd line of points ist the other driving direction.

Is there a way to get a trendline of these both point lines also as a point line?

Thanks!!


texts: Updating Fusion Table layers from csv files using python

texts:
I'm trying to use Python to take csv files and use them to update (or replace) tables in Fusion Tables.  I found the "Python client library for Fusion Tables" on googles' site and downloaded them as well as the two required dependencies:

python-oauth2 library: http://github.com/simplegeo/python-oauth2
python-oauth2 dependency, httplib2: http://code.google.com/p/httplib2/

I've installed all three on python 2.7 like required and just can't seem to figure out how to actually get it to work. I've even tried to get their samples to work with no luck either.  At this point let me say that I'm not a programmer and don't know python.  However, I can grasp it enough (and look at and modify code) pretty well though and used python to write the scripts that export from SDE to shapefiles, then convert the shapefiles to KML using GDAL and then convert the KML to CSV using elementtree.

Does anybody know of a good resource to look at for help in doing what I'm trying to do?  All I can find online is a few people saying they did it, with no explanation, and a ton of links to the google code (which has no documentation)

Thanks in advance,
Adam


texts: PostGIS - syntax error inserting GeoJSON

texts:
I am trying to insert a GeoJSON into PostGIS, but I am encountering a syntax error.  I can't seem to see where the problem is - I'm sure it is an easy issue to resolve, but my inexperienced eyes can't seem to see where the problem is.

Insert statement (PHP):

$conn = pg_connect("host=localhost dbname=postgis user=$user password=$password");

$sql = "INSERT INTO carto_doodle (poly_geom) VALUES
    (ST_SetSRID
        (ST_GeomFromGeoJSON
            ({"type":"Feature",
              "properties":{},
              "geometry":{"type":"Polygon",
                          "coordinates":
                              [[[-74.54635620117188,40.773261878622634],
                                [-74.69741821289062,40.61082491956405],
                                [-74.30740356445312,40.61603737424185],
                                [-74.54635620117188,40.773261878622634]]]
                          }
              })
          ),
     4326)"

$success = pg_query($conn, $sql);

However, this returns an error for the $success line:

Warning: pg_query(): Query failed: ERROR: syntax error at or near &quot;{&quot; LINE 1: ...(poly_geom) VALUES (ST_SetSRID(ST_GeomFromGeoJSON({&quot;type&quot;:&quot;F... ^

Here's what ended up working for me...

$sql = 'INSERT INTO carto_doodle (point_geom) VALUES
       (ST_SetSRID(ST_GeomFromGeoJSON(\'
            {"type":"Point",
             "coordinates":
                          [-48.23456, 20.12345]
            }\'
        ), 4326))';


texts: Can anyone help me find the link to download ArcGIS Engine 10? I have tried to find it, but unfortunately I have not found it.


texts:
I have a 3D network of streams represented using linestrings with Z-coordinates. I would like to rasterise this feature, where the raster values represent the linearly interpolated Z-values between 3D coordinates from the linestring verticies. I have all open source and Esri tools available, but I don't see any simple solution.

Did I miss an available tool to do this conversion, or do I need to make one myself? (I have a few ideas floating around my head, but it would probably take the better part of an afternoon to code out.)


texts: "ERROR: Unable to load GDAL library" occurs when trying to GRASS functions in QGIS

texts: I'm new to QGIS. I'm using GRASS operations to analyze raster layers and getting ERROR: Unable to load GDAL library in QGIS. Maybe the path is wrong which should point the GRASS operations to GDAL. Any tips?


texts: QGIS Server does not publish PostGIS, Oracle or SQLITE Layer

texts:
I want to publish a Postgis Layer which is defined in a QGIS Project with the QGIS Server. The capabilities document of the QGIS Server does not include the layer description of the  Layer. If i try this for example with a local shp Layer it will work. Is there a solution about this.

My Version of QGIS is 2.6 for Windows 64 bit and i installed the QGIS Server over the OSGEO setup on the same system.

The REQUEST=GetMap return a blank image


texts: Getting tile number of sinusoidal MODIS product from lat/long?

texts:
I'm using this tool to get the tile number (horizontal/vertical) from lat/long. Is there any code,preferably python, to get the tile number ? I tried using this tool. But it gives wrong tile number. Here's my code:

import numpy as np
data = np.genfromtxt('sn_bound_10deg.txt', skip_header = 7, skip_footer = 3)
lat,lon = -25.439538,149.053418
in_tile = False
i = 0
while(not in_tile):
  #print data[i,0],data[i,1],data[i, 4],data[i, 5],data[i, 2],data[i, 3]
  in_tile = lat >= data[i, 4] and lat <= data[i, 5] and lon >= data[i, 2] and lon <= data[i, 3]
  i += 1

vert = data[i-1, 0]
horiz = data[i-1, 1]
print('Vertical Tile:', vert, 'Horizontal Tile:', horiz)

The output tile from the above code is: h30v11

Expected output and output from the online calculator is: h31v11

Here's the sn_bound_10deg.txt file.


texts:
I'm trying to find a shortest path between 2 nodes using djikstra method but it keeps saying

  ERROR:  Error, query must return columns 'id', 'source', 'target' and 'cost'

this is my query:

SELECT seq, id1 AS node, id2 AS edge, cost
FROM pgr_dijkstra('SELECT gid, source,target, ST_Length(geom2d) as cost FROM  public.canada_rail_net', 1, 200, false, false);

is there anything wrong with column names?


texts: Feature/label visibility based on spatial relationship to index feature using data-driven pages

texts:
I'm creating a series of maps using ArcGIS data driven pages. I know how to display features that share an attribute with the index layer, but is it possible to display features (and similarly, labels), based on their proximity (i.e. within a certain distance) of the current index layer feature? I know I can do this using the Atlas Generator in QGIS but I'd like to achieve the same result in ArcGIS.

Here I have an example of my index layer in ArcMap Desktop (polygon with dashed outline). Suppose the point data shown in the map does not share any attributes with the index layer - can I display only the points within the polygon (or close to it, for example)?

I can do this in QGIS 2.12 (above) using rule-based styling and tests between the geometry of the point feature and the geometry of the current atlas feature. And I can display labels based on these tests. Is there a similar alternative in ArcGIS?

I'd rather not have to create a new layer using a spatial join to add the index attribute to the other layer(s), and a spatial join seems to be limited to only features within the polygon.


texts: Query Overpass Turbo for postal code boundaries

texts:
I'm looking to query Overpass Turbo for postal code boundaries, specifically all within the USA. Here's what I found so far and where I'm stymied:

This query for postal code boundaries in Liepzig works:

area[name="Leipzig"][admin_level=6][boundary=administrative]->.myarea;
rel(area.myarea)["boundary"="postal_code"];
out geom;

This query for boundaries of the State of Oregon works:

area[name="Oregon"][admin_level=4][boundary=administrative]->.myarea;
rel(pivot.myarea);
out geom;

But this query for postal code boundaries inside the State of Oregon does not:

area[name="Oregon"][admin_level=4][boundary=administrative]->.myarea;
rel(area.myarea)["boundary"="postal_code"];
out geom;

Thoughts on why this last query doesn't work while the first two do?


texts:
Running psql (13.1, server 12.4) trying to install PostGIS.
OS X 10.13.6.
postgis installed via brew install postgis.
CREATE EXTENSION postgis;
ERROR:  could not load library "/usr/local/lib/postgresql/plpgsql.so": dlopen(/usr/local/lib/postgresql/plpgsql.so, 10): Symbol not found: _CachedPlanAllowsSimpleValidityCheck
  Referenced from: /usr/local/lib/postgresql/plpgsql.so
  Expected in: /usr/local/opt/postgresql/bin/postgres
 in /usr/local/lib/postgresql/plpgsql.so

While the source code referenced can be found for this symbol, I have no idea what is going on.
How to rectify it?


texts:
I use the following code in order to change the fillColor of a marker within my leaflet-client, as suggested on this post: https://gis.stackexchange.com/a/124928/20297

var marker = marker = L.circleMarker(latlng, myMarkerOptions).addTo(map);
marker.valueOf()._icon.style.fillColor = 'green'

However when I execute this, I get the following error:

  Cannot read property 'style' of undefined

I also tried to directly change the markers color:

marker.fillColor = 'green'

which simply does absolutely nothing.

My browser is an Opera 66. However I also checked FireFox to no avail.


texts:
I asked a question here:
https://stackoverflow.com/questions/14790463/qgis-copy-layer-in-c-plugin
and I was suggested to find the library which provides QgsOgrProvider::createEmptyLayer() function.

Can anyone help me with this?


texts: Re-scale and move vector layer manually in QGIS

texts: I have imported a dxf file (site layout drawings for proposed development), but its scaling is too small and it does not line up with the google street maps. I want to re-scale the layer manually based a common feature on the aerial photo and correlating the length of the same feature on the dxf import. I then want to move the layer manually to line up with a particular point on the map. Any advice would be appreciated.


texts: How to crop Geotiff image in C++ (GDAL or other possibilities)?

texts:
I am developing a code to process a geotiff image. But firsly I need to crop a geotiff image preserving its metadata (projection...) in C++. My idea it's to provide the corners and the geotiff image to a function inside the code and then it return the subset. I thought about using GDAL but I don´t know if there is any code or function.

Anybody know about this?


texts: ListFields (arcpy) not reading from layer

texts:
I have a very simple python script that opens an mxd from a folder and then loops through each layer, and finds a ObjectID field and changes the alias name.

#setup
import arcpy
from arcpy import env

#Go into folder where mxd's are
env.workspace = r"D:\Temp\Test_Mxd"
#for each mxd
for mxd in arcpy.ListFiles("*.mxd"):
    print "Working on map document: " + mxd
    mapdoc = arcpy.mapping.MapDocument("D:\\Temp\\Test_Mxd\\" + mxd)
    #   for each dataframe
    df = arcpy.mapping.ListDataFrames(mapdoc,"Layers")[0]
    ##layers = arcpy.mapping.ListLayers(mapdoc,"*",df):
    layers = arcpy.mapping.ListLayers(mapdoc,"*",df)

    #   for each layer
    for layer in layers:
        print "Working on layer: " + layer.name
        **fieldList = arcpy.ListFields(layer)**
        for field in fieldList:
            if field.name == "OBJECTID":
        # Set the new alias name
                field.aliasName == "IGNORE::OBJECTID"
    mapdoc.save()
del mxd

my problem is on the line fieldList = arcpy.ListFields(layer), I get a error in the interpreter (pythonwin) which says that the layer does not exist.
Has anyone got any idea on where I am going wrong?


texts: Geoserver ImagePyramid/ImageMosaic with Time Dimension

texts:
I am having trouble to set the time range of image mosaics granules. The problem is, that it would be necessary to set a valid_from and valid_till for each granule - as far as I could see, there is no solution to this problem. It is just possible to display the picture 1 day - i tried every setting and changed the attributes, but it remains 1 day.

How to set a time range for a image mosaic?
If it is not possible: How to set a time range for the granules of a image mosaic?

The documentation mentions a reference value and a range in this format:
2001-01-01/2002-03-23

But I couldnt find a way to set it for raster data. Has anybody an idea?


texts: Converting PostGIS table to Feature Set using ArcPy

texts:
I am writing an ArcPy script (ArcGIS Pro) to export tile in the area of interest with help of arcpy.ExportTileCache_management(). I am getting the area of interest from a table in PostGIS. As far as I read in the documentation are the area of interest must be Feature Set data type.
I tried FeatureSet.from_json to convert my table in PostGIS but it was not successful. I kept getting the following error
File "C:\Program Files\ArcGIS\Pro\bin\Python\envs\arcgispro-py3\lib\site-packages\IPython\core\interactiveshell.py", line 3326, in run_code
    exec(code_obj, self.user_global_ns, self.user_ns)
  File "<ipython-input-2-90516b90cd60>", line 1, in <module>
    runfile('C:/Users/ncipa/PycharmProjects/arcpy/test.py', wdir='C:/Users/ncipa/PycharmProjects/arcpy')
  File "C:\Program Files\JetBrains\PyCharm Community Edition 2018.3.5\helpers\pydev\_pydev_bundle\pydev_umd.py", line 197, in runfile
    pydev_imports.execfile(filename, global_vars, local_vars)  # execute the script
  File "C:\Program Files\JetBrains\PyCharm Community Edition 2018.3.5\helpers\pydev\_pydev_imps\_pydev_execfile.py", line 18, in execfile
    exec(compile(contents+"\n", file, 'exec'), glob, loc)
  File "C:/Users/ncipa/PycharmProjects/arcpy/test.py", line 518, in <module>
    areaofinterest = FeatureSet.from_json(polygon)
  File "C:\Program Files\ArcGIS\Pro\bin\Python\envs\arcgispro-py3\lib\site-packages\arcgis\features\feature.py", line 694, in from_json
    return FeatureSet.from_dict(json.loads(json_str))
  File "C:\Program Files\ArcGIS\Pro\bin\Python\envs\arcgispro-py3\lib\json\__init__.py", line 348, in loads
    'not {!r}'.format(s.__class__.__name__))
TypeError: the JSON object must be str, bytes or bytearray, not 'dict'

How can I convert the record in my table in PostGIS to Feature Set?
# Connect and query tables in Postgresql
try:
    connection = psycopg2.connect(user='xxx',
                                  password='xxx',
                                  host="xxx",
                                  port="5432",
                                  database="xxx")
    cursor = connection.cursor()

    # Fetch records
    postgreSQL_select_Query1 = "select pgnr,  ST_AsGeoJSON(geom) from xxx"
    cursor.execute(postgreSQL_select_Query1)
    records = cursor.fetchall()

except (Exception, psycopg2.Error) as error:
    print("Error while connecting to PostgreSQL")

finally:
    # closing database connection.
    if (connection):
        cursor.close()
        connection.close()
        print("PostgreSQL connection is closed")

for i, row in enumerate(records):

    polygon = json.loads(row[1])

    cacheSource = r"E:\xxx_cache"
    cacheTarget = r"E:\_xx"
    cacheName = "Target"
    cacheType = "TILE_CACHE"
    storageFormat = "Compact"
    scales = "100000"
    areaofinterest = FeatureSet.from_json(polygon)

    arcpy.ExportTileCache_management(cacheSource, cacheTarget, cacheName,
                                     cacheType, storageFormat, scales, areaofinterest)

    if i == 0:
        break


texts: How to create new attributes in a feature from the value of another attribute

texts:
I have an Excel reader in my workspace that have 2 attributes (CurrentName and NewName). The value of "CurrentName" shows a list of attribute names in the another reader that already are exist in my work space and the value of "NewName" shows a list of attribute names that must be create in the workspace.

I want to create new attributes that their name comes from "NewName" value and the value of new attributes comes from the another attribute that the name of it is similar to "CurrentName" value , for example:

I have this Excel file:

and add this file to my FME workspace as a dataset (reader1). In my Workspace I have another reader (reader2) that has "id" and "technology" attributes.
So I want to create "NAME" and "MODEL" attributes in reader2 with the exact value of "id" attribute for the value of "NAME" attribute and the value of "technology" attribute for the value of "MODEL" attribute.

How can I do that in FME? Is there any transformer to do that?


texts: I have 1000 or so points on a map, (I have GPS coordinates. I would like to use google maps to find the shortest distance using roads, between pairs of points. I would like to measure all pairs within a distance of 100km. Is this possible?


texts: Geojson + Ajax to openlayers 3 repeated loading of data

texts:
I'm loading Geojson-data (from Postgis) to Openlayers 3 via bounding box stategy. Everything looks good.

The problem ist, that the data is reloaded all the time even if no change of the boundig box occurs.

Here is my code:

var geoJsonFormat = new ol.format.GeoJSON();
var vectorSource = null;

function loaderFunction(extent, resolution, projection) {
  var extent = ol.proj.transformExtent(extent, projection.getCode(), ol.proj.get('EPSG:4326').getCode());
  var url = "./_mapdata/n2k_dh_geojson.php?bbox=" + extent.join(',');
  $.ajax({
    url: url,
    success: function (data) {
      var features = geoJsonFormat.readFeatures(data, {
        dataProjection: 'EPSG:4326',
        featureProjection: 'EPSG:3857'
      });
      vectorSource.clear(true);
      vectorSource.addFeatures(features);
    }
  });
}

vectorSource = new ol.source.Vector({
  format: geoJsonFormat,
  loader: loaderFunction,
  projection: 'EPSG:4326',
  strategy: ol.loadingstrategy.bbox
});

.... later I add the layer to a layers array

      new ol.layer.Vector({
      source: vectorSource,
      style: new ol.style.Style({
            fill: new ol.style.Fill({color:'rgba(100,250,0,0.1)'}),
                stroke: new ol.style.Stroke({
                    color: 'green',
                    width: 2
                })
            })
   })

How can I prevent, that the data is reloaded all the time?


texts: Earth Engine Plugin in QGIS to use existing vector

texts:
I success to run this script in Earth Engine Code Editor (with Java). My vector shape is shareable to you.
var pemotong = ee.FeatureCollection("users/muhammadichsan/mencoba_clip");
var dataset = ee.Image('USGS/SRTMGL1_003');
var elevation = dataset.select('elevation');
var slope = ee.Terrain.slope(elevation);
Map.setCenter(108.5550607304925137, -6.4026603600933889);
Map.addLayer(slope, {min: 0, max: 255}, 'slope');
Map.addLayer(slope.clip(pemotong));

Now, I have QGIS and EarthEnginePlugin with Python API. I want to do same to QGIS. So, I write this in built in python editor in QGIS
import json
import ee
from ee_plugin import Map

pemotong = ee.FeatureCollection("????")

dataset = ee.Image('USGS/SRTMGL1_003')
elevation = dataset.select('elevation')
slope = ee.Terrain.slope(elevation)
Map.setCenter(108.5550607304925137, -8.2129304795692413, 10)
Map.addLayer(slope.clip(pemotong))

So the question is, what is "????" if I want use existing vector layer in QGIS?.
For example, I have active layer named 'xyz'
What I have done, but still error
import ee
from ee_plugin import Map

xyz = iface.activeLayer()
pemotong = ee.FeatureCollection(xyz)

dataset = ee.Image('USGS/SRTMGL1_003')
elevation = dataset.select('elevation')
slope = ee.Terrain.slope(elevation)
Map.setCenter(108.5550607304925137, -8.2129304795692413, 10)
Map.addLayer(slope.clip(pemotong))

Output:
Traceback (most recent call last):
  File "G:\QGIS3~1.10\apps\Python37\lib\code.py", line 90, in runcode
    exec(code, self.locals)


texts: How to get 2010 Census Block Housing Unit Data for the entire US for 15-digit FIPS blocks

texts:
Possible Duplicate:
  Where can I find the finest granularity of 2010 US census population information?

I'm working on a little data munging project and I need to find 2010 US Census data for the entire country down to 15 digit FIPS numbers / census block level. I can find some summaries but nothing granular enough. Is this data out there? And how would one get it?


texts: Removing pixels in raster map using ArcGIS for Desktop?

texts:
I made a raster. I wonder how can I remove pixels with values lesser than an specific value? for example How is possible to set the value of classes from 0.92?

i have searched it and i found this answer to Removing pixel values below 0 in Digital Elevation Model in ArcGIS for Desktop?:

  You can use tool Raster Calculator. Then type this expression:

SetNull("raster_name"<0.92, "raster_name")

should i type the "SetNull" in the command?
i am really confused....


texts: I am kind of new with statistics. I have two raster with categorical values and I want to get chi-square. How can I do this using R? Rasters have the same region and  same resolution. Is chi-square the right choice for the computation of variance by ranks?


texts: Map Graticule in UTM projections (openlayers)

texts:
I'd like to know if there is a way to display the Graticule labels on a different projection like UTM.

Visit this link and search for EPSG: "31982" (UTM 22s SIRGAS2000) projection to understand the issue.

See the EPSG:3857 degrees coordinates at the image above.

This is what I want to display:

The decimal values for each line that comes from the UTM projection.

My code is here: https://codesandbox.io/embed/reprojection-by-codegraticule-74zd4?fontsize=14&hidenavigation=1&theme=dark


texts: What does graphic.resLayer mean in ArcGIS API for JavaScript

texts: There is a property called reslayer in Graphic of ArcGIS JavaScript API. What does it mean? I have gone through the documentation and there wasn't any information about it.


texts: Using QGIS to find and calculate parts of polylines within polygon?

texts:
I've already got QGIS (2.18.2) to

add an ArcGIS shapefile with polylines
create a polygon which is a buffer region around a subset of the polylines in the ArcGIS file
add a KML file with one polyline
reprojected the KML file into the same projection as the ArcGIS layer.

What I'd like to do is:

split the polyline from the KML file into two layers: one with the parts of the polyline that lie within the buffer polygon, and another with the parts that are outside the buffer polygon,
so that I can get the total lengths of the polyline within the buffer polygon and outside.

I tried to use the QGIS Geoprocessing > Intersection tool:

with the buffer polygon as Input layer and polyline as Intersect layer, I get a polygon virtual layer with a bunch of features that don't make sense, and
with the polyline as Input layer and buffer polygon as Intersect layer, I get 1005 See log for more details in the Log tab.

I'm concluding from this that the Intersect tool is meant for finding the intersection between polygons, not polylines and polygons.

I found a similar question that involves the Process Modeler, but am new to QGIS and going through the Training Manual (in lesson 12 now — looks like the Process Modeler is covered in lesson 17), so the Process Modeler approach is beyond me as I ask this question.

If there's a straightforward way to do this with tools other than QGIS+GRASS, I'd be happy to hear about alternatives, too.

Looking ahead (in case there are alternatives which have trade-offs down the line), I'd like to be able to automate this analysis so that

instead of manually comparing one KML file's polyline against one polygon to find the polyline lengths inside and outside,
I'd like to compare a set of hundreds of KML files' polylines against a couple of regions to find the lengths inside and outside the regions.


texts:
We are looking for a free or paid version of a Geocoding API that will allow us to use it for commercial purposes that is as good as google. Google's api starts at $10,000 a year which is too much. We expect to do 1000-5000 geocodes per month.

I have tried Cloudmade, but have had problems geocoding some addresses that google can do.

Does anyone have any recommendations based on experience?

Here is how I am using the Cloudmade API:

http://beta.geocoding.cloudmade.com/v3/API_KEY/api/geo.location.search.2?format=json&source=OSM&enc=UTF-8&limit=10&locale=en&q=%5Bcountry%3DUnited%20States%5D%5Bstate%3DNY%5D%5Bzip%3DZIP%5D%5Bcity%3DCITY%5DAddressRemoved


texts:
I use R (and sf, mainly) and I was wondering if it's possible to create Voronoi polygons starting from LINESTRING objects instead of POINT. For example, given a sfc such as

# package
library(sf)
#> Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3

# example linestring
my_linestring_sfc <- st_sfc(
  st_linestring(matrix(c(-2, -2, -1, -1, 0, 0), ncol = 2, byrow = TRUE)),
  st_linestring(matrix(c(2, -2, 1, -1, 0, 0), ncol = 2, byrow = TRUE)),
  st_linestring(matrix(c(0, 0, 0, 1, 0, 2), ncol = 2, byrow = TRUE))
)

par(mar = rep(0, 4))
plot(my_linestring_sfc, reset = FALSE)
plot(st_cast(my_linestring_sfc, "POINT"), add = TRUE, pch = 16)

I can create the Voronoi polygons associated with its points using the following code:

my_voronoi <- st_voronoi(do.call("c", st_cast(my_linestring_sfc, "POINT"))) %>%
  st_collection_extract()

plot(my_voronoi, reset = FALSE)
plot(my_linestring_sfc, add = TRUE)
plot(st_cast(my_linestring_sfc, "POINT"), add = TRUE, pch = 16)

Created on 2020-04-11 by the reprex package (v0.3.0)

and I can merge the polygons according to their intersection with the original linestrings. The problem is that this tessellation does not correspond to my "intended" result since the middle polygon is shared between the three lines. The ideal result should be as follows (which I created manually):

Created on 2020-04-11 by the reprex package (v0.3.0)

I read a few SO posts on this topic. This is a solution with Mathematica software while this is another solution with ArcGis, and they both use the Euclidean Allocation. In this other question they suggest "Segment Voronoi Diagrams" and CGAL. This looks like the most promising solution using PostGIS, but I have no idea how to translate that code into R. Is it possible to adapt those ideas to my problem?

Note: Check here since maybe it's possible to adapt etiennebr approach to the previous PostGIS answer.


texts: CRS problem between QGIS and fiona / OGR

texts:
I have a really annoying problem working with QGIS for visualization and python (fiona / OGR) for the actual operations.
I am, by example, doing a simple dissolving operation with fiona.
The input is a EPSG : 26931 defined by :

  +proj=omerc +lat_0=57 +lonc=-133.6666666666667 +alpha=323.1301023611111 +k=0.9999 +x_0=5000000 +y_0=-5000000 +gamma=323.1301023611111 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

But then my output systematically switch to "Generated CRS" 10000X defined by :

  +proj=omerc +lat_0=57 +lonc=-133.666666667 +alpha=323.130102361 +k=0.9999 +x_0=5000000 +y_0=-5000000 +no_uoff +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

So obviously fiona ignore the gamma value. My data being really huge I would like to avoid to correct it with gdalwarp every time as it is really long.

Here is an example of code that modify the CRS :
def dissolve(base_watershed):

from shapely.geometry import shape, mapping
from shapely.ops import unary_union
import fiona ,sys, os
import itertools
from osgeo import ogr
watersheds = "watersheds.shp"
#from http://gis.stackexchange.com/questions/149959/dissolve-polygons-based-on-attributes-with-python-shapely-fiona
with fiona.open(base_watershed) as input:
    # preserve the schema of the original shapefile, including the crs
    meta = input.meta
    with fiona.open(watersheds, 'w', **meta) as output:
        # groupby clusters consecutive elements of an iterable which have the same key so you must first sort the features by the 'HUC_8' field
        e = sorted(input, key=lambda k: k['properties']['HUC_8'])
        # group by the 'HUC_8' field
        for key, group in itertools.groupby(e, key=lambda x:x['properties']['HUC_8']):
            properties, geom = zip(*[(feature['properties'],shape(feature['geometry'])) for feature in group])
            # write the feature, computing the unary_union of the elements in the group with the properties of the first element in the group
            output.write({'geometry': mapping(unary_union(geom)), 'properties': properties[0]})

Some other part of my code using OGR does the same thing.
Is there an easy solution or do I really have to do the full process and reproject everything at the end?
I am wondering if overwriting the 10000X CRS definition by the correct 26931 would be an issue since it looks like there is no transformation between the two.


texts: Monte Carlo test of Space-Time clustering in R

texts:
I have several point patterns representing the distribution of territories of a bird species. Each point pattern represents the distribution of one year. I already succeeded in computing Ripley's L function for each year seperately, but now I would like to test whether there is space-time clustering present over all years.

Some research leads me to a Monte-Carlo test of space-time clustering which should be possible with the R-package SPLANCS.
I did not manage to find some examples on this (with R-scripts) and the manuals (http://cran.r-project.org/web/packages/splancs/splancs.pdf) are too advanced for me. Is there anyone who can provide an example, or can explain how exactly my data should look like before I read them in in R?

I assume that my points for all years should be put together in one shapefile. But how do I link that with a 'vector of times'. And most of all, what does a vector of times looks like?


texts: Get XY in different projection without changing data frame using python addin

texts:
I been searching for a way to use this http://forums.arcgis.com/threads/49273-getting-xy-from-mouse-click-using-python tool but print the XY in different coordinate system like WGS 1984 (decimal degrees). I want to do this without changing my data frame and keep it in its original stateplane system. I have created a python addin to display the coordinates in DD but this required the data frame to convert to WGS 84 and back again to its original projection (stateplane). This option does work but cannot be used when editing because it locks the mxd. So what is the python code to find the XY in DD from the mouse click regardless of what projection the data frame is set?

here is what i have so far:

class Sel_Property(object):
    """Implementation for Notifications_addin.selectproperty (Tool)"""
    def __init__(self):
        self.enabled = True
        self.shape =3
    def onMouseDownMap(self, x, y, button, shift):
        mxd = arcpy.mapping.MapDocument("Current")
        dataframe_string = mxd.activeView
        dataframe= arcpy.mapping.ListDataFrames(mxd,dataframe_string)[0]
        oldspatialref = dataframe.spatialReference
        spatialref = arcpy.SpatialReference(4326)
        dataframe.spatialReference = spatialref
        dataframe.spatialReference = oldspatialref
        pythonaddins.MessageBox("Your mouse clicked at " + str(x) + " , " + str(y),"My Coordinates:")

Here is the new code if anyone is interested.

    mxd = arcpy.mapping.MapDocument("CURRENT")
    pt = arcpy.PointGeometry(arcpy.Point(x, y), mxd.activeDataFrame.spatialReference)
    Point = pt.projectAs(arcpy.SpatialReference(4326))
    partObj = Point.getPart(0)
    partObj.X
    partObj.Y
    pythonaddins.MessageBox("Your mouse clicked at " + str(partObj.X) + " , " + str(partObj.Y),"My Coordinates:")


texts: Function to change parts of data source instead of replacing path?

texts: Using local files on a remote desktop server is giving me some issues. Using a locally made mxd in ArcView on the remote desktop will keep the C:... start to the filepath, but the server accesses my local files with a C$... start to the file path. Is there an arcpy function that allows me to just replace the ":" with a "$"? I can only seem to find functions that entirely replace the path.


texts: Calculating azimuth for shortest line from point to polyline in QGIS

texts:
How can I calculate the azimuth for the shortest line from a point (start) to a polyline (end) in QGIS?

NNJoin calculates only the distance and joins the attribute from the line to the point attributes. Does QGIS have a function or plugIn for this?


texts: Uploading mxd to ArcGIS Online with Python?

texts:
I need to create a Python script to upload an MXD to AGOL. I have created an script in arcpy to create a service definition file (*.sd) from the mxd. I have not been able to upload the sd to AGOL unless ArcMap is open. I use the following code:

arcpy.SignInToPortal_server("AGOL_USER", "AGOL_PASS", "http://www.arcgis.com/")
arcpy.UploadServiceDefinition_server(sd, "My Hosted Services", service)

I get the following error:

arcgisscripting.ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000732: Server: Dataset My Hosted Services does not exist or is not suppoted
Failed to execute (UploadServiceDefinition).

I have been able to upload the service document file using ArcGIS API for Python:

gis = GIS("https://XXXX.maps.arcgis.com", "AGOL_USER", password)
item = gis.content.add({}, "XXX.sd")
published_item = item.publish(overwrite=True)

a) How can I upload the service definition file to AGOL using arcpy? The script needs to be run scheduled, so I cannot have ArcMap open

b) Would be posible to merge arcpy code with ArcGIS API for Python?


texts:
I have polygon featureclass (BookmobileLocations) which has one field (Marker, data type Text). The field Marker has values (Site 1, Site 2, Site 3...). I want to do query using Select_Analysis.The first code works fine. The second code has problem. It should be the where_clause (SQL expression) issue. What's wrong with the second one "Marker = \'+siteName+\'". Where to search right SQL expression rules?

siteCursor = arcpy.da.SearchCursor("BookmobileLocations","Marker")
for row in siteCursor:
    siteName = row[0]
    arcpy.Select_analysis("BookmobileLocations",r"C:\EsriPress\GISTPython\MyExercises\Scratch\Temporary Storage.gdb\SiteTemp3",'"Marker" = \'' + siteName + "\'")

siteCursor = arcpy.da.SearchCursor("BookmobileLocations","Marker")

for row in siteCursor:
     siteName = row[0]
     arcpy.Select_analysis("BookmobileLocations",r"C:\EsriPress\GISTPython\MyExercises\Scratch\Temporary Storage.gdb\SiteTemp3","Marker = \'+siteName+\'")


texts: How to add marker dynamically when clicked on map using OpenLayers 3?

texts:
I'm trying to add point (marker) on map click with custom style using OpenLayers 3. I tried this code, but nothing happens when user click on the map:

  var map = new ol.Map({
            layers: [
          new ol.layer.Tile({
              source: new ol.source.MapQuest({layer: 'sat'})
          }),
          new ol.layer.Vector({
              source: vectorSource
          })
      ],
            target: 'map',
            view: new ol.View({
              center: [0, 0],
              zoom: 2
            })
          });

    var styles = {
            'square': new ol.style.Style({
                image: new ol.style.RegularShape({
                    fill: fillColor,
                    stroke: stroke,
                    points: 4,
                    radius: 10,
                    angle: Math.PI / 4
                })
            }),
            'triangle': new ol.style.Style({
                image: new ol.style.RegularShape({
                    fill: fillColor,
                    stroke: stroke,
                    points: 3,
                    radius: 10,
                    rotation: Math.PI / 4,
                    angle: 0
                })
            })
    };

    var stroke,
        fillColor;

    map.on("click", function(evt) {

        stroke = new ol.style.Stroke({color: 'black', width: 1});

        var coords = evt.coordinate;

        var addedMarker =  new ol.geom.Point( [
        ol.proj.transform(coords, 'EPSG:4326', 'EPSG:3857')]);
            fillColor = 'blue';

         var featurething = new ol.Feature({
                name: "Marker 01",
                geometry: addedMarker
            });

            featurething.setStyle(styles['triangle']);

            vectorSource.addFeature( featurething );

    });


texts: Export composer to single PDF file using PyQGIS

texts:
I am using PyQGIS to create automated maps.

My script exports atlas features to new JPG images. But I need now to export my features to single PDF file all my features like option where have qi's composer in atlas define.

Here to code for export to my features:

for i in range(0, myAtlas.numFeatures()):
    myAtlas.prepareForFeature(i)
    jobs = r"C:/Users/username/Desktop/test/"
    output_jpeg = jobs + str(i)+ "_BMS_plan.jpg"
    myImage = myComposition.printPageAsRaster(0)
    myImage.save(output_jpeg)
myAtlas.endRender()

How should I change this code?


texts: Creating parallel lines in Geographic Coordinate System using ArcPy?

texts:
I have a polyline feature class that is using GCS_WGS_1984 (code 4326). I am trying to create lines that are parallel to a polyline from this feature class at distance d.
From some searching, I've only been able to find out how to create parallel lines here:
Creating parallel line in ArcGIS Pro using ArcPy?
but it only works for 2D:
import arcpy, math
infc=r'..\SCRARCH\clone.shp'

def CopyParallel(plyP,sLength):
    part=plyP.getPart(0)
    lArray=arcpy.Array();rArray=arcpy.Array()
    for ptX in part:
        dL=plyP.measureOnLine(ptX)
        ptX0=plyP.positionAlongLine (dL-0.01).firstPoint
        ptX1=plyP.positionAlongLine (dL+0.01).firstPoint
        dX=float(ptX1.X)-float(ptX0.X)
        dY=float(ptX1.Y)-float(ptX0.Y)
        lenV=math.hypot(dX,dY)
        sX=-dY*sLength/lenV;sY=dX*sLength/lenV
        leftP=arcpy.Point(ptX.X+sX,ptX.Y+sY)
        lArray.add(leftP)
        rightP=arcpy.Point(ptX.X-sX, ptX.Y-sY)
        rArray.add(rightP)
    array = arcpy.Array([lArray, rArray])
    section=arcpy.Polyline(array)
    return section

with arcpy.da.UpdateCursor(infc,("Shape@","Width")) as cursor:
    for shp,w in cursor:
        twoLines=CopyParallel(shp,w)
        cursor.updateRow((twoLines,w))

I have tried first projecting to spatial reference code 3857, then running the above code, then converting back. However, there is some distortion (~8% error).
Is there a solution that would work for geographic coordinate systems?


texts: Is it possible to label pie/bar charts in MapInfo Pro v10.0?

texts:
What I have in mind looks something like this...

or like this...

The idea is to get values/percentages shown next to the pie/bar chart.

If it is not possible to do it in MapInfo, is there some other GIS software that has this option?


texts: I have a fairly simply question about something I just don't know how to do. I am aware you are able to simplify polygons in QGIS based on distance tolerance but my questions is, is there a way to simplify based on a set maximum number of desired points? For example, if the number 1000 was entered, the polygon would be trimmed to 1000 points regardless of any distance thresholds. Does anyone know of any ways to do this?


texts: Adding feature class to TOC using ArcPy?

texts:
This question has been asked multiple times with proper answers:

How to add feature class to MXD with ArcPy (Python)?
Adding shapefile or feature class as layer in ArcGIS Desktop using Python/ArcPy?
Adding feature class into ArcGIS 10.3 mxd using python script

but I still cannot get it right.

I am using ArcGIS 10.2.2.

I want to add via a Python standalone script a Featureclass from a file-gdb to the TOC.
If I try it in the ArcGIS Python window it works fine. However in a standalone script nothing happens. There is no error message or anything. I already thought it might be a bug, but it seems to work for everyone else.
I also tried it with a .lyr and it didn´t work either. Someone got a hint for me?

Here is the code (I also tried the one from the ArcGIS help site, and from all the links above).

import arcpy
mxd = arcpy.mapping.MapDocument(r"C:\GIS\KUP\Check.mxd")
df = arcpy.mapping.ListDataFrames(mxd, "*")[0]
addLayer = arcpy.mapping.Layer(r"C:\GIS\KUP\Ergebnisse.gdb\Potenziale_singlepart_final")
arcpy.mapping.AddLayer(df, addLayer, "BOTTOM")

I also tried :

addLayer = arcpy.mapping.Layer(r"C:\GIS\KUP\Potenziale_singlepart_final.lyr")

What am I doing wrong? I am trying for a week now. I`m getting crazy with this one. Seems so easy.


texts: Connecting Unconnected Multiple Vertical Lines in QGIS

texts:
So I have these multiple vertical lines as shown below.

So what I am aiming for is to connecting these unconnected lines into one lines automatically. What I've been done so far is to connected it manually using join multiple lines plugin. The result is given below,

When using this plugin I can only select one vertical lines that I want to merge, I can't really like select all the unconnected lines and running the plugin because it will give me result like this.

any idea as to how to do this task automatically using QGIS ?


texts: Zonal statistics tool excluding some features in ArcMap?

texts:
In Arcmap 10.6, while using Zonal Statistics as Table tool, some of my feature zone data is being excluded from my output table.
Here is the scenario. I'm using a polygon layer of buildings, with a raster layer above it that represents radio coverage.  Each pixel is an integer value representing signal strength at that area.  My intention is to add attributes to the building polygons that represent the amount of signal power in their are. I'm using Zonal Statistics to Table tool to with buildings as my feature zone data and FID as the Zone field.  I'm using radio coverage raster as my Input Value raster.  The output is a table that includes a row for each polygon/building, ( except those that are inexplicably missing  ), and each row has a value for based on the raster signal strength values in that polygon's area.  As an example, the polygon in light blue is one of the polygon/buildings being left out.  It's well within raster coverage as can be seen by the Pixel Inspector tool.  Also the pixels are smaller that the building, 10' x 10'. The surrounding features are being included.  Out of 48400 polygon/buildings, about 80 seemingly randomly located building are being excluded from the output.
Can someone please tell me why?
Also as a side note. I originally had this problem using 1100' x 1100' pixels.  I read somewhere that the pixels should be smaller than the polygons.
Why is that?
If a single pixel covers several polygons, shouldn't all the polygons inherit the value of the pixel that's completely covering it?


texts: Combining several polygon objects in one polygon object with GeoTools

texts:
I have several polygon objects and I need to "merge (or combine or dissolve)" them toghether into one polygon object. With a GIS tool like QGIS, I do this with processing option "Combine selected spatial objects" and the result is what I expect but doing programmatically with Java and GeoTools library, I'm not able to get this behaviour.

To give you an idea, I have two polygons in one layer like this:

Each square is an independent polygon in my shapefile (an independent object with its attributes).

Now I want to transform them into one object. I tried with "union" method of Geometry class in this way:

    //polygon1 and polygon2 are Polygon objects (the squares in the above image)
    Geometry combine = polygon1.union(polygon2);

And this is the result when I visualize the union in the GIS:

As you can see, a line is added between the two figures (although the union is made succesfully and I have only one entity as I wish) but I want this without the line connecting them.  The polygons are well formed because by doing with QGIS the operation is made right.

Anyone know if is it possible?

I'm looking for some time and I'm concluding that maybe my approach to achieve this is bad.


texts: Add UTM?Xy /Lat long to attribute table for point data in QGIS

texts: I am running V QGIS 3.2.3.  I have found some answers on here how to do this but in previous versions. Is there a simple way to calculate/Add Xy or utm coordinates to the attribute table?


texts:
I want to create a map that shows the catchment area of download mirrors. What I found so far is only a manual option where you have to draw an area around that point. But I don't know exactly where the boundaries are for the mirror, the algorithm always chooses the nearest one.
I created a draft, it would probably look something like this:

I have not found a good way to create a map like this just by entering lat/long or the addresses of the servers and not manually drawing the affected areas.

Are there tools already out there that I missed?

Or how could I get at least some points of the border to draw the area?


texts: Full Easting and Northing coordinates values cannot be encoded in GRASS r.water.outlet in SEXTANTE-QGIS

texts: The Easting or Northing coordinate parameter accepts only 5 digits.  The Easting and Northing values for my outlet are 6 and 7 digits long (integer).  How can they be encoded in r.water.outlet?  Can real numbers be encoded as well?


texts:
I am running several WMS layers from a GeoServer (2.2 RC3 at the moment of writing). All layers are ImagePyramids with sizes ranging from 100GB to 1TB in size.

I have read the GeoServer on Steroids several times and I'm probably missing half of it still. But there's a few things I'm wondering about.

My pyramids consists of 9-10 levels of GeoTIFF files, each level having 1/2 resolution of the previous. Each GeoTIFF has a file size with an average at around 8MB, a defined CRS and a size aspect ratio of 1:1.

Is my TIFF file size inefficient?
Will I get a performance boots by adding overview to my TIFF's?
Will inner tiling increase performance?
What could help me increase performance?
What have you guys learned from your experiments with ImagePyramids in GeoServer?

1,3) In my view these two comes down to how a TIFF file is being handled in general and especially by GeoServer and I'm clueless here.

2) If I added overviews to my TIFF files my guess is that it corresponds to jump a level in the pyramid and thus I would only waste disk space?

Edit: Changed pyramid resolution from (the wrong) 1/4 to (the correct) 1/2. The pyramids are following the powers of 2 guidelines.


texts: Calculating distance of individual points along line in PyQGIS?

texts:
I have a vector layer and need the distance between the points. This is an intermediate step for my further calculation and programming with Python. The starting point is stored in another layer (brown point). The order along the line is not known and must be determined from the distances between the points. Only the starting point is known and the order depending on the length of the point to the starting point must be determined. This step is still missing.

How exactly can I calculate the length of the points along the line with Python or QGIS and store it in an attribute table?


texts: Draw graph in QGis on selection (of multiple points)

texts:
The issue is that I have a geopackage with some points. One of the attribute columns stores the date information. My rather naive question would be, without actually thinking about any kind of implementation, if it is possible to interactively select some points and draw some kind of graph to visualize quickly the dates of those points?
I know there are the Qgis Actions and I already some examples a while ago. But can I somehow combine this with the current selection (of potentially multiple points) ?


texts:
Would like to start my leaflet map without any features displayed on start. The feature should not be displayed before a location is found via my search.control.

I started with

// first "invisible" icon
var smallIcon = new L.Icon({
    iconUrl: 'mypicture',
    iconRetinaUrl: 'mypicture',
    iconSize:    [0.00, 0.00],
    iconAnchor:  [0, 0],
    popupAnchor: [1, -34],
    shadowSize:  [0, 0]
});

// second "visible" icon
var smallIcon2 = new L.Icon({
    iconUrl: 'mypicture',
    iconRetinaUrl: 'mypicture',
    iconSize:    [25, 41],
    iconAnchor:  [12, 41],
    popupAnchor: [1, -34],
    shadowUrl: 'mypicture',
    shadowSize:  [41, 41]
});

// adding geojson and display feautures with "invisible" marker
var marker;
var json_AlleJSON = new L.geoJson(myjson, {
    onEachFeature: pop_All,
    pointToLayer: function (feature, latlng) {
        marker = new L.marker(latlng, {icon: smallIcon});
        return marker;
    }
});

And to change the icon when location is found:

// changing the icon to "visible"
searchControl.on('search_locationfound', function(e) {
    sidebar.hide();
    e.layer.setIcon(smallIcon2);
});

This actually works BUT the problem I got: even by reducing the iconSize to 0.0, the map interacting gets so slow because of loading an image for every feature (my json contains about 100 entries).

I tried to solve this by leaving the iconUrl empty but this didn't work.

Any ideas?


texts:
What can I do with the given information: I don't know where I am in the World having some particular tile. It's only an ocean (or sea).

"https://services.arcgisonline.com/ArcGIS/rest/services/World_Topo_Map/MapServer/tile/14/5210/8141"

I know, the viewing scale is 36111.909643 and current zoom is 14. Am I able, using the naming convention to guess what would be it's absolute coordinates? I'm not sure about the projection but may choose something later if only I knew what's the region of these waters. It's a guess game.

How to take this approach? How to calculate it? Here I have info about the systematics of this "World_Topo_Map" service.

https://services.arcgisonline.com/arcgis/rest/services/World_Topo_Map/MapServer

How can I get to lon/lat for center of a particular tile?

Here's similar for google: Calculate lat lon bounds for individual tile generated from gdal2tiles
but where to put the tile IDs in that formula?

EDIT

I've found a formula for R (and many other languages) but how can I be sure the precision is correct for this particular service? Looking at the globe it's quite there.

# tile formula
rad2deg <- function(rad) {(rad * 180) / (pi)}
num2deg <- function(xtile, ytile, zoom) {
  n = 2.0 ** zoom
  lon_deg = xtile / n * 360.0 - 180.0
  lat_rad = atan(sinh(pi * (1 - 2 * ytile / n)))
  lat_deg = rad2deg(lat_rad)
  return (list(lat=lat_deg, lon=lon_deg))
}

num2deg(8141,5210,14)

I mean I'd need to have a chance of including the projection Spatial Reference: 102100  (3857) for this service in some way.


texts: Syncing database with ArcMap and ArcSDE?

texts:
I am using SQL in Oracle to fill and empty feature class table, basically something like this:

INSERT INTO DAMO_W.BRUG (OBJECTID,CODE,NAAM,GLOBALID,SHAPE,CREATED_USER,CREATED_DATE,LAST_EDITED_DATE,LAST_EDITED_USER)
SELECT  OBJECTID, CODE, NAAM, GLOBALID,SHAPE,CREATED_USER,CREATED_DATE,LAST_EDITED_DATE,LAST_EDITED_USER
FROM WS_BRUG;

The SQL runs and the table in the database gets filled, but in arcMAP when opening the feature it is totally empty, therefore there is out of sync between ArcSDE and Database.

I have restarted  ArcGIS connection and ArcMAP and no success, probably there is some trigger on ArcSDE or metadata that needs to be updated. Any tips on what could I do?

ArcGIS version: 10.6.1
Database: Oracle


texts: Apply a definition query to a Layer Group

texts:
I'm trying to apply a definition query to multiple layers within a Layer Group by using the following script where I would simply change the values for my 'SUBCODE' field:

import arcpy

#Variables to form defintion query
field = '"SUBCODE"'
values = "'1152','1153'"
#concatenate query syntax
queryStr = str(field) + "=" + str(value)
#Specify the aprx project (CURRENT), dataframe (Layers)
p = arcpy.mp.ArcGISProject("CURRENT")
m = p.listMaps("Map_3D")[0]
#Apply defintion query to specified layer group
for lyr in m.listLayers( "3D Layers")[0]:
    if lyr.supports("DEFINITIONQUERY"):
    lyr.definitionQuery = queryStr
arcpy.RefreshActiveView()
del aprx

This is the error I'm receiving when I run this in Arcpy:

Traceback (most recent call last):
  File "<string>", line 12, in <module>
IndexError: list index out of range

There's a similar question posted on SE that applies the same logic of applying a definition query to a Layer Group, however I'm using ArcGIS Pro and the post refers to ArcMap.  There are some migration from arcpy.mapping to ArcGIS Pro I'm having trouble with - http://pro.arcgis.com/en/pro-app/arcpy/mapping/migratingfrom10xarcpymapping.htm

Here's an illustration of my Layer Group:


texts: WFS from geosever very slow to load in openLayers3

texts:
I originally posted this question about trying to load a geoJSON with openLayers3: geojson very slow to load in openLayers3.

Based on a suggestion given in the post I'm now loading the layer in as a WFS from geoserver. However the response time is still really slow. It will eventually load in firefox but won't even load in Chrome - I get the 'He's dead Jim' message. I'm wondering about limiting the zoom level at which the data draws (only once you get into to a fairly zoomed in view), and also what else I might be able to do to improve the response time. Here is my code:

sourceVector = new ol.source.Vector({
    loader: function(extent) {
        $.ajax('http://52.24.37.126:8080/geoserver/wfs',{
            type: 'GET',
            data: {
                service: 'WFS',
                version: '1.1.0',
                request: 'GetFeature',
                typename: 'Iada_workspace:parcels',
                srsname: 'EPSG:3857',
                outputFormat: 'application/json',
                bbox: extent.join(',') + ',EPSG:3857'
                },
            }).done(loadFeatures);
        },
        strategy: ol.loadingstrategy.tile(new ol.tilegrid.createXYZ({
            maxZoom: 19
            })),
    });

window.loadFeatures = function(response) {
    geoJSON = new ol.format.GeoJSON();
    sourceVector.addFeatures(geoJSON.readFeatures(response));
    };

var layerVector = new ol.layer.Vector({
    source: sourceVector,
    style: new ol.style.Style({
        stroke: new ol.style.Stroke({
            color: 'rgba(0, 0, 255, 1.0)',
            width: 2
            })
        })
    });

var map = new ol.Map({
    layers: [new ol.layer.Tile({source: new ol.source.OSM()}), layerVector],
    target: document.getElementById('map'),
    view: new ol.View({
                //projection:projection
                center: ol.proj.transform(
                    [-116, 42], 'EPSG:4326', 'EPSG:3857'),
                zoom: 6
        })
    });


texts: Adding per-row creation time and modification time in GeoPackage geometry table?

texts:
I need to add ctime and mtime to a GPKG sqlite database in table with geometry.

Apparently logic solution: add ctime column with DEFAULT CURRENT_TIMESTAMP and add trigger ON UPDATE to set mtime = CURRENT_TIMESTAMP.

Problem: this works perfectly for non-geometric tables.
With geometric tables QGIS complains "Warning 1: geometry column '[ctime|mtime]' of type 'TIMESTAMP' ignored" and moreover when I add a geometric feature using qgis the trigger on mtime seems ignored (probably due to gpkg driver)

Is there a correct way to perform this at db level?


texts: Summing numeric attributes of field during Spatial Join using ArcPy?

texts:
Just performed a spatial join operation on a polygon layer (made of 500 grids cells) and point layer (1 million points). The aim was to find the sum of a  field in the points layer (field AREA) during the spatial join process. However, the result of the spatial join output shows 0 for all the rows. Below is the code used.

import arcpy

arcpy.env.overwriteOutput = True

point = r"C:\Example\Scratch.gdb\trees"
grid = r"C:\Example\Scratch.gdb\grid"
grid_join = r"C:\Example\Scratch.gdb\grid_join"

#Make feature layers
arcpy.MakeFeatureLayer_management(point, "point_lyr")
arcpy.MakeFeatureLayer_management(grid, "grid_lyr")

arcpy.SelectLayerByAttribute_management("point_lyr", "NEW_SELECTION", """ "TYPE" = 'Mangrove' """)
arcpy.FeatureClassToFeatureClass_conversion("point_lyr", "in_memory", "tree_new")
tree_new = r'in_memory\tree_new'

arcpy.SpatialJoin_analysis(grid, tree_new, grid_join, "JOIN_ONE_TO_ONE", "KEEP_ALL", "AREA \"AREA\" true true true 8 Double 6 7 ,Sum,#,C:\\Example\\Scratch.gdb\\in_memory\\tree_new,AREA,-1,-1", "COMPLETELY_CONTAINS", "", "")


texts: Zoom to feature by clicking cell in the Attribute Table

texts:
I have a layer which has FID field representing the unique ID of features. In QGIS, there is a Zoom to selected rows button on top of the Attribute Table. It has a shortcut (Ctrl + J) and it doesn't work (at least in my computer). (In ArcGIS, when you double-click a row number in the Attribute table, it zooms to related feature)

Sometimes, it's a waste of time to click the zoom button if considering hundreds of features and I need to speed up the zooming process to save time.

I can display the attribute table using the lines below and get the reference of it.

lyr = QgsProject.instance().mapLayersByName("LAYER_NAME")[0]
attr_table = iface.showAttributeTable(lyr)

Is there any programmatic way to zoom to a feature by clicking FID cell in the Attribute Table in QGIS? How can I add this capability to the Attribute Table?


texts:
Does the pgr_astar heuristic function of postgreSQL verify that the proposed path is optimal on time and length despite what mesure we choose as cost on the query?

Example:
When asking this query

SELECT a.seq, a.id1 AS node, a.id2 AS edge, b.source, b.target, b.cost, b.reverse_cost, b.geom_way
 FROM pgr_astar
('SELECT id, source, target, fuel as cost, x1, y1, x2, y2, reverse_cost FROM San_francisco', 8,895, true, true)as a LEFT JOIN San_francisco as b
ON a.id2 = b.id
order by seq

Does pgr_astar verify that the result is optimal on length and time in addition of the cost(fuel)?


texts: ArcGIS Online: Adding Screenshots to my content items

texts:
Within ArcGIS Online, is it possible to associate an array of screenshots with a content item?

I see a property in the ArcGiS Online Search API response that includes "screenshots":[]


texts: Using Mapfish Print v2 yaml file for Mapfish Print v3


texts:
I have az update cursor loop, in which I'd like to use the function written earlier.

Function:
def TL ():
    addlayer = arcpy.mapping.Layer(r"gk25.shp")
    arcpy.mapping.AddLayer(df, addlayer, "AUTO_ARRANGE")
    arcpy.SelectLayerByLocation_management("gk25", 'intersect', "ActualLayer")
    cursor = arcpy.SearchCursor("gk25") # cursor loops through feature class
    for row in cursor:
        rows = arcpy.SearchCursor("gk25") # defines rows variable, which comes from SearchCursor function
        fields = arcpy.ListFields("gk25", "", "ALL") # defines fields variable as a list, from all the fields in the table
        for row in rows: # loops through all rows
            for field in fields: # loops through all fields
                if field.name == "T_LAP": # gets the number of the last registered object from joined table
                        tlp = (row.getValue(field.name))
                        arcpy.AddMessage(tlp)
                        return tlp
    arcpy.Delete_management("gk25") # deletes the "gk25" layer from the map document

And here is the update cursor:

fc = "Kutak" # defines fc variable as the "Kutak" layer
    fields = arcpy.ListFields(fc, "", "") # gets all fields to a list for the loop
    for field in fields: # loops through all fields
        with arcpy.da.UpdateCursor(fc, field.name) as cursor: # loops through the fields
            for row in cursor:
                if field.name == fieldName10:
                    row [0] = TL
                    cursor.updateRow(row)
                    arcpy.AddMessage(TL)

I get this error message:

TypeError: value #0 - unsupported type: function

Any ideas to get rid with this?


texts:
I am doing indoor multi-floor routing using pgRouting.

At the moment, I have a dummy network system created to analyze shortest routes between locations with the results served as a WMS layer in GeoServer.  My initial test of the system only included a simple network created on the same z-plane, and everything works perfectly.

At that point, I created some mock interior floor plans to create indoor networks on two separate floors.  For simplicity, I simply assigned the z-elevation of the bottom floor as ‘0’ and the upper floor as ‘10’.  I then created the vertical corridors (e.g. elevators and stairs) in ArcScene to make sure the vertical connections existed, resulting in a fully-connected 3D network.

Using the SQL code found here https://github.com/mapcentia/geocloud2/blob/master/app/scripts/sql/pgr_createTopology3D.sql, I created the modified pgr_createTopology3d and pgr_pointToIdZ functions to build topology on the 3D network.  Following that, I followed normal procedure for building the SQL View to serve the queried results as a layer to GeoServer.

Now that I have the 3D network in the routing system, I tried running some test queries across different floors to ensure that it routes appropriately; unfortunately it does not.  Instead, the system continues to operate on the same z-plane, and seems to ignore the z-values in the vertices, meaning a route is returned, but only on the plane with z-levels assigned as ‘0’.  I’m guessing this has to do with the SQL Statement that’s entered when creating the SQL View in GeoServer.  Right now, the statement is as follows:

SELECT ST_MakeLine(route.geom) FROM (
    SELECT geom FROM pgr_fromatob('ways', %x1%, %y1%, %x2%, %y2%
) ORDER BY seq) AS route

That statement above is passing to/from lat/longs to the table ‘ways’ to produce the resulting shortest path; however, there’s nothing in there to give the query any sort of z-information.

Is that where things are going wrong and not correctly routing across multiple floors?

If so, can anyone provide insight on how to correct the SQL Statement to consider z-values?

I’m no SQL expert and been trying a number of different modifications, but nothing working as of yet.


texts:
I am following a tutorial to install pgrouting. However, I didn't get the part when it says

  Start pgAdmin and create a new database based on your PostGIS template

Where can I get this postgis template? I know I have installed PostGIS using homebrew. But I couldn't find the template. I am running macos now.


texts: Error in Mapnik installation: "%1 is not a valid Win32 application."

texts:
While installing Mapnik in Windows 10 64bit and with Python version 2.7.5(64bit) I have encountered this error:

  import mapnik Traceback (most recent call last):   File "",
  line 1, in    File

  "C:\Python27\lib\site-packages\mapnik__init__.py", line 69, in

      from _mapnik import * ImportError: DLL load failed: %1 is not a valid Win32 application.

Any suggestion?


texts: More stable spatial sort using Shapely's within?

texts:
I was trying to use a modified version of the spatial sorting method illustrated in shapely's documentation, using only polygons exteriors (so that I can detect if a polygon lies inside another polygon hole), and it didn't always work.

To investigate this, I went back to the documentation example:

from shapely.geometry import Point, Polygon

class Within(object):
    def __init__(self, o):
        self.o = o

    def __lt__(self, other):
        return self.o.within(other.o)

# a is inside b that is inside c; d is far away
a = Point(2, 2)
b = Polygon([[1, 1], [1, 3], [3, 3], [3, 1]])
c = Polygon([[0, 0], [0, 4], [4, 4], [4, 0]])
d = Point(-1, -1)

sorted_feat = [d, c, c, b, a]
assert sorted_feat == sorted([c, a, d, b, c], key=Within, reverse=True)

so far so good.
But when I tried to change the order of the initial list, the result changed:

shuffled = [b, a, c, d, c]
# I was expecting the same sorted list...
assert sorted_feat != sorted(shuffled, key=Within, reverse=True)
# ... but I got this instead.
assert [c, b, a, d, c] == (sorted(shuffled, key=Within, reverse=True))

I understand that this is because of how sorted() works, that is guaranteed to be stable (items with the same key appears in the original order) and that the within() method obviously returns false if the two geometry are apart from each other, but I was expecting that at least the "relative order" would be kept (that is, to always have [c, c, b, a] order and d randomly inserted).

Is there a solution to this? How can I group the geometries and get the order that I need?


texts:
In QGis 2.x, I could rotate a point symbol towards the nearest geometry of a certain kind by using:

geomnearest('featuretype', 'BEARING')

This however does not work in QGIS 3.x, and I cannot find out what rotation command to use.


texts: Feature Coordinate Transform & Translate? "copy"

texts:
Using ArcGIS Pro, I have a 2D layer in California. The feature is a line-split-at-vertices w/ Z points.

I exported the feature to a .MAPX file and then imported into the Mt. Egmont (New Zealand) 3D tutorial.

How do I copy the feature from the imported California layer over to some location on Mt. Egmont?

I've been trying using copy/paste but keep getting "Paste Features failed".


texts: Search and replace in print composer of QGIS

texts:
I made a print layout with servales pages with the QGIS print composer and added some text to the layout, like the project name etc.
Is there a fast way to change the project name in all pages?
Like a "search and replace" or do I have to change it manually in all pages?


texts: Transferring from Windows to Linux OS for ESRI ArcReader

texts:
I have a .pmf file and I would like to open it under ubuntu. However all the layers are directing to a C:/ source. As Linux OS does not have such drive, how can I change the layer source path? Thanks!

Edit: I only have ArcReader with me..... Is it a must that I need to use a editor to output a new .pmf again for linux OS?


texts: Convert Shape Length Area Data to Latitude Longitude

texts:
I am attempting to import publicly available shapefile data from http://gis.utah.gov/data/boundaries/zip-codes/ into Tableau for a data visualization project.

Unfortunately, the .shp file and .gdb database that are available for download include only the shape length / area for each of the Utah zip codes.

My goal is to convert these from the current NAD 1983 UTM 12N (EPSG: 26912) format into latitude / longitude coordinates that can be read into Tableau.

I am using Quantum GIS (QGIS 2.0.1) and would appreciate any and all help.


texts: Calculating areal distortion outside UTM zone?

texts:
One of my colleagues is working with data that is spread across two UTM zone. The majority of the data is in one zone, with a few outliers in another zone. He would like to know what the area distortion of those outliers would be if they were in the main UTM zone.

Is there a formula to calculate the areal distortion knowing how far into the other UTM zone the features were?


texts:
I am trying to create a simple ArcGIS Server SOAP consumer by C# like below snippet  but I am getting error on NAService_NAServeras illustrated in attached image. Can you please let me know which namespace I have to use?

I am already using this namespaces

using ESRI.ArcGIS.esriSystem;
using ESRI.ArcGIS.Server;
using ESRI.ArcGIS.Geometry;
using ESRI.ArcGIS.Geodatabase;
using ESRI.ArcGIS.Carto;

but they are not doing the job!

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using ESRI.ArcGIS.esriSystem;
using ESRI.ArcGIS.Server;
using ESRI.ArcGIS.Geometry;
using ESRI.ArcGIS.Geodatabase;
using ESRI.ArcGIS.Carto;

namespace ConsoleApplication1
{
    class Program
    {
        static void Main(string[] args)
        {
            private void SolveRoute()
        {
            // Connect to server
            NAService_NAServer naService = new NAService_NAServer();
            naService.Url = "http://localhost/ArcGIS/services/NetworkAnalysisService/MapServer/NAServer");
            //Get Route layers
            string[] routeLayers = naService.GetNALayerNames(esriNAServerLayerType.esriNAServerRouteLayer);

            // Get the default settings for the first route layer
            NAServerRouteParams naServerRouteParams = (NAServerRouteParams)naService.GetSolverParameters(routeLayers[0]);

            // Create 2 stops as property sets
            List<PropertySet> propertySets = new List<PropertySet>(2);
            propertySets.Add(CreatePropertySet("My starting point", -117.195905, 34.057783));
            propertySets.Add(CreatePropertySet("My ending point", -117.180943, 34.056286));

            // Set the property sets on a NAServerPropertySets object and set on the routeParams object
            NAServerPropertySets naServerPropertySets = new NAServerPropertySets();
            naServerPropertySets.PropertySets = propertySets.ToArray();
            naServerRouteParams.Stops = naServerPropertySets;

            // specify to return compact directions
            naServerRouteParams.ReturnCompactDirections = true;

            // Solve
            NAServerRouteResults naServerRouteResults = (NAServerRouteResults)naServer.Solve(naServerRouteParams);

            // Output directions
            NACompactStreetDirections streetDirections = naServerRouteResults.CompactDirections[0];
            StringBuilder sb = new StringBuilder();
            sb.AppendFormat("DriveTime = {0}", streetDirections.Summary.TotalDriveTime);
            sb.AppendLine();
            foreach (NACompactStreetDirection streetDirection in streetDirections.Directions)
            {
                sb.AppendLine(streetDirection.Text);
            }
            Console.WriteLine(sb.ToString());
        }

    }
    }
}


texts: Moving annotations, or texts, without making changes in enterprise geodatabase?

texts:
I'm working in a company, producing cadastral maps in ArcGIS where some clients prefer the maps being delivered in a .dwg-format. We're using a tool which preserves the feature class styles when importing them into AutoCAD.

The problem is how to handle som texts in the cartographic visualization. As we're not really allowed to move around objects in the enterprise database but still want some adjacent annotation classes, eg. administrative property information, to be visualized in the layout and the product, I can't really find a solution.

Converting the text to graphic would label the layer 'Graphic' in the DWG.
Exporting the layer to a file geodatabase and make the edits solves the problem. But might be time consuming.
Entering edit mode, move the annotations, make the export and stop editing without saving seems to be the fastest practice but at the same time a bit risky if the user by accident saves the edits.

So what I basically want is a method of making changes in the text, that doesn't store them in the enterprise database but still preserves the attributes of an annotation class.

Does anyone have a solution for the issue?


texts: Selecting several layers programmatically using PyQGIS?

texts:
In my plugin I want to select some layers programmatically. But with

virtual QList<QgsMapLayer *> QgsLegendInterface::selectedLayers(bool inDrawOrder = false)

I only receive the actually selected layers.

I've tried in QgsLayerTreeView the method void QgsLayerTreeView::setCurrentLayer(QgsMapLayer * layer) But I only can select one Layer, and I need to select more layers.

Any idea what could I try?


texts: Determining Appropriate DBDB-V Subdataset

texts:
I am using DBDB-V v6.2 and GDAL, specifically the dbdbv6_level0c.h5 HDF5 file and the GDAL 2.3.3 C# binaries.

Using these lines of code I can easily find the full list of Subdatasets used in the raster:

var path = @"F:\test\dbdbv6_level0c.h5";
var raster = Gdal.Open(path, Access.GA_ReadOnly);
var metadata = raster.GetMetadata("SUBDATASETS");

This produces an array of strings, pairs of which contain subdataset names and descriptions, for example:

SUBDATASET_150_NAME=HDF5:"F:\test\dbdbv6_level0c.h5"://geo/0.05000/229_543/depth
SUBDATASET_150_DESC=[600x600] //geo/0.05000/229_543/depth (32-bit floating-point)

This then allows me to apply Gdal.Open(...) to any of Subdataset NAMEs and open up that particular raster for processing.

It can also be correlated with the basic metadata, which has information such as:

geo_0.05000_229_543_depth_compile_date=20091230
geo_0.05000_229_543_depth_datatype=depth
geo_0.05000_229_543_depth_ellipsoid_offset=0
geo_0.05000_229_543_depth_units_of_measure=meters
geo_0.05000_229_543_depth_vertical_datum=MSL
geo_0.05000_229_543_depth_zoo_value=0
geo_0.05000_229_543_east=272
geo_0.05000_229_543_north=25
geo_0.05000_229_543_security_key=100
geo_0.05000_229_543_south=24.5
geo_0.05000_229_543_west=271.5

Note that 229_543 is the shared identifier between the Subdataset Name and the metadata, which provides the geographic coordinates of interest.

My question is, can I start with known geo coordinates (e.g., 24.75N 271.75E) and somehow determine that the correct Subdataset to open is the one identified with 229_543?

Please note, I am able to open the Subdataset just fine, but what I want is to be able to determine which Subdataset to open based on geo coordinates of interest, without having to iterate the entire metadata collection for every new coordinate of interest... ideally, some sort of algorithm that takes in geo coordinates and spits out DBDB-V Subdataset identifiers.


texts: How to check whether a feature has been split in GeoTools?

texts:
In my program I am looking for a best place for a golf course. Three of the conditions are:

The course must be on a hill,
it must be at least 100 meters away from roads and
must not be larger than x.

Some of the generated courses are divided into two pieces by a road in the middle.

My question is:

How do I get those pieces as separate features or geometries to check their area? Is that even possible? Can I  check if a certain geometry is whole, not split?

Edit: I am using geotools platform. (Feeling stupid that I forgot to mention)


texts: Solutions for generating grade/slope estimate heatmaps from a topography contour layer

texts:
Given a shapefile layer of 5-ft interval topography contour lines overlaid on google earth, I could, using the Ruler tool and some simple arithmetic, calculate average percentage grade numbers for a sample of points in an area.

What are some ways I could do calculations that like using a GIS program to ultimately generate a heatmap of grade estimates for an area? In terms of functionality, I suppose it would segment out the spaces between contour lines by some predefined length (or other method based on variation...), calculate the average grade between those lines, and output a number, which would then be used to color code each given polygon.

I have a fair amount of GIS scripting experience (ruby/rgeo/geos) and could probably cook up a custom solution for this, but given my experience in the GIS space by now, I have come to learn that it is best to see what's out there already in terms of pre-packaged solutions and/or established algorithms (for example, Hexagonal tessellations rather than cubes) before breaking out the custom scripting.


texts: Disappearing labels on large maps' edge - Easy Custom Labeling Plugin - QGIS 2.8.3

texts:
I've been using the Easy Custom Labeling plugin to fine-tune placement of labels. However, I have been having some labels disappear when they are on the edge of large maps. For example, I have an A0 sized map where labels on the very right side disappear. Once I pan the map in the print composer towards the left, these labels appear little by little.

Anyone else experience this? Any way to fix it? For now, I am working around it by adding a second, smaller map on top of the area that has missing labels to sort of "mask" it.


texts: openlayers3 Draw interaction vertexadd event

texts:
I am trying to implement an undo/redo functionality during drawing a fetaure using the draw interaction. For the modify interaction I was able to to do it using the 'modifyend' event. like so:

modifyInteraction.on('modifyend',function(e){
}

while for the draw interaction no such events exist.
'drawend' event fires when finishing my drawing. But I need to trace all the vertices drawn in the meantime.
For the time being I came up with this:

myapp.Editor.editDraw = new ol.interaction.Draw({
    source: myapp.Editor.PseudoEditLyr.getSource(),
    type: type,/** @type {ol.geom.GeometryType} */
    style: drawStyle,
    condition : function(e){
        //fires the mousedown event while drawing
        console.log("e.browserEvent.type",e.browserEvent.type);
        if (myapp.Editor.editSketch ){
        myapp.Editor.checkAndPushUndo(myapp.Editor.editSketch);
        }
        return true;//always return true but need to listen for the mousedown event
    }
    });

Which does the job but this is not elegant, as it trace the 'mousedown' event and not the 'vertexadd'.

Do any of you have any idea?


texts: How to overlap a base map with a png map (format .pgn) on QGIS-Pisa?

texts: How can I overlap a relief base map taken from maps-for-free.com  and this png-formated map  so that I can get a resulting map with relief of the 1st map and boundaries of the 2nd map ?? I'm a beginner. I don't know if QGIS has any functionality like that ? I use QGIS version 2.10.1-Pisa and my operating system is Windows 7 x64bit.


texts: Importing ArcPy taking many seconds and code is really slow for GeoProcessing Tools 64-bit (ArcGIS 10.1 SP1)

texts:
I've installed ArcGIS10.1 SP1 and whenever I import arcpy (on Pywin console or Eclipse PyDev) it takes aprox. 15-20seconds. Is this normal? I don't think so...

I had Python 2.7.3 64-bit previously installed on my Windows 7 (and still use it) and now I also have Python 32-bit and Python 64-bit that came with ArcGIS.
I think I might have the paths confused, because the 64bit version of Python 2.7 in the folder C:\Python27\ArcGISx6410.1was supposed to be 2.7.2, but when I check the version, inside that path, it gives me:

C:\Python27\ArcGISx6410.1>python.exe
Python 2.7.3 (default, Apr 10 2012, 23:24:47) [MSC v.1500 64 bit (AMD64)] on win
32

My %PATH% has this, among other paths:

...
C:\Python27
C:\Python27\Scripts
C:\Program Files (x86)\ArcGIS\Desktop10.1\Bin
C:\Program Files (x86)\ArcGIS\EsriProductionMapping\Desktop10.1\Bin
...

And I didn't need to set PYTHONPATH. My ArcMap Python shell works fine and so does the 64-bit version for the backgorund scripting.

My code is very simple: basically I have to find how many times does a node (NO) from a node feature class appears in other line feature class. It works, but for #6290 nodes, it takes more than 10min!

'''

# Import system modules
import arcpy, sys
from arcpy import env

print "Running against: {}".format(sys.version)

def do_search_by_cursor(table,searchString):
    count = 0
    count += cursor_search_count(table,searchString)
    #print "done by cursor: " + str(count)
    return count

def cursor_search_count(table,searchstring):
    rows = arcpy.SearchCursor(table,searchstring)
    count = 0
    for row in rows:
        count += 1
    return count

# Set the workspace environment to local file geodatabase
env.workspace = "S:/workspace/MSc_git/ArcGis/estagio_db.gdb"

# Set local variables
fc_nodes = "Rede_Cbr_MM_node_proj"
fc_links = "Rede_Cbr_MM_2008_WithoutMetro"
# fieldName1 = "LEGS"
# fieldAlias = "numLegs"

fieldname="TONODENO"
# Create field name with the proper delimiters
delimitedfield = arcpy.AddFieldDelimiters(fc_links, fieldname)

# Create search cursor:
rows = arcpy.SearchCursor(fc_nodes, "", "", fields="NO")

for row in rows:
    searchStr= delimitedfield +" = " +str(row.NO)
    myvar= do_search_by_cursor(fc_links, searchStr)
    print "Node %d has %d legs" %(row.NO,myvar)

My PyDev is configured to use the specific Python Interpreter that came with ArcGIS processing Tools 64-bit and looks like this (I've checked all...):

If you could give a hint or a link (I've already explored, but nothing seems to work), I'll appreciate it.


texts: What's the difference between Geographic coordinate system 4019 and 4121?

texts:
I have a DEM with cell size 100,100 with geographic coordinate system GRS_1980_Transverse_Mercator, Datum: D_GRS_1980. I created some extra shapefiles with projected coordinate system: Greek_grid, projection: Transverse_Mercator and geographic coordinate system: GCS_GGRS_1987 and datum:D_GGRS_1987.
Adding them in ArcMap I have the message that my layers use different spatial reference which is reasonable because the shapefile datum was different from the DEM I am using.
I checked online and found in http://edndoc.esri.com/arcims/9.2/elements/gcs.htm the geographic coordinate systems I am using. They are the 4121 and 4019. They look to have the same characteristics. Does anyone know why they exist with different names? Or am I missing something?

4121    GCS_GGRS_1987
GEOGCS["GCS_GGRS_1987",DATUM["D_GGRS_1987",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]

4019    GCS_GRS_1980
GEOGCS["GCS_GRS_1980",DATUM["D_GRS_1980",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]

Thanks in advance,
chch


texts: Assign Z value to start & end of PostGIS line

texts:
I have a line table in my PostGIS database and a table with invert levels for the start and end of each line.

My goal is to create create 3D lines in my table, from my invert data.  I know I can use st_makeline to recreate the line with a 3D point at the start and end but I was hoping to be able to run a update just on the Z value for the start and end.

Can this be done or is it easier to just recreate the line.


texts: gdal_translate Resample different min max values

texts:
I tried to resample a Sentinel2 GeoTiff with a 20m resolution to 10m resolution using gdal_translate:

result = gdal.Translate(outFile,resampleBandPath,format="GTiff",width=10980,height=10980,resampleAlg="nearest",outputType=gdal.GDT_UInt16)

I would like to have one raster value in the 20m repeated four times in the 10m resolutuion output Tiff.

Computing statistics on both images:

resample.GetRasterBand(1).GetStatistics(0,1))
old.GetRasterBand(1).GetStatistics(0,1))

gives me:

[27.0, 8337.0, 799.74556483486, 763.06787190687]
[29.0, 9257.0, 801.11717304189, 768.93310461989]

Why is the output image "smoothed"?


texts: Adding PDF (especially GeoPDF) to ArcMap as data?

texts:
I receive a lot of PDF maps (vector, image, with georeference and without) and I am having to convert them to an image file before I can add them to ArcMap and georeference.

I can't see any obvious out-of-the-box utility that would allow me to add and georeference a PDF directly without having to convert to an image file first.

Any ideas?

I would prefer to add a PDF to map directly without any prior "pdf to image" conversion


texts: Spatial interpolation of 30 weather stations to other areas

texts:
I have monthly (and sometimes daily) data on temperatures from about 30 Swedish weather stations in the mid-19th century that I want interpolate to the whole country. I'm using ArcGIS and Stata.

For later years, I have data on additional weather stations (around 100), so I've thought of using the later data as to calibrate a model using the old stations to predict weather at the new stations. Since I know the actual weather at the new stations, I could calibrate to model to achieve the best possible fit. But I'm not sure what a good method for going about achieving a suitable fit (don't want to risk overfitting for example).


texts: Constructing a Voronoi diagram using a complicated travel time metric

texts:
I have in my possession four sets of shapefiles:
1. Locations (points), around 1000
2. Roads (lines), each road has a speed attribute
3. Terrain (polygons), each zone also has a speed attribute
4. Rivers (lines), each river has a crossing delay time attribute

How can I construct a Voronoi diagram for my location using travel time as a distance metric? I could have used pgRouting if I had roads only, but I am not sure how to approach different off-road terrain and rivers. I'm willing to code this, but I've spent enough time trying to come up with an algorithm, but to no avail.

My current plan is to dump it all into a sufficiently large raster and create 1000 distance maps out of it, then use them to calculate the Voronoi cells, but I have a suspicion that a sufficiently large raster will take too long to process. Is there a better vector-based algorithm?


texts: Loading point features into Network Analyst without losing the point's attribute fields

texts:
I have points that represent schools (facilities) and points that represent homes (incidents). I need to do a closest facility analysis to calcualte how long the residents need to reach the schools. However, in order to calculate an average traveltime per block of houses, I need the traveltime to be linked to the homes and not to the routes.

The problem: ArcGIS deletes every attribute field of the point layers (including the ID Field), when loading them into the Network Analyst. That makes linking the traveltime from the routes onto the homes very tedious (e.g. via joining the routes to the incidents using their ID and then spatially joining the incidents with the homes).

Is there a way to keep the attribute fields of the original point layer when loading them into the network analyst?

The same problem appears when loading points into a service layer analysis.


texts: Shapefile disappeared after changing CRS

texts:
I change the coordinate system (CRS) of a point data but it disappeared. I tried zooming to layer but its not visible.
I am using QGIS 2.14


texts:
I have a Garmin 60csx, about 2 years old and I have attempted to connect it to my laptop/ArcView (GPS toolbar) using the Detect GPS Port to no avail.  There is no port number being displayed and that is what it states in the popup box.

the baud rate, parity etc... are accurate.  Should the Comm Port have a number displayed upon connection if it is working?

Is there a patch or a specific usb cable to use or specific Garmins will only function properly in ArcVIEW?

now it is saying the USB is not recognized.  I'm thinking a new USB cable as one solution.


texts:
For debugging it would be great to figure out which EE object type an object is.

Example code which adds a buffer to a point in a random area:

var polygon = ee.Geometry.Point([0,51.477]).buffer(10000); // a circle in London

var random_points = ee.FeatureCollection.randomPoints(polygon.bounds(), 100); // 100 random points in London
var point_list = random_points.geometry().geometries();
print(point_list); // prints -- List (100 elements)
var one_point = point_list.get(10); // The 11th point.
print(one_point); // prints -- Point (-0.03, 51,52)
var small_circle = one_point.buffer(100);
print(small_circle); // Error -- one_point.buffer is not a function...

Since the console prints Point..., I would think that one_point is of type ee.Geometry.Point. However, according to the docs ee.Geometry objects are supposed to support the .buffer() function, which this one does not support. Hence I assume it's a different kind of object but I have no clue which. When I code in Earth Engine, I run into these kind of problems very frequently.

Is there a way to figure out which type of EE object an object is? Something like the fantasy .EE_OBJECT_TYPE() function here:

print(polygon.EE_OBJECT_TYPE()); // should print -- ee.Geometry.Polygon
print(one_point.EE_OBJECT_TYPE()); // should print -- ee.WhateverObjectThisIs

Non-solutions:

print(my_object) often indicates the type of object but, as the example above demonstrates, this is not always the case (or correct).
my_object.type() only specifies the kind of ee.Geometry and doesn't work on any other EE objects according to the docs.
my_object.typeOf() always returns object for EE objects.


texts: Can't add layers to new PostgreSQL database

texts:
I made a new PostgreSQL database and am trying to add layers/tables to it in QGIS using the DBManager. And I am failing.

Here is the error message I'm getting:

However, I can add layers to other PostgreSQL databases that previously existed. I thought it might be some problem in the way I created the new database, using pgAdmin 4, but as far as I can tell the new database is exactly the same as the others.

I tried some workarounds, including creating a new schema in an existing database and couldn't read from that new schema.

Any suggestions?


texts: Flattening multi-band raster to single-band raster in ArcMap?

texts:
I need to convert a multi-band PNG image into a single-band raster, for analysis and visualisation.

How can I do this in ArcGIS Desktop 10.2?


texts: Sort feature class in Modelbuilder and populate new ids

texts:
I would like to sort my feature class in a model based on coordinates stored in its attribute table. The sorting must be done on Y_coord descending, then X_coord ascending. The answer to Sorting feature class to calculate sequential ID field using ArcGIS Field Calculator? looks like what I need, except I only have an ArcEditor Licence, and free versions of Xtools Pro and ET Geowizards. My model performs several other processes as well, so I need the sorting to be included in it. How can I perform the sort with the above limitations?

However, looking at the options in the Sort tool, I've realised that what I actually want to do is sort my features spatially, starting from the upper left, but to sort based on the Shape field requires an ArcInfo licence. The sort tool from ET Geowizards does not allow sorting on this field either. Are there other tools which could do this?


texts: How Can I Merge Geocoded Addresses With a Shapefile and Export One Table?

texts:
I currently have two layers in my QGIS (v1.7.3 for Mac) project -- one being a Pennsylvania State House shapefile and the other being 166 geocoded addresses (for which I used 'Add Text Delimited Layer' with XY latitude and longitude columns as coordinates). The mapping of these two worked great:

But is there a way for me to get the data out of QGIS and back into an Excel file, with regards to which state house district each geocoded address falls into?

Ideally, the data would come out similar to:

    Latitude Longitude District
    40.14524 -75.21491 131

Thanks!


texts: Attributing stream segments based upon road intersections and stream confluences

texts:
I am currently working on attributing nodes along some stream polylines. The attribution is supposed to include:
1)length of the stream from it's reference point
2) the location of the node in relation to its reference point - upstream/downstream)
3) either a road crossing or a stream confluence, whichever is closest. This will be the reference point.

I am trying to find a way to automate this process (if possible) because I need to attribute somewhere around 2,000 points.

My question is, how would you go about doing this?

Here is what I have conceptually worked out so far:

I need to attribute each stream node to the closest either road or confluence point. For this I would do a Spatial Join and set the Match Option to "Closest"
Once the segments are attributed then the Locate Features Along Routes tool can be applied to gather the length of the segment to its initial point/vertex (determined by the initial attribution)
The direction (upstream/downstream) can be determined by the directionality of the way the line was drawn

The main problem I see with this process is that I will be doing a Spatial Join to the node and it will be joining the closest node and not to the confluence point on the river or closest road, which will not give me the correct length once I run the Locate Features Along Routes tool.

Here is why I don't think my concept is working yet:
When looking at Rock Creek Tributary 12, the node at 9th Street should be referencing 9th Street, since it is closer. But the node right by the confluence should be referencing the confluence point with Rock Creek. I am not entirely sure how to be sure that the nodes will be attributed to the closest either confluence or road.


texts: How to calculate Spatial Autocorrelation of grouped points in ArcMap 10.2.2

texts:
I have 32 points uniformly distributed with 7 attributes that comprises continuous values. I used the Grouping Analysis tool in ArcMap 10.2.2. I set 2 groups that have clear differences and are also proposed as best number of groups by the tool.
I also didn't use spatial constraints on purpose.

My goal is to show that the pattern of those grouped points isn't random. In other words, I'd like to apply something like Incremental Spatial Autocorrelation (global Moran's I) tool on the grouped points... I searched in here, but didn't come up with something definite.

Notably, I applied global Moran's I on the grouped points using the "SS GROUP" attribute as an Input Field, both with and withoute row standardization. The results showed no errors, as well as positive but weak autocorrelation.


texts:
I have assigned one method to combobox like this

def run(self):

    """Run method that performs all the real work"""
    GetAllLayers(self)
    #attach index changed event / passing parametrized method
    self.dlg.cbLayerNamesAll.currentIndexChanged.connect(lambda arg=self:  LayersValueChange(arg))

getting error here

def LayersValueChange(self):
    print "abc"
    #selectedLayerName = str(self.dlg.cbLayerNamesAll.currentText())
    #selectedLayerName = self.dlg.cbLayerNamesAll.currentText()

    layers = self.iface.legendInterface().layers()//here error

And Error is:

  layers = self.iface.legendInterface().layers()
  AttributeError:
  'int' object has no attribute 'iface'

self is object but it getting it like int.


texts: ModelBuilder - How to calculate Field using part of path name

texts: I am trying to calculate a field using part of the feature class path name (specifically, the last three characters of the geodatabase name). Is this possible? I've tried using the parse path tool and including the value as an inline variable subsitution with the VB string right function. Perhaps a more complicated Python script is neccessary?


texts: Finding GIS consultancy/project work (rather than jobs)?

texts:
The post What are useful websites for scouting GIS related jobs? listed resources for finding a GIS job (as in a full-time or contract position).

What resources are available for finding GIS projects or consultancy work? eg:

a small project to integrate a map into an existing site;
incorporate an address-search into a workflow;
migrate legacy GIS code to a new format;

These are one-off tasks which require a GIS professional, but which don't entail a full-time position. How do you find the right people for the task, or find the task if you're the right person? All the online resources I've seen are aimed at "jobs".


texts: Generating route from user's current location to desired location on web map created in QGIS?

texts:
I am developing a GIS in QGIS and I want to get the user's current location, then according to the user's second point selection, I want to generate a route between these two points. If possible, I want routing to be based on walking, not car.

I want to use OpenStreetMap as the base map.

Any idea how I can achieve this?

I want this to be done on a web map after creating it using the qgis2web plugin. I have looked at OSRM but I can't find out how exactly to use it on this system.


texts: Preserving line segment ObjectIDs in point output from line Intersect?

texts:
I am using the Intersect tool in ArcGIS for Desktop 10.1 to identify all the intersections in a polyline shapefile.  The tool generates a point shapefile, but when I inspect the data it does not preserve the ObjectIDs from the original polyline file.

How can I preserve the two line segment ObjectIDs that create the intersection?


texts: Convert Microstation dgn v8i files to shapefiles in Linux

texts:
I am looking for some program which would allow me to convert Microstation .dgn v8i files to OGR supported format (preferably ESRI shapefile) in Linux.

I have tried GRASS and QGIS, but both support only dgn version 7 or lower.


texts:
I have some commuting data showing people who cycle to work. The data has residence, workplace and unique ID numbers for each commuter. I have georeferenced the home and work locations so I have two point shapefiles, work and home. I also have a shapefile of the road and cycle route network.

I have done some experimenting and the shortest path algorithm in QGIS 3 allows shortest route analysis of shapefile to point OR point to shapefile (which is super), but I need to calculate a route for EVERY commuter. Is there a way to automate this using the ID number as the base?

From looking into this online it seems that this may be possible using PGRouting in PostGIS, but it would be nice to find a solution avoiding PostGIS because my colleagues do not have access to this.

This stackexchange question alludes to a similar problem, but, again, this user is only using one point. I have not found a stackexchange thread which answers this problem.
QGIS 3 - Shortest path layer to layer, not point!


texts: Reshaping polygon to match other polygon/polyline without using Reshape/Trace or Align to Shape in QGIS

texts:
I am trying to edit the geometry of one polygon to match another polygon, same as in this question:
Snapping polygons while reshaping edges?
Unfortunately, using the Reshape and Trace editing tools is achingly slow and tedious and I have a bunch of locations spread across a bunch of polygons that require this editing.  I also tried the Enable Tracing tool in QGIS as this thread suggests: Adjust polygons by edition (freehand and reshape tools used together) to no avail.
Are there any other approaches that have worked for others to do this efficiently?
Just doing a rough reshape and then clipping or erasing (as suggested in the first thread) gets complicated since the lines intersect in both directions (see image).  I want to align the low-res black polygon to the high-res blue polygon in specific places.