texts: Single Positional Indexer is out of bounds

texts:
I am playing with OSMNx and I would like to extract street blocks,

I have this code which should display the blocks

import osmnx as ox
import geopandas as gpd
import shapely
import matplotlib.pyplot as plt

print("Start")
point = (49.224251, 16.578077)

streets_graph = ox.graph_from_point(point, distance=500, network_type='drive')
streets_graph = ox.project_graph(streets_graph)
print("Finished graph")
streets = ox.save_load.graph_to_gdfs(streets_graph, nodes=False, edges=True,
                                     node_geometry=False, fill_edge_geometry=True)

point = streets.unary_union.centroid
print("Found Point")
polygons = shapely.ops.polygonize(streets.geometry)
polygons = gpd.GeoSeries(polygons)
print("Made Polygons")
target = polygons.loc[polygons.contains(point)]

target_streets = streets.loc[streets.intersection(target.iloc[0]).type == 'MultiLineString']
print("Beginning of plotting")

ax = polygons.plot()
gpd.GeoSeries([point]).plot(ax=ax, color='r')

plt.show()

print("Done")

But when I am trying to define target_streets I got an IndexError that single positional indexer is out-of-bounds, The error depends on a given GPS point, it works with some points but most of the time it does not work.
for example when I use this point:

point = (40.742623, -73.977857)

it does exactly what I want it to do

How can I fix this error?
Thank you very much


texts: Extracting German restaurants from planet.osm?

texts:
I want to get the data out of http://wiki.openstreetmap.org/wiki/Planet.osm

I have just come back from the site http://wiki.openstreetmap.org/wiki/Planet.osm

which has got some files that belong to the openstreetmap-Project. i downloaded a similar file that is derived from a German Mirrorsite - this one here:

http://ftp5.gwdg.de/pub/misc/openstreetmap/download.geofabrik.de/

i downloaded the file called  germany.osm,bz2

according  to this instruction

  type bzip2 -d planet.osm.bz2; or your OS may support double-click
  unpacking. See Wikipedia's list of compression programs.

i did the following - but without success-,

martin@linux-wyee:~/poi_abzug/geofabrik> bzip2 -d germany.osm,bz2;
bzip2: Can't open input file germany.osm,bz2: No such file or directory.

Can anybody explain why i did not succeed?

further questions:

Is this a XML-File ?

I only need the data for the Pois in Germany: can i use this file here - http://ftp5.gwdg.de/pub/misc/openstreetmap/

As this might be of general interest i think that many users can benefit from an answer to this question.

So my question is:

How to read the above mentioned dataset - how to read it - on a OpenSuse-Linux-machine!?
Should i do it in  the terminal?

Note - i want to store all in a csv file - and i look for all the restaurants... What is the easiest way to do this?


texts: Surface and Volume from TIN read through text file?

texts:
I am trying to create a script that will iterate through a list of elevations in a text file and calculate surface and volume using the 3D analyst surface volume tool. I have success for the first line of a file only. It seems to stop after it's calculated the first value and I can't for the life of me figure out why it will not continue. Looking at the second elif statement with the height_file.

#Purpose: To calculate area and volume on a surface.

# Create the Geoprocessor object
import arcpy
from arcpy import env
import exceptions, traceback, os, sys

#Check out the 3D Analyst extension
arcpy.CheckOutExtension("3D")

try:
    # Set the workspace (to avoid having to type in the full path to the data every time)
    # Select the 3D Analyst Toolbox
    env.toolbox = "3D"
    tin = sys.argv[1]
    output = sys.argv[2]
    reference_plane = sys.argv[3]
    z_factor = sys.argv[4]
    zmin = int(sys.argv[5])
    zmax = int(sys.argv[6])
    base_z = sys.argv[7]             # For specified height value, conflicts with height text file and height interval, inactive value=-9999
    interval = float(sys.argv[8])    # For specifying an interval, conflicts with specified height value and height text file input, inactive when zmin&zmax=0
    height_file = sys.argv[9]       # For height text file input, conflicts with specified height value and height interval inactive when no input
    print "Base Z = " + str(base_z)

    if base_z != '-9999':
        # Process: Calculating areas and volumes on a TIN above a reference plane
        print "Calculating single base value: " + str(base_z)
        arcpy.SurfaceVolume_3d(tin, output, reference_plane, base_z, z_factor)
    elif zmin != 0 and zmax != 0 and interval != 0 and zmax > zmin:
        print "Processing intervals"
        elev = float(zmin)
        while elev < zmax + interval:
            print "Height: " + str(elev)
            arcpy.SurfaceVolume_3d(tin, output, reference_plane, elev, z_factor)
            elev = elev + interval
    elif height_file != "":
        print "Processing elevation file"
        fo = open(height_file, 'r')
        while True:
            elevInString = fo.readline()
            if not elevInString:
                break
            a = elevInString.strip()
            if a == "":
                continue
            elev = float(a)
            print "Height: " + str(elev)
            arcpy.SurfaceVolume_3d(tin, output, reference_plane, elev, z_factor)
        fo.close()

except arcpy.ExecuteError:
    print arcpy.GetMessages()
except:
    # Get the traceback object
    tb = sys.exc_info()[2]
    tbinfo = traceback.format_tb(tb)[0]
    # Concatenate error information into message string
    pymsg = 'PYTHON ERRORS:\nTraceback info:\n{0}\nError Info:\n{1}'\
          .format(tbinfo, str(sys.exc_info()[1]))
    msgs = 'ArcPy ERRORS:\n {0}\n'.format(arcpy.GetMessages(2))
    # Return python error messages for script tool or Python Window
    arcpy.AddError(pymsg)
    arcpy.AddError(msgs)


texts: Drawing polygon using parts of existing polyline in ArcGIS editing?

texts:
I use AcrGIS 10.2, and I met a problem about drawing a polygon "effectively".

There is an existing line (see above). I want to draw a polygon whose part of the border is this polyline. I can catch every vertex of the line to make the polygopn match with the polyline. But this method may become impossible (and tiring) when this existing line has many turning points.

Is there some useful method to do this?


texts: How to Enable/Disable Topological Editing Programatically

texts:
Is there a way to enable or disable topological editing with python in QGIS?
I'm trying to set a bunch of rules, one of them being the editing, when starting up a project, but I can't seem to find how to enable this.
For that matter, there's even no
iface.SnapppingToolBar()


texts: Error regarding missing feature property when training/classify image in Earth Engine

texts:
I am trying to run a land cover classification in GEE. I have a csv file which I converted and exported as a shapefile in GEE. The problem is when I specify the column name containing the land use class, I get the error below.

  classification: Layer error:
  Property 'land_use_5' of feature '00000000000000000011_0' is missing.

What am I doing wrong? Here is my code:

https://code.earthengine.google.com/15a705dc1133e6978f09e5113bbbcb99

var rwanda = ee.FeatureCollection("users/sbrightaboh/edem");
print(rwanda)
//var count = rwanda.aggregate_count("_land_use_category_alias_class")
//print(count)
//var lookupIn = ee.List([0,1,2,3,4,5,6])
//var lookupOut = ee.List([6,5,4,3,2,1,0])
//var remaped   = rwanda.remap(lookupIn, lookupOut, "land_use_c,C,254")
//print(remaped)
Map.addLayer(rwanda,{color:"red"},'Rwanda');

//var values = ee.Dictionary

// Input imagery is a cloud-free Landsat 8 composite.
var l8 = ee.ImageCollection('LANDSAT/LC08/C01/T1');

var image = ee.Algorithms.Landsat.simpleComposite({
  collection: l8.filterDate('2018-01-01', '2018-12-31'),
  asFloat: true
});

// Use these bands for prediction.
var bands = ['B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B10', 'B11'];

// Load training points. The numeric property 'class' stores known labels.

// This property of the table stores the land cover labels.
var label = 'land_use_1';

// Overlay the points on the imagery to get training.
var training = image.select(bands).sampleRegions({
  collection: rwanda,
  properties: [label],
  scale: 30
});

// Train a CART classifier with default parameters.
var trained = ee.Classifier.cart().train(training, label, bands);

// Classify the image with the same bands used for training.
var classified = image.select(bands).classify(trained);

// Display the inputs and the results.
Map.centerObject(rwanda, 11);
Map.addLayer(image, {bands: ['B4', 'B3', 'B2'], max: 0.4}, 'image');
Map.addLayer(classified,
             {min: 0, max: 2, palette: ['red', 'green', 'blue']},
             'classification');


texts: Can polygons marked by more than one values reflect the average of them? (Carto)

texts:
The color and pop-up info of my polygons (zip codes) depend on a variable that takes many values for each polygon (One for each combination of few categories and dates).

I can see that the range of colors (and histograms with autostyle) takes all the values of the variable. The average is well done, but the color and the pop-up value of the polygons is not the average of all the values from each polygon. In the image with two values it takes the minimum value in both the color and the pop-up, but not always is neither the minimum value the chosen one.

I cannot make a new dataset with the mean values because they change depending on the filtering with other widgets

Is it possible that the color and the pop-up of the polygons reflect the average values of the variable they are showing?


texts: Boolean field (attribute type) in PostGIS with OGR Python

texts:
In my project I'm using ogr (python bindings) for inserting spatial data in PostGIS.

The spatial data comes in all different formats (CSV, Shapefile, GeoJSON) etc and before creating new tables in PstGIS I'm reading the properties (or attributes) of the features and do some processing.

I've ran into an issue when I want to create a boolean field in the PostGis table. I could not find a boolean type in the OGR variable list that can be used for creating a field definition as I do with string and integer types:

    # For string type
        ogr.FieldDefn(field_name, ogr.OFTString)

    # For integer type
        ogr.FieldDefn(field_name, ogr.OFTInteger)

Does anyone know how to create the Postgresql boolean type with OGR Python?
Have I missed something in the OGR Python docs ?


texts:
i was create a simple export code using pyqgis3 and atlas print composer map layout on QGIS3.

my question is about map scale.

if i use this line :

mapScale = map.scale()

i take the scale from map and i can use this number on grid interval as i want.

but after generate atlas mode and using this settings for atlas scale

layout.referenceMap().setAtlasScalingMode(QgsLayoutItemMap.Auto)

then the map scale change.

i want to take with some way to take the scale from the atlas features after generate or how to use expression on define IntervalX like this .

expression123 = "$scale*10"
qgs_layout_item_map_grid = QgsLayoutItemMapGrid("map 1", map)
#qgs_layout_item_map_grid.intervalY(50)
#qgs_layout_item_map_grid.setIntervalX(50)

i try to use this but without success.

qgs_layout_item_map_grid.dataDefinedProperties().property(QgsLayoutItemMapGrid.intervalY).setExpressionString(expression


texts:
I'm trying to reproject a LAZ file from EPSG 2226 to EPSG 6418.
The reprojection works fine (I followed this: Can I re-project an LAS file in LidR) but I cannot store the information onto the LAS header and subsequently, as well as I save my LAS file (and get out of R environment), the CRS is lost.
Download this LAZ file
I used the code below:
library(lidR)
las_tmp <- readLAS("./USGS_LPC_CA_Sonoma_2013_SOCO_0054_108_LAS_2017.laz")
las_tmp_reproj <- spTransform(las_tmp, sp::CRS("+init=epsg:6418"))

#> Warning messages:
#> 1: In showSRID(uprojargs, format = "PROJ", multiline = "NO") :
#>   Discarded datum NAD83_National_Spatial_Reference_System_2011 in CRS definition
#> 2: In sp::proj4string(x) : CRS object has comment, which is lost in output
#> 3: EPSG code not found: header not updated. Try to use the function epsg() manually to ensure CRS will be written in file.
#> 4: In showSRID(uprojargs, format = "PROJ", multiline = "NO") :
#>   Discarded datum Unknown based on GRS80 ellipsoid in CRS definition

As said by warnings() number 3, the epsg code is not found and header not updated. I followed the tips and tried the following command:
epsg(las_tmp_reproj) <- as.integer(6418)

#> Warning message:
#> In showSRID(SRS_string, format = "PROJ", multiline = "NO") :
#>   Discarded datum NAD83_National_Spatial_Reference_System_2011 in CRS definition

Unfortunately, after this step, the access to the epsg information epsg(las_tmp_reproj) throws me the following error:
epsg(ctg_tmp_reproj)
#> Error in if (tags[[i]]$key == 3072) return(i) :
#>   l'argument est de longueur nulle

Is it related with the datum itself? The error message is not really explicit and I'm pretty stuck


texts: How can I combine multiple raster datasets with multiple unique values in ArcMap 10.6?

texts:
I have to combine plenty of raster datasets. The problem is there are to many unique values. I guess that’s because there are layers where areas have geometry and there are not geometry on the other layers on the same areas. So, I have to combine raster layers with a little bit different geometry. I tried Combine function but it did not work because of the unique values. The resolution is the same on all layers.

The attribute table should look like something like this:

Here is an example for two layers:


texts: Reprojection not generating tiles properly

texts:
I am executing this query with cartodb.js. I am moving the markers randomly to generate a new layer.

  cartodb.createLayer(map, {
      user_name: 'mickemr',
      type: 'cartodb',
      sublayers: [{
        sql: "SELECT  cartodb_id, user_email,  apo_about_me,apo_headshot_url,apo_member_number, first_name, user_nicename,  ST_Transform(ST_Project(the_geom, 300*random(), random()*(3*pi()/2))::geometry , 3857) AS the_geom_webmercator   FROM users",
        cartocss: "#users  {marker-fill: red;marker-width:20px; } "}]
    })

My problem is that the tiles generated are not created properly, they are cut in the middle:
Tile: http://screencast.com/t/5WFESgZkrG3
Chunk of map: http://screencast.com/t/g7RVmIujvye

If I don't use the random function it is working properly.

Anybody know what could be happening?

Is the SQL query wrong?


texts: Need less significant digits returned in ArcGIS REST Services

texts:
In ArcGIS REST Services Directory when querying a service the values that are set as floats or doubles (elevation, lat, long) return with too many significant digits if the format is JSON. For example if an elevation is 189.9 (and set to have only two decimal places in the attribute table of the original) GET JSON returns 189.9000000001. This seems to only happen where the last digit is a 0. The same thing happens with lat, long. In order to share this service with an outside user they need less significant digits returned. Is there a way to set the significant digits regardless of the 0 in the hundredths place when returning XML?

Returns:{
 "displayFieldName": "County",
 "fieldAliases": {
  "OBJECTID": "OBJECTID",
  "County": "County",
  "USGS_____F": "USGS_____F",
  "Date": "Date",
  "Time": "Time",
  "Depth_to_W": "Depth_to_W",
  "Water_Leve": "Water_Leve",
  "Method_of": "Method_of",
  "Accuracy_o": "Accuracy_o",
  "Lat": "Lat",
  "Long": "Long"
 },
 "geometryType": "esriGeometryPoint",
 "spatialReference": {
  "wkid": 102100,
  "latestWkid": 3857
 },
 "fields": [
  {
   "name": "OBJECTID",
   "type": "esriFieldTypeOID",
   "alias": "OBJECTID"
  },
  {
   "name": "County",
   "type": "esriFieldTypeString",
   "alias": "County",
   "length": 254
  },
  {
   "name": "USGS_____F",
   "type": "esriFieldTypeString",
   "alias": "USGS_____F",
   "length": 254
  },
  {
   "name": "Date",
   "type": "esriFieldTypeDate",
   "alias": "Date",
   "length": 8
  },
  {
   "name": "Time",
   "type": "esriFieldTypeString",
   "alias": "Time",
   "length": 254
  },
  {
   "name": "Depth_to_W",
   "type": "esriFieldTypeDouble",
   "alias": "Depth_to_W"
  },
  {
   "name": "Water_Leve",
   "type": "esriFieldTypeString",
   "alias": "Water_Leve",
   "length": 254
  },
  {
   "name": "Method_of",
   "type": "esriFieldTypeString",
   "alias": "Method_of",
   "length": 254
  },
  {
   "name": "Accuracy_o",
   "type": "esriFieldTypeString",
   "alias": "Accuracy_o",
   "length": 254
  },
  {
   "name": "Lat",
   "type": "esriFieldTypeDouble",
   "alias": "Lat"
  },
  {
   "name": "Long",
   "type": "esriFieldTypeDouble",
   "alias": "Long"
  }
 ],
 "features": [
  {
   "attributes": {
    "OBJECTID": 1110,
    "County": "HOLMES",
    "USGS_____F": "051N0010",
    "Date": 1247097600000,
    "Time": "UNKNOWN CST",
    "Depth_to_W": 189.90000000000001,
    "Water_Leve": "Feet",
    "Method_of": "Tape Down",
    "Accuracy_o": "±0.2 ft",
    "Lat": 33.102849999999997,
    "Long": -89.915080000000003
   },
   "geometry": {
    "x": -10009301.557499999,
    "y": 3908964.3601000011
   }
  }
 ]
}


texts:
How we can check if feature has intersections or not with other features in one shapefile data?

I want to make sure if in that file there is no overlaps between polygon, I have read this approach "An Approach for Checking Overlaps and Gaps in Polygons using Geopandas" and I know that work, but is there any better solutions? With these way the process take time because the loop.

data_overlaps=gpd.GeoDataFrame(crs=data_temp.crs)
for index, row in data_temp.iterrows():
    data_temp1=data_temp.loc[data_temp.ID!=row.ID,]
    overlaps=data_temp1[data_temp1.geometry.overlaps(row.geometry)]['ID'].tolist()
    if len(overlaps)>0:
        for j in overlaps:
            temp_area=gpd.overlay(data_temp.loc[data_temp.ID==row.ID,],data_temp.loc[data_temp.ID==row.ID,],how='intersection')
            temp_area=temp_area.loc[temp_area.geometry.area>=9e-9]
            if temp_area.shape[0]>0:
                data_overlaps=gpd.GeoDataFrame(pd.concat([temp_area,data_overlaps],ignore_index=True),crs=data_temp.crs)


texts: Connecting QGIS to ODBC database (MS Access *.accdb)

texts:
Has anybody managed to successfully connect to open a MS Access *.accdb database from QGIS using ODBC?

All the Q&A, tutorials, guides and videos I've found work with *.mdb files and ESRI geodatabases.  I tried *.accdb but QGIS gets frozen.

Details:

Testing environment

Windows 7 64bit
MS Access 2016 32bit
ODBC drivers used: MS Access Database Engine 2016 32 bit (*.mdb, *.accdb)
Using ODBC 32bit config in Windows (C:\Windows\SysWOW64\odbcad32.exe)
Using QGIS 2.14.16 32bit

In QGIS, I use:

Add vector layer > Database > Type ODBC > New
Type ODBC  >  Name (name1)  >  Database Name (nameA)
Test Connect:  "connection to x was successful"
Press OK  -->  QGIS gets frozen

I can read the same ODBC database from other applications (SQLite, MapInfo), so there is no problem with the ODBC side of it.  It looks like something's not right on the QGIS side.

Does anybody know if this connection should work / why is it not working?


texts: ArcMap Attribute Table not seen when second screen disconnected from Main Computer?

texts: I am working with ArcGIS Map and using two computer screens. I use the second screen to see my attribute table clearly. But once I disconnect my main computer from the second screen and want to see the attribute table I am not able to see it. Even I restart the computer and try, it is opened somewhere. Again, when I connect my main computer with the second screen, I can see my attribute table there. It creates problem when I have to show my work to team on projector far from my second screen.


texts: Invalid result in coordinate transformation

texts:
I want to convert local coordinate (SWEREF99) to  Global coordinate(WGS 84). I tried 3 different ways and all produce almost same result which is not correct.
From here I got actual result: https://www.lantmateriet.se/sv/Sjalvservice/enkel-koordinattransformation/

Input: 6589150 1677500

Expected output: 59.21888889 18.92916667

EGSG codes:

EPSG:3011 : SWEREF99 18 00

EPSG 4326 : WGS 84

But what I got form different options:

1. https://epsg.io/

Output: 68.45944444 9.87361111

2. proj4js

Output: 68.45962374786117, 9.873651976751825

      var firstProjection = '+proj=tmerc +lat_0=0 +lon_0=18 +k=1 +x_0=150000 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs';
      var secondProjection = '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs';
      var output = proj4(firstProjection, secondProjection, [6589150, 1677500]);

2. GDAL

Output: 68.8978437634193 9.65669736119052

      gdaltransform -s_srs EPSG:3011 -t_srs EPSG:4326
      6589150 1677500


texts: Calculate the Intersect Area of Polygon Features and Circles of Equal Area Based on Feature Centroid

texts:
For this task I am using ArcGIS Pro 2.1.  I am trying to fist create a process then automate the process to do the following:

Create circular polygons of equal area to and with the same centroid as the polygon features.
Determine the intersection of the original features and the new circular polygons.
Solve for Bottemma's formula for each feature and write that solution to the original feature class in one field.  Bottema's formula is here:

I have been able to calculate the centroid of each feature.  Next, I have created buffers using a field which calculates the radius of each buffer based on the area of each feature.  My process begins to break down when I run the intersect tool.  As seen below, the resulting intersect is a multipart feature.  The highlighted polygon represents one of the original polygon features.  Additionally, the buffers often overlap each other which creates problems with the intersections.

It seems like every feature needs to be run through this process individually in order for this to work.

Are their any built in tools to ArcGIS that might also simplify this process?


texts: (OSMNX, Folium) Plotting and Isochrone on a Folium Web Map

texts:
I am trying to plot an isochrone map like the one in this example onto a Folium Web Map.

I understand how to set up the colors for time/edges/nodes/polygons... I just don't understand how I could get any of that onto a Folium map...

I am having trouble understanding which object(s) to pass into the following two functions:

osmnx.plot.plot_route_folium()
osmnx.plot.plot_graph_folium()

My code, currently, is simply following the above cited example, but I an stuck on how to incorporate Folium.


texts: Geowebcache is not using bounding box parameter when seeding

texts:
I have an issue where when I am seeding tiles for a layer with Geowebcache, the tiles are seeded to the full extent of the layer despite my use of the bounding box parameter. My goal, of course, is to constrain where the tiles are generated to a specific extent.
I have seen the bounding box parameter specified in two ways (I am using JSON format for both):

This page (https://docs.geoserver.org/latest/en/user/geowebcache/rest/seed.html) suggests using following format:

'bounds':{'coords':{ 'double':['-124.0','22.0','66.0','72.0']}}

The web interface (ie. https://server.org/geoserver/gwc/rest/seed/<layer_name>) seems to send the bounding box parameter is a different format:

         "minX": -1100000,
         "minY": 5000000,
         "maxX": -9500000,
         "maxY": 6600000,

Neither of these have worked for me.
Could someone please advise me on how to specify the bounding box when seeding tiles?
EDIT
I should clarify that the tiles are generated fine. But they are generated outside of the bounding box parameter I specify. Indeed, they seem to obey the default extent (or bounding box) of the layer itself.
So, how can I constrain the tile generation to a specific extent?


texts:
I have to calculate the decline in snow in a mountain area of Peru (between two points = two mountains). This is the link to what I have. https://code.earthengine.google.com/?scriptPath=users%2Fromeoliekens_GIS%2FEx6%3AEXTRA
Next I would like to calculate the pixels that account for snow (white pixels) from the two layers named Snow Cover2000 Class and Snow Cover2020 Class. I tried to mask the layers first but I was unable to to this. Can someone help me calculating the area snow covered in both years (2000 and 2020)?
I used MOD10A1.006 Terra Snow Cover Daily Global 500m as image collection.
var POI1 = ee.Geometry.Point(-77.3162, -9.5734)
var POI2 = ee.Geometry.Point(-77.3141, -9.3953)

Map.addLayer(POI1, {color:'red'}, 'Point 1' )
Map.addLayer(POI2, {color:'blue'}, 'Point 2' )
////
var dataset = ee.ImageCollection('MODIS/006/MOD10A1')
                  .filter(ee.Filter.date('2020-01-01', '2020-12-31'));
var snowCover2020 = dataset.select('NDSI_Snow_Cover');
var snowCoverVis = {
  min: 0.0,
  max: 100.0,
  palette: ['black', '0dffff', '0524ff', 'ffffff'],
};
print(snowCover2020)
Map.addLayer(snowCover2020, snowCoverVis, 'Snow Cover 2020')

/////////////////////////////////////////
/////////////////////////////
var dataset = ee.ImageCollection('MODIS/006/MOD10A1')
                  .filter(ee.Filter.date('2000-01-01', '2000-12-31'));
var snowCover2000 = dataset.select('NDSI_Snow_Cover');
var snowCoverVis = {
  min: 0.0,
  max: 100.0,
  palette: ['black', '0dffff', '0524ff', 'ffffff'],
};
print(snowCover2000)
Map.setCenter(-77.3103, -9.566, 9);
Map.addLayer(snowCover2000, snowCoverVis, 'Snow Cover2000')
//

var clipImage2000 = function(snowCover2000){
  var im_clipped=snowCover2000.clip(geometry)
  return im_clipped
}

var snowCover2000Clipped = snowCover2000.map(clipImage2000)
Map.addLayer(snowCover2000Clipped, snowCoverVis, 'Snow Cover2000 clipped')
/////////////
var clipImage2020 = function(snowCover2020){
  var im_clipped=snowCover2020.clip(geometry)
  return im_clipped
}

var snowCover2020Clipped = snowCover2020.map(clipImage2000)
Map.addLayer(snowCover2020Clipped, snowCoverVis, 'Snow Cover2020 clipped')
//Classify
var snowCoverVis2 = {min: -1, max: 1, palette:['black',  'white',] }
var calcSnow = function(snowCover2020Clipped){
  var snow = snowCover2020Clipped.gt(0)
  return snow
}
var snowCover2020Class = snowCover2020.map(clipImage2020).map(calcSnow)
Map.addLayer(snowCover2020Class, snowCoverVis2, 'Snow Cover2020 Class')

var calcSnow = function(snowCover2000Clipped){
  var snow = snowCover2000Clipped.gt(0)
  return snow
}
var snowCover2000Class = snowCover2000.map(clipImage2000).map(calcSnow)
Map.addLayer(snowCover2000Class, snowCoverVis2, 'Snow Cover2000 Class')


texts: Finding Clusters of N Points Within X Distance using ArcGIS Desktop?

texts:
We have tried two methods of cluster analysis:

Cluster and Outlier
Hot Spot

... but neither meets our specification which has these 3 parameters:

a Set of Point features
N points within a Cluster
X distance band to look in

For example "find me the points where each cluster has 10 points within 50 meters" or something to that effect.

ArcGIS 10 has the two tools I’ve linked above but neither seems to address the N points within a cluster part.

Can our requirements be met with other tools in ArcGIS?
Is what I am calling cluster analysis known as a particular clustering type?

We did have a version of this in another application but I really wanted to use the tools within ArcGIS.


texts: Offsetting points from lines using SQL Server Spatial?

texts: I'm new with SQL spatial data processing, so I'm not sure how to solve my task. I have two tables - table with lines and table with points that are located on lines. Lines and points are related by ID value. I need to create another one table with points whee location of points is shifted for specified distance off corresponding line to the left or right (defined as attribute of a point). I.e. like it is shown below.


texts: Is there a way to Query what's "NOT THERE"?

texts:
If you have land tenure mapping, and you're only interested in the areas for which there are no polygons in a region, is there a way to query these areas out and create a shapefile from them.

A tangible example using Mining Claim Mapping:

Various forms of land tenure are represented for an area.  These areas are captured with polygon geometry.  Areas include Mining Claims, and various forms of Patented Mining Lands (ie - Surface rights only, surface & mineral rights, 11 year lease, etc, etc).

In this example, If I'm ONLY interested in the areas that are not captured (which would be Crown Land), how do I go about extracting that data (data which is "NOT THERE") without manually digitizing?  Is there a way to do this with a select by location query?  The only thing I can come up with is to use the bpoly command in AutoCAD, but this would be a bit time consuming and would also lend itself to the possibility of oversights/errors on the part of the user.


texts: Which OSM2PO method returns segment speed?

texts: In OSM2PO, RoutingResultSegment has a getH() method - what values does it return? I was looking for a method that returns kmh, similarly to what is stored in the DB, but it doesn't seem like getH() returns that...


texts: Recording the ArcPy script execution method

texts:
I have composed an ArcPy script which is run via a windows scheduler.The same script is loaded into a script tool so a user can run the process manually. I've used: get parameters as text, with or's and not's, to hard-wire the standard variables, if they are not speicifed.

ReportFolder = arcpy.GetParameterAsText(0)
if ReportFolder == '#' or not ReportFolder:
    ReportFolder = "C:\\Data\\GIS"

The process runs and during so writes to a text file log, for example:

txtFile.write("=    For ArcGIS 10.3.1: Date: "+str(timed)),txtFile.write ('\n')

I'd like to record what method was used to execute the script; was it via the windows scheduler, or by the script tool, or by a python client like PyScripter.

Is anyone aware of some form of os environment thingy that can be called by ArcPy?


texts: Renaming ArcGIS Pro project, APRX, file geodatabase and toolbox?

texts:
Using ArcGIS Pro 1.3.1 I created a project called TestProject in its own folder.

As expected an examination of the TestProject folder in Windows Explorer showed that it contained TestProject.aprx, TestProject.gdb and TestProject.tbx (as well as Index and ImportLog subfolders).

I proceeded to do considerable work in this project before deciding to rename my project from TestProject to ProperProject.

Renaming the project folder from TestProject to ProperProject using Windows Explorer worked as expected but left the old names of TestProject.aprx, TestProject.gdb and TestProject.tbx within it, as expected.

I was hoping to rename these as ProperProject.aprx, ProperProject.gdb and ProperProject.tbx so I tried to do this two ways:

Using Windows Explorer to rename the three items but when I then opened ProperProject.aprx it resulted in a project with TestProject.gdb (empty), TestProject.tbx (empty), ProperProject.gdb and ProperProject.tbx
Using Windows Explorer to rename TestProject.aprx to ProperProject.aprx and then opening ProperProject.aprx to use the Project pane to try and perform the remaining two renames:

a. Renaming TestProject.gdb to ProperProject.gdb seemed to work as would be expected from the help that covers that in Rename project items

b. No option to rename TestProject.tbx to ProperProject.tbx was available.  I cannot find this mentioned in the help.

c. When I closed (without saving) and re-opened ProperProject.aprx it once again created a new TestProject.gdb (empty)

Is it possible to rename an ArcGIS Pro project in a way that also renames the project's *.aprx, *.gdb and *.tbx to the new project name?


texts: How to make a uniformly transparent layer of overlapping features?

texts:
In a map used for bicycle navigation, I wanted to include highways because they are important landmarks/barriers. But I wanted to de-emphasize them relative to the street network so I made them partially transparent. The problem is that overlapping links of these roads brighten as per the image below, when I would like them to just be uniformly neutral. To clarify, Martha Road is bikeable, and therefore at 100% alpha, while a series of tunnels and overpasses directly to its north are unbikeable, and therefore transparent. The two types of road are in the same layer, with a rule-based styling.

What symbology properties do I have to modify so that overlapping transparent features blend into each other?


texts: How to calculate turn rate from a GPS log?

texts:
I'd like to calculate a rate of turn from a GPS log (GPX file).

Each trackpoint in a GPX track has location and time, and in addition to speed (a common calculation and sometimes embedded in a tracklog) I'd like to calculate the rate at which the aircraft or vehicle is turning, for each trackpoint. That said, I suspect the 'rate of turn' is more accurately reflected by a moving average over several points, and would be interested to see if anyone has examined that.

Sample GPX file here.

Open-source solutions in R or QGIS preferred!


texts: Clipping a raster layer with a shapefile on mac

texts:
I am trying to clip a raster layer using the raster clipper function.
However, for some reason this is not working for me and I haven't been able to find the answer anywhere.

I choose my raster layer as the input file, my shapefile as the mask layer and creates an output file, but after pressing OK nothing happens.

I am using QGIS 2.6.0 and following the steps presented at: http://www.qgistutorials.com/en/docs/raster_mosaicing_and_clipping.html

Does anybody have an idea of what I might be doing wrong?


texts: Trouble Converting from a shp. to gpx file and up loading to GPS

texts:
Currently I have a series of random points within a shapefile that I would like to upload to my GPS. Every time I try to upload the points via "Save file" it come up with an error message (posted below).

I have tried turning Force_gpx_Track=Yes and it still does not work.
I have also tried unclicking the skip attribute creation and nothing works.

It worked one time and then come up with a select file class with waypoints, lines, and tracks so I clicked on waypoints and it ran the file.

Once I try to add the file to my GPS (Garmin GPSMAP 64) using GPS tools it only comes up with tracks under Upload to GPS any help would be great because it have reached a road block and need to have this inventory done.

Error Message

Export to vector file failed.
Error: creation of data source failed (OGR error:Failed to create GPX file Inventory_Points .gpx.

Still having trouble converting SHP to GPX, had it work once and never again. I created the name field (seen below) and copied the folder into my Garmin 64 via copy and past but it shows up blank. The one time it worked was with random points, now I am trying to get it to work with a grid

In the post above I checked skip attribute creation and forced GPX tracks with still no prevail. There is also no blanks at the end of the inventory_points

Was also wondering if I could use a program called basecamp or DNR_gps to export the file that way, if so could someone post the directions?


texts: TIGER/Line GeoID - Census Tract map to Lat/Long

texts:
I'm working in Python with ACS data downloaded from here:  https://www.census.gov/geo/maps-data/data/tiger-data.html

I was able to convert the geodatabase format to a pandas dataframe where every row is a GeoID representing a census tract.

I can't figure out how to understand the GeoID and to connect it to latitude/longitude.  My end goal is to build a function that takes in a lat/lon and returns what census tract it is in.

For example, this library: https://pypi.org/project/censusgeocode/ will take a lat/lon and return the census tract information.  But it returns (as an example):

'GEOID': '42079216600'

whereas the TIGER/Line data I downloaded has the GeoID in a different format:

GEOID: '15000US170010001001'

Can someone help me understand how to interpret these GeoID's so I can join on them? I've looked here: https://www.census.gov/geo/reference/geoidentifiers.html but it doesn't explain down to the census tract level.


texts:
I'm going through some of the javascript tutorials - build your first application. It works fine when I create an html file that holds the style and script tags but when I separate the files out to html, js and css the map can only be seen over half the screen. I've tried a few things but I can't get the map to full extent like it should be...should I tweak the css?

html, body, #map {
       height: 100%;
       width: 100%;
       margin: 0;
       padding: 0;
     }
     body {
    background-color: #FFF;
    overflow: hidden;
    font-family: "Trebuchet MS";
  }


texts: Extract OSM data (xml) from large polygon (up to 1GB of data)

texts:
I need to extract a large area from OSM to import it into my local PostGIS-DB (using osm2pgsql).

For smaller areas I was using http://extract.bbbike.org/ to download OSM-XML-Data. It's possible to select the area by rectangle or by polygon. The maximum size of extract can be 512MB. But now I need to download a larger area - maybe about 1GB - too large for bbbike.org.

QuickOSM seems to be also no option, cause of the size of the area.

The area of interest is a 30km buffer of a 900km-track and runs across 2 Countrys.  I could download both Countrys, but that would by a way too much data for my available resources.

How can I get to my data?


texts: Calculating Length of Line features using GDAL, Shapely libraries of Python

texts: I have a line shapefile named "river" which has 385 features. I would like to calculate the length of each feature using Python. I am currently using GDAL, Shapely, and Fiona libraries.


texts: Is it possible to save a KML layer's style? My layer is set to 50 % opacity and I want to save that setting into the KML-file. The plugin qgistokml can´t yet do it, and only work with points.


texts: Selecting features from shapefile where field="value" using Geopandas?

texts:
I don't understand the difference between OGR and Geopandas librairies.
I have checked the documentation of GDAL as well as http://geopandas.org/ but I still don't get it.
I just started learning about them and for now I'm trying to write some scripts in python using Geopandas.

My goal is to select features from a shapefile where field="value", then export the selected features to a shapefile.
The second step would be to add another field to the output shapefile(with only the selected features from the first shapefile) and set it to a certain value.
Doing the second part with OGR seems easy. This is how I've done it :

from osgeo import ogr
driver = ogr.GetDriverByName('ESRI Shapefile')
dataSrc = driver.Open("myshp.shp", 1)
fld = ogr.FieldDefn('City', ogr.OFTString)
lyr = dataSrc.GetLayer()
lyr.CreateField(fld)
for feature in lyr:
    feature.SetField("Not Found", fld)
    layer.SetFeature(feature)
dataSrc= None

However I don't see how I can do the same thing using Geopandas. I might be missing the point though that's why I'm confused between the use of the two librairies.

I know that this is easily done using QGIS or ArcGIS manually, but as of now I'm trying to learn how to use python with Geopandas.

How can I achieve my goal explained above using Geopandas?


texts: Polygons with WKT extension load extremely slow

texts:
I am quite new to QGIS. I have temporal data per month for the period of 8 years for all Indian districts (+/- 650 districts). This means I have approx. 65.000 objects in my attributes table. Once I upload my file as a layer (a text delimited file) in QGIS using WKT-text, the polyons load extremely slow.
Is there anything I can do about this? I would ultimately like to use the time-manager to visualize my temporal data on India.


texts: How to pin FEATURE ATTRIBUTE FORM window in QGIS?

texts:
When I identify a feature in QGIS the window of IDENTIFY FEAtURE is open and the FEatURE ATTRIBUTE FORM window too.

The form of the attribute table is open in the middle of the QGIS window and I want it to be open somewhere in the corner.

Every time I click to identify new feature - the attribute form is opening in a new window. How can I prevent this?


texts: Problem with SPIT (error while importing layers to POSTGIS)

texts:
I was attempting to import a polyline roads layer into postgis, using both SPIT in QGIS 1.7 and postGIS GUI for this purpose. I have received errors through both processes as described here:

SPIT

  The database gave an error while executing this SQL:
  INSERT INTO "public"."highways"("STREET","FROMLEFT","TOLEFT","FROMRIGHT","TORIGHT","STREETCODE","OBJID","USERID","FJUNCTION","TJUNCTION","LENGTH","ROA... (rest of SQL trimmed)
  The error was:
  ERROR:  value too long for type character varying(31)

postGIS GUI

What can be the problem and how can I fix it?

Thanks for your response.

Dor


texts: Updating field of selected results in PyQGIS?

texts:
I'm trying to write a script that will update a field with 2 values (Wrong/Ok) when specific conditions are met.
This is what I've come up so far, but I can't get it to work..

from qgis.core import *
from qgis.gui import *

layers = QgsMapLayerRegistry.instance().mapLayersByName('Mukono_Polygon')
layer = layers[0]
it = layer.getFeatures()
idx = layer.fieldNameIndex('PlotArea')
idy = layer.fieldNameIndex('Pr')
idz = layer.fieldNameIndex('CTRL')

layer.startEditing()
for feat in it:
  if feat.attributes()[idx] <> '888.8888' and feat.attributes()[idx] <> '999.9999' and feat.attributes()[idy] > '120':
    layer.changeAttributeValue(feat.attributes()[idz], 'Wrong')
  elif feat.attributes()[idx] <> '888.8888' and feat.attributes()[idx] <> '999.9999' and feat.attributes()[idy] > '80':
    layer.changeAttributeValue(feat.attributes()[idz], 'Ok')
  else:
    print 'Cant find specific entries'

layer.commitChanges()

Traceback (most recent call last):
  File "", line 23, in
IOError: [Errno 9] Bad file descriptor

If remove change the line 23 to print '', the script doesnt do anything.

I chose to get the fields by name instead of their number index because it's different from shapefile to shapefile.


texts: Select columns to keep from attribute table

texts:
I have lots of feature classes that I'm working with, but the attribute tables are really messy because I got the information from different sources. So I decided to clean them up.

I made a table for every criteria one row. I want to use the table with an iterator, so all the features are imported, renamed and sorted out at once.
I added ONE column (called "Keep_it") with the name of the column that I want to keep from the old feature - all the other columns should be deleted.

It's working perfectly if I only want to keep one column from every feature, but for some of my features I want to keep two or three columns.

I tried to add a second column, but it was not working, because it was empty for most of my features. I also tried to add more than one column name in the "keep_it_2" column... also not working.

Here you can see my Table:

So for ObjectID 1 I want to keep the Field named "WSGNR" and "WSG_Name". Just keeping "WSGNR" from "keep it" works - but to combine it with a second field is not working.

Any ideas?


texts: Extracting information from overlapping shapefiles

texts:
I'm working with two shapefiles (geometries): Census tract data and middle school attendance zones. While both shapefiles obviously have geometries, the census tract data also has substantial amount of zipcode-level income information (things like median family income).
These middle school attendance zones are completely irregular and do not follow any zipcode boundary lines, which means that each middle school attendance zone can have different portions of multiple zip codes comprising the attendance zone. Looking at a certain county, which has a number of middle schools and many zip codes, I would like to have an idea of the average wealth of any particular middle school family.
I need to use the zipcode-level Census data on family income but will need to approximate the average family income of a middle school family by taking a weighted average of the different zipcodes that feed into that middle school attendance zone.
Thus, I want to overlay the two shapefiles and extract the proportion of the every zipcode that is present in every middle school attendance zone. I have attempted to use the GeoPandas overlay function with intersect but since this is a tricky problem and I'm relatively new to this, I didn't get to the answer I was looking for...or maybe I did but I'm just too confused to understand it.


texts: Find pairs of near features with geopandas, fiona, GDAL etc

texts:
I have a feature class of points, a few thousands I have now, but I expect to deal with hundreds of thousands of points.  I would like find pair of points that is within a specific distance.  The points I have is typically spread apart spatially so I don't expect to find many pairs.  How I can find those pairs.

Right now the points are originally from ESRI shapefile, and read into geopandas geodataframe.  I tried to tbl = df['geometry'].distance(df['geometry']) see if I can get distance of all pairs.

But I got the error "python: GeometryComponentFilter.cpp:34: virtual void geos::geom::GeometryComponentFilter::filter_ro(const geos::geom::Geometry*): Assertion `0' failed.".  Beside

But, I don't want to get full set of distance as I am only interested in points that are close to each other.

I have a version of code which uses ArcGIS, and GenerateNearTable_analysis is what I used for the purpose.

I want to run code with open source code now.

My solution has to be python and should run on Linux.


texts: Looking for system specifications for a modern GIS workstation

texts:
My current employer has given me the opportunity to design my own GIS workstation. The problem is I have no idea what is important when it comes to hardware! I have been given a datasheet from Fujitsu which includes the latest CELSIUS M720 series workstations but the configuration includes about 30 different processors, graphics cards, and hard drives.

I will be working a lot with the ESRI desktop suite (spatial analyst and 3D analyst) and doing a fair amount of raster processing (DEMs, visibility analysis, etc.), but I also use the GDAL library, SAGA-GIS, Quantum GIS and others.  Does ArcGIS 10 actually use multiple cores?

I have the choice from 4 cores with 8 threads up to 8 cores with 8 threads - 16 threads?

Is SSD SATA III better than HDD SATA III?

The Celsius R920 supports up to 512 GB RAM. But what do I need?

It will probably be the last time I get a new workstation for the next 5 years, so I am curious what others have ordered or feel will be necessary for the coming years. Obviously the cloud isn't taking over as quick as they say!


texts:
I am running GeoServer 2.4.1 (Recently Upgraded)

Server: Windows 2008R2, 2Gb Ram (Has been running perfectly until now)

My WMS Rendering has slowed down significantly. I am unsure if this is a problem with the Upgrade or the configuration. We added a number of additional layers in Layer Groups in the last few days.

The problem is intermittent, some tiles will render, some will not. I have noticed that the Server goes to 100% CPU and 90% Ram use and this coincides with the tile rendering failure.

Data is being pulled form a PostGres Db on another Computer, Network load is insignificant.

Question: Is this a configuration issue?
Question; How can I better tune the GeoServer Configuration

This is the Relevant Error from the Geoserver logs:

2013-11-07 10:56:23,744 INFO [geoserver.wms] -
Request: getServiceInfo
2013-11-07 10:56:24,384 WARN [geotools.rendering] - Assuming rendering buffer = 1, but estimation is not accurate, you may want to set a buffer manually
2013-11-07 10:56:34,806 INFO [geoserver.wms] -
Request: getServiceInfo
2013-11-07 10:56:34,837 WARN [geotools.rendering] - Assuming rendering buffer = 4, but estimation is not accurate, you may want to set a buffer manually
2013-11-07 10:56:34,962 ERROR [geotools.map] - Call MapContent dispose() to prevent memory leaks
2013-11-07 10:56:34,962 INFO [geoserver.wms] -
Request: getServiceInfo
2013-11-07 10:56:35,165 WARN [geotools.rendering] - Assuming rendering buffer = 4, but estimation is not accurate, you may want to set a buffer manually
2013-11-07 10:56:47,848 ERROR [geotools.rendering] - JDBC DataSource not available after dispose() has been called
java.lang.NullPointerException: JDBC DataSource not available after dispose() has been called

I have a Max Rendering Memory of 100000 (when this error was generated) but have tried from 16000 and steps in between with no differences.

Geoserver Configuration as follows:

JVM config as follows (Using Jdk1.7.0_40)

# Java Classpath (include wrapper.jar)  Add class path elements as
#  needed starting from 1
wrapper.java.classpath.1=wrapper/lib/wrapper.jar
wrapper.java.classpath.2=start.jar
wrapper.java.classpath.3=lib/*.jar

# Java Library Path (location of Wrapper.DLL or libwrapper.so)
wrapper.java.library.path.1=wrapper/lib

# Java Additional Parameters
wrapper.java.additional.1=-Djetty.home=.
wrapper.java.additional.2=-DGEOSERVER_DATA_DIR="%GEOSERVER_DATA_DIR%"
wrapper.java.additional.3=-XX:MaxPermSize=128m
wrapper.java.additional.4=-XX:SoftRefLRUPolicyMSPerMB=36000
wrapper.java.additional.5=-Xmx1024M
wrapper.java.additional.6=-Xms512M
wrapper.java.additional.7=-XX:+UseParallelGC
wrapper.java.additional.8=-XX:NewRatio=2
wrapper.java.additional.9=-server

# Initial Java Heap Size (in MB)
wrapper.java.initmemory=512

# Maximum Java Heap Size (in MB)
wrapper.java.maxmemory=1024

# Application parameters.  Add parameters as needed starting from 1
wrapper.app.parameter.1=org.mortbay.start.Main
wrapper.app.parameter.2=etc/jetty.xml


texts:
I was inspired by this wikipedia image from the mercator projections page to hack around with leaflet. The idea was to interactively have your cursor move a continent and that it would be corrected for the projection.

See the following JSFiddle for the status sofar;

http://jsfiddle.net/cantdutchthis/NrAEL/1/

I am a bit new in the projections world and although the demo works well enough for the latitude. When I compare the equirectangular and the mercator projection on wikipedia I get the impression that the projection's bias is only on the 'y-axis' when you compare these two projections.

I get the impression that I am not correcting for the longitude due to the first wikipedia picture. I would like to support this but I don't know what the best way to do this is.

Currently the way I am able to show the bias of the mercator projection is by shifting the coordinates to the mousecursor everytime the mouse moves;

map.on("mousemove", function(d){
    Object.keys(myLayer._layers).map(function(d){ myLayer.removeLayer( myLayer._layers[d] ) })

    var plotCoords = [];

    // green is the coordinate array for the landmass
    var green = [ [-46.76373672, 82.62796021],... ]

    for ( i in green ) {
       // shift the landmass to where the cursor is
       var x = d.latlng.lng + green[i][0] + 40
       var y = d.latlng.lat + green[i][1] - 80
       circleCoords.push( [ x,y ] );
    }

    window.coords = plotCoords

    window.geojsonFeature = {
      "type": "Feature",
      "properties": {},
      "geometry": {
          "type": "Polygon",
          "coordinates": [coords]
      }
    };

    myLayer.addData(geojsonFeature);

});

What would be the best way to have this interactive leaflet plot also show the effect as shown on wikipedia.

-- Edit!

For more examples see my blog; http://koaning.com


texts:
I am trying to write my own layerswitcher because I want a dropdown box (and later a second one)

My Webinspector says: Uncaught TypeError: Cannot read property 'value' of undefined
but I cannot find what is wrong. In my opinion it is defined.
Thanks to everyone trying to help me.

Here is the code:

<select onchange="auswahlLayer1()">
        <option value="Einwohner">Einwohner</option>
        <option value="Energieverbrauch">Energieverbrauch</option>
</select>


texts: Converting NAD83 (epsg:4269) to WGS84 (epsg:4326) using pyproj

texts:
I'm trying to convert US census data (NAD83) to WGS84.
I know the difference between the two is quite minor, but trying to apply the transformation using pyproj does not do anything. e.g.:

from pyproj import Proj, transform
p1 = Proj(init='epsg:4269')
p2 = Proj(init='epsg:4326')
transform(p1, p2 , -80.00001, 40.00001)

yields:

(-80.00001, 40.00001000000001)

Any idea why? Shouldn't I expect small shifts of ~1 meter?


texts: On the calculus of the variable "autocov" from spdep package

texts:
I’m facing a spatial autocorrelation problem  and I want to calculate the «  autocov» covariate through the function autocov_dist from spdep package. The problem is when I defined the neighborhood radius as twice the minimum distance to the nearest neighbor bw2as follow :
library(spded)
        knn <- spdep::knn2nb(spdep::knearneigh(xy, k=1)) # xy are in UTM
        bw <- max(unlist(spdep::nbdists(knn, xy)))
        bw2 <- 2*min(unlist(spdep::nbdists(knn, xy)))
        AUTOCOV <- spdep::autocov_dist(Data_Species[,"P_A"], xy = xy, nbs = bw2, style = "B", type = "inverse")

I got this error message :

Error in nb2listw(nb, glist = gl, style = style, zero.policy = zero.policy) :   Empty neighbour sets foundDe plus : Warning message:In spdep::autocov_dist(Data_Species[, "P_A"], xy = xy, nbs = bw2,  :  With value 135.529724518632 some points have no neighbours

and I was obliged to use the maximum distance to the nearest neighbor bw to avoid this error but the radius here is too large. Is there any way how to manage the situation where some points have no neighbors.


texts:
I have some polygons and multipolygons created using latitudes and longitudes using shapely(there is no shapefile). I want to plot them using basemap. However they don't show up. When I use matplotlib without using basemap, they show up just fine. But they don't work alongside basemap:

fig = plt.figure()
ax = fig.add_subplot(111)
w, h = maxx - minx, maxy - miny
ax.set_xlim(minx - 0.2 * w, maxx + 0.2 * w)
ax.set_ylim(miny - 0.2 * h, maxy + 0.2 * h)
ax.set_aspect(1)

em = Basemap(
    projection='merc',
    llcrnrlon=minx - 0.2 * w,
    llcrnrlat=miny - 0.2 * h,
    urcrnrlon=maxx + 0.2 * w,
    urcrnrlat=maxy + 0.2 * h,
    lat_ts=0,
    resolution='h')

m.drawcoastlines()
m.drawmapboundary()

patches = []
# Initialize patches to array of shapely Polygons using latitude and longitude

ax.add_collection(PatchCollection(patches, facecolor= 'm', edgecolor='k', linewidths=1., zorder=2))
plt.show()


texts: Snapping leaves small gap between features in ArcGIS for Desktop?

texts:
I'm using classic snapping in ArcGIS 10.1.

At the scale (1:5,000) that I'm using to snap the edge of one rectangular feature to the edge of another rectangular feature, snapping seems to work fine. But when I zoom into the "shared" line, I begin to see a gap between the features. I thought snapping was supposed to create topologically adjacent lines.

Is this just some kind of drawing problem in ArcMap and not actually a data problem?


texts:
I'm following the Switch2OSM guide to the tee. I have not had any trouble I couldn't figure out until now. I've typed the following according to instruction:

sudo mkdir /var/run/renderd

sudo chown username /var/run/renderd

sudo -u username renderd -f -c /usr/local/etc/renderd.conf

I get this error:

  renderd[8641]: An error occurred while loading the map layer 'default': Could not create datasource for type: 'shape'  encountered during parsing of layer 'land-low' in Layer at line 48 of '/usr/local/share/maps/style/OSMBright/OSMBright.xml'renderd[8641]: An error occurred while loading the map layer 'default': Could not create datasource for type: 'shape'  encountered during parsing of layer 'land-low' in Layer at line 48 of '/usr/local/share/maps/style/OSMBright/OSMBright.xml'renderd[8641]: An error occurred while loading the map layer 'default': Could not create datasource for type: 'shape'  encountered during parsing of layer 'land-low' in Layer at line 48 of '/usr/local/share/maps/style/OSMBright/OSMBright.xml'


texts: Measurements across UTM zones, how to apply spatial parameters to a customize a projection?

texts:
I have a survey grid that spans 3 UTM zones. Using http://www.epsg-registry.org/ I can get an indication of what projections are available. Where can I find more info on the patterns of distortion within projections? And where can I find information on how I can modify the parameters of these projections (i.e. aspect and parallels)?

Ok So I gather the necessary information for a bounding box: i.e North Latitude, South Latitude, etc...

Let say that I then choose Africa Lambert Conformal Conic as my projection. If I then want to change the parallels of the projection using proj4 - how do I do it?

    +proj=lcc +lat_1=20 +lat_2=-23 +lat_0=0 +lon_0=25 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs

I would want to change this in order to bring the parallels closer to contain the survey grid i.e.:

    +proj=lcc +lat_1=5 +lat_2=-5 ?

Where can I find out if this is appropriate?

What do I need to consider with the rest of the proj4 string?

   +lon_0=25 etc?


texts: I have high resolution images captured by a drone using a , they are RBG and NIR images. My question is how can I calculate the reflectance for the images, the values are currently in digital numbers (DN), I do not know how to proceed. I believe I would need the information from solar radiation (using a sensor) and also from a reflective material with known refluence values, but I'm not sure.


texts: Why does GeoServer copy my PostGIS tables?

texts:
I have some geospatial data in a PostGIS database.  I want to serve that data to a web map app with GeoServer.  But every time I add a table as a layer in GeoServer, it creates a duplicate table in PostGIS.  For example, when I try to add the table "district_boundary", I end up with a new table, named "district_boundary0" that contains all the data from the original table, duplicated for no apparent reason.

I don't want this!  There's no reason to have multiple copies of the same data - and there are lots of reasons not to!  So how can I make GeoServer use the table I tell it to and not make its own?

Thanks in advance!

Note:  This is probably not a matter of my tables not having the correct structure to be served up.  I added the GeoServer-created "district_boundary0" to check this, and just ended up with "district_boundary00"...


texts: skimulator - Python error message on Windows 10 - load_entry_point

texts:
I am running python code called skimulator on Windows 10 to model Arctic ice drift and sea/wave velocities for a satellite mission in 2025. Its installed using Anaconda navigator with a python 3, Numpy, Scipy and NetCDF4 environment. Can anyone help with the first line error message?:

File "C:\Users\pclink\Anaconda3\Scripts\skimulator-script.py", line 11, in <module>
load_entry_point('skimulator==1.0', 'console_scripts', 'skimulator')()

Or even the whole error message?

skimulator-script.py code is shown below.

Error message in Anaconda prompt:

(base) C:\Users\pclink\Anaconda3\Scripts>skimulator C:\Users\pclink\skimulator\example\params_2017v2_23W.py
Traceback (most recent call last):
  File "C:\Users\pclink\Anaconda3\Scripts\skimulator-script.py", line 11, in <module>
    load_entry_point('skimulator==1.0', 'console_scripts', 'skimulator')()
  File "C:\Users\pclink\Anaconda3\lib\site-packages\pkg_resources\__init__.py", line 484, in load_entry_point
    return get_distribution(dist).load_entry_point(group, name)
  File "C:\Users\pclink\Anaconda3\lib\site-packages\pkg_resources\__init__.py", line 2707, in load_entry_point
    return ep.load()
  File "C:\Users\pclink\Anaconda3\lib\site-packages\pkg_resources\__init__.py", line 2325, in load
    return self.resolve()
  File "C:\Users\pclink\Anaconda3\lib\site-packages\pkg_resources\__init__.py", line 2331, in resolve
    module = __import__(self.module_name, fromlist=['__name__'], level=0)
  File "<frozen importlib._bootstrap>", line 983, in _find_and_load
  File "<frozen importlib._bootstrap>", line 967, in _find_and_load_unlocked
  File "<frozen importlib._bootstrap>", line 668, in _load_unlocked
  File "<frozen importlib._bootstrap>", line 638, in _load_backward_compatible
  File "C:\Users\pclink\Anaconda3\lib\site-packages\skimulator-1.0-py3.7.egg\skimulator\__init__.py", line 62, in <module>
  File "C:\Users\pclink\Anaconda3\lib\site-packages\pkg_resources\__init__.py", line 1140, in resource_filename
    self, resource_name
  File "C:\Users\pclink\Anaconda3\lib\site-packages\pkg_resources\__init__.py", line 1647, in get_resource_filename
    return self._extract_resource(manager, zip_path)
  File "C:\Users\pclink\Anaconda3\lib\site-packages\pkg_resources\__init__.py", line 1668, in _extract_resource
    timestamp, size = self._get_date_and_size(self.zipinfo[zip_path])
KeyError: 'skimulator\\share\\VERSION.txt'

skimulator-script.py   :

#!C:\Users\pclink\Anaconda3\python.exe
# EASY-INSTALL-ENTRY-SCRIPT: 'skimulator==1.0','console_scripts','skimulator'
__requires__ = 'skimulator==1.0'
import re
import sys
from pkg_resources import load_entry_point

if __name__ == '__main__':
    sys.argv[0] = re.sub(r'(-script\.pyw?|\.exe)?$', '', sys.argv[0])
    sys.exit(
        load_entry_point('skimulator==1.0', 'console_scripts', 'skimulator')()
    )


texts: OpenLayers WFS no FID being passed from TinyOWS

texts:
I have a webapp - OpenLayers<->tinyows<->postgis. With 1 layer I have a problem.

On localhost all works well - here is a snip of the XML I get back:

 <gml:featureMember>
   <tows:trad gml:id="trad.34">
   <tows:geom><gml:Point srsName="EPSG:4326"><gml:pos>12.01865 57.666245</gml:pos></gml:Point></tows:geom>
   etc...

but on my web server, the XML I get back is:

 <gml:featureMember>
   <tows:trad>
   <tows:ogc_fid>34</tows:ogc_fid>
   <tows:geom><gml:Point srsName="EPSG:4326"><gml:pos>12.01865   57.666245</gml:pos></gml:Point></tows:geom>
   etc...

The layer displays fine, but this quirk means that I am unable to edit the layer:

tinyows error:
    [Tue Mar 18 16:33:07 2014] [ERROR] Element '{http://www.opengis.net/ogc}FeatureId': The attribute 'fid' is required but missing.

/etc/tinyows.xml:

<tinyows online_resource="http://localhost/cgi-bin/tinyows"
     schema_dir="/usr/share/tinyows/schema/"
     log="/var/log/tinyows.log"
     log_level="15" >

<pg host="localhost" user="postgres" password="xxxx" dbname="test" port="5432"/>

<metadata name="TinyOWS Server"
       title="treeMapp TinyOWS Server - Service" />

<layer retrievable="1"
    writable="1"
    ns_prefix="tows"
    ns_uri="http://www.tinyows.org/"
    name="trad"
    table="trad"
    title="Träd" />

<layer retrievable="1"
    writable="1"
    ns_prefix="tows"
    ns_uri="http://www.tinyows.org/"
    name="tradtest"
    table="tradtest"
    title="Träd test" />

</tinyows>

Other layers on the same set up on the web server work fine.


texts: Automatically installing plugins on QGIS, deploying QGIS plugins on an IT infrastructure

texts:
My company is preparing a default QGIS configuration file (qgis_global_settings.ini) for deploying on a great number of MS Windows workstations.

We have selected a few QGIS plugins we want to be available on each worktations.

Is there a way to automatically install these extra plugins when installing QGIS?
By example, can the API be used for that?
Can a network directory be shared to share pre-installed plugins between all workstations?


texts:
I have a tabular data set that contains eastings and northings and LAT and LONG  data. I check some points manually with the ArcGIS "go to xy" point tool and they al appear to be in the correct place on a base map. However, after numerous attempts at bringing the data into arc and using "display xy data" to place them on the map I still end up with the data being projected in the wrong place and to the wrong scale.
What could be causing this, and how do I fix it?

I am operating with ArcGIS, Displaying in Decimal Degrees, WGS 1984 WM (Data frame coord system)

Projected Coordinate System:
  Name: WGS_1984_World_Mercator

Geographic Coordinate System:
  Name: GCS_WGS_1984


texts:
I have downloaded script "Save features filtered by expression" in Processing toolbox ->  Scripts (It is saved in /home/user/.QGIS/processing) and I would like to call this script in standalone PyQGIS script. It works nice if used in normal Ubuntu with QGIS.

I call it with this:

# select by attribute (select water)
input_file=sys.argv[1]
output_file=sys.argv[2]
expression='DN=1'
general.runalg("script:Savefeaturesbyfilteredexpression",input_file,expression,output_file)

But when I try use it with my Docker Image it is not working anymore because that directory does not exist. If I have the file downloaded where I should script copy to and how to call it?

I am asking in general how to use external script. Other processing tools (grass, gdal) works in my Docker Image, but this script has a error: Algorithm not found.

EDIT:
I have copied the script (in docker image) to the same location as in Ubuntu. But it is still not working. Error: Algorithm is not found.

Is there a default folder to put external script and PyQGIS will automatically find it there and will use it or script can be anywhere and I am supposed call script by some function to make it working?


texts: Case studies of Participatory Spatial Planning

texts: I am wondering if you know any well implemented geoportals for Participatory Spatial Planning. I am looking for the one related to environmental protection or connected somehow with the topic (but of course not only- if you have any just put it here).


texts: Getting 999999 error from ArcPy script to increment one field when another field contains 0?

texts:
Before I begin I should emphasise that my python knowledge is pretty limited; generally self taught or inspired by snippets I have found on here.

I have a .dbf table, and I need to use python to increment values (row above + 1) in a field ONLY when a value in another field is 0. For all other rows the value is 1. To better illustrate this, when it is complete the table should look something like this:

FID   Value   Increment
 0      1        0
 1      1        0
 2      1        0
 3      1        0
 4      0        1
 5      0        2
 6      0        3
 7      0        4
 8      1        4
 9      1        4
 10     1        4
 11     1        4
 12     0        5

Using the code from the top answer to Computing new attribute based on changes in another attribute using ArcGIS Desktop with Python?

I edited the code to look like this:

import arcpy

fc = arcpy.GetParameterAsText(0)

rows = arcpy.UpdateCursor(fc)

increment = 0
for row in rows:
    if row.VALUE == 0:
        increment += 1
        row.INCREMENT = increment
        rows.updateRow(row)

    while row.VALUE == 1:
        row.INCREMENT = increment
        rows.updateRow(row)
        rows.next()

    row.INCREMENT = increment
    rows.updateRow(row)

del row, rows

I run it, and it takes quite a while (I have about 11,000 records in the table), before I get this error:

Traceback (most recent call last):
  File "C:\scratch\incrementLine_IncreFieldEverytimeIncrementChanges.py", line 16, in <module>
    rows.updateRow(row)
  File "c:\program files (x86)\arcgis\desktop10.1\arcpy\arcpy\arcobjects\arcobjects.py", line 102, in updateRow
    return convertArcObjectToPythonObject(self._arc_object.UpdateRow(*gp_fixargs(args)))
RuntimeError: ERROR 999999: Error executing function.

The thing is though, the code DOES work, it gives me the result in the table that I want, but due to this error I can't slot this script into my ModelBuilder flow, which is rather frustrating!

So my question is, what is causing this error? As I said, the script is editing the table how I want it, so if there is way I can just stop the error so I can use it in a flow that would be ideal! I imagine my limited knowledge in python (and scripting in general) is making me miss something entirely obvious, and also the fact that I'm basically ripping off a script from 2012 could be introducing out of date techniques?


texts: Create a simple Natural Earth topographic map with Mapbox

texts:
I want embed a Mapbox map on a website and put some placemarks on it (I already did so, but that was a year ago and everything was easier back then ;)). I am a bit confused with the new Mapbox Studio software, one has to use to create custom styled maps with Mapbox now.

What I want is actually quite simple: a map that looks more or less exactly like the first one in this blog post:

https://www.mapbox.com/blog/eight-new-terrain-baselayer-maps-available-mapbox/

So just the Natural Earth terrain baselayer and nothing else. I probably would change the water colour, but that's it.

How do I do that?


texts: Importing polylines and polygons into QGIS from CSV

texts:
I want to import a large number of polylines and/ or polygons into one one or more QGIS layers in order to plot the data.  Each will only have simple attributes (e.g. name).  I need to know:

The required format of the CSV file to enable step 2.
The command or plug in required to create a polyline or polygon layer in QGIS.

For information, the input CSV files are written by bespoke Excel VBA or by FORTRAN codes, so I have full control of these output formats. Essentially I'm looking for the polyline and polygon equivalent of the delimited text plugin.


texts: Clip only works on some areas, and not on others, despite being from the same dataset

texts:
I have a shapefile of provinces (nxprovincias), and another shapefile of small, disconnected areas ("protected areas") that overlaps with the provinces. I want to easily calculate the total area of the protected areas within each province.

I am using QGIS 2.8.2 and all shapefiles are in the same CRS.

I have tried to "join attributes by location" - this gives an empty shapefile.

So because this didn't work, I then split the provinces shapefile up so that I have a different shapefile for each province. Then I used "clip" to clip the protected areas by province. This works perfectly for some provinces, and creates empty shapefiles for the others. Why would this happen?

Shapefiles used can be found here: https://www.dropbox.com/sh/thxvhkkoswyrldg/AAB5Lxh16uhj1WWOKVYtBj29a?dl=0

(In case anyone does try, the provinces that do not work are 4, 9, 11, 13 and 24, which all happen to be in the West of the country, for reasons unknown to me)

My question is:
1. Why does the clip only work sometimes?
2. Is there a better way to do this? I just want to know the area in each province, which seems like it should be easy to calculate?


texts: Why does QGIS 2.16 have longer load times than 2.8, especially working with tables?

texts: I was previously running an LTR version of QGIS (2.8.9 I think). I installed 2.16 and found that opening tables takes a very long time. Deleting columns in tables now takes hours as opposed to seconds. It also seems to load and delete columns one by one instead of all at once. I tried the 2.14 LTR and am experiencing the same problem.


texts: Copying shapefile then making feature layer using ArcPy gives ERROR 000732?

texts:
The only thing I need is an exact copy of an existing shapefile in the same directory with a different name.

I have tried:

arcpy.CopyFeatures_management (ORLinks, MyLinks) #ORLinks: old links, MyLinks: new links

but for some reason, I cant make it a feature layer using:

arcpy.MakeFeatureLayer_management(MyLinks, 'MyLinks')

Error Encountered:

Error: Failed to execute. Parameters are not valid.
ERROR 000732: Input Features: Dataset
M:\RAIL\Rail2.0.shp does not exist or is not supported
Failed to execute (MakeFeatureLayer).


texts: How to calculate the Oriented Bounding Box (OBB) of a polygon or point set?

texts:
I'd like to calculate the "Oriented Bounding Box" (OBB) which is the enclosing rectangle aligned to the longest extent of a point set (which is not necessarily aligned with the axes of the co-ordinate system).
The solution could be similar to Finding minimum-area-rectangle for given points? but I'm looking for a robust solution in PostGIS SQL (and PL/Python if necessary).


texts: Setting new size for printLayout using PyQGIS

texts:
I am trying to create an A0 layout from PyQGIS. I tried the following code but it doesn't work.

from qgis.core import *

manager = QgsProject.instance().layoutManager()
layouts_list = manager.printLayouts()
for l in layouts_list:
    l.clear()

layout = QgsPrintLayout(QgsProject.instance())
layout.setName("TEST")
layout.initializeDefaults()
page = QgsLayoutItemPage(layout)
pageSizes = QgsPageSizeRegistry()
pageSize = pageSizes.find("A0")[0]
orientation = QgsLayoutItemPage.Landscape
valid = page.setPageSize(pageSize.displayName, orientation)
print(valid, page.pageSize().width(), page.pageSize().height(), page.orientation())

manager.addLayout(layout)

When I execute this code the layout is created and I get the following message in the console:

True 1189.0 841.0 1.
True means that the setPageSize was a success and the sizes are correct.

But when I open the layout in the interface I always get an A4 page.

Am I doing something wrong?


texts: How to convert geometry coordinates into Google lat/lons?

texts: For storing geospatial data, the geometry field is used in SQL Server or Oracle Spatial. I want to know is there any way where we can convert this geometry field into a set of latitudes and longitudes.


texts: Including two different maps in one QGIS Print Composer?

texts:
How do I use print composer in QGIS 2.0.1 to produce a final product with two different maps on the same sheet?

I can't figure out how to use the 'add a map' button to add anything other than a duplicate copy of what I already have.


texts: Raster processing "Tool is not licensed"

texts:
I have a piece of python code written that works fine on a small test area but when I blow up the size of the test area I get a "tool is not licensed" error.

arcpy.CheckOutExtension("Spatial")
spatial_reference = arcpy.SpatialReference(Environment.srid)
arcpy.env.outputCoordinateSystem = spatial_reference
arcpy.env.overwriteOutput = True
arcpy.env.snapRaster = Environment.dem
process_boundary = Environment.boundary

boundary = Environment.output_location + '/background/boundary'
f_ = 'FID'
arcpy.PolygonToRaster_conversion(process_boundary, f_, Environment.scratch_workspace + '/boundary', cellsize='1')
ras = Int(((Raster(Environment.scratch_workspace + '/boundary') + 1) / (Raster(Environment.scratch_workspace + '/boundary') + 1)) - 1)
ras.save(boundary)

The polygon to raster works fine and produces the expected result.
The scratch_workspace is a file geo-database.

Error:

File "G:\Code\models\background.py", line 30, in process
    ras = Int(((Raster(Environment.scratch_workspace + '/boundary') + 1) / (Raster(Environment.scratch_workspace + '/boundary') + 1)) - 1)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.3\ArcPy\arcpy\sa\Functions.py", line 4180, in Plus
    in_raster_or_constant2)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.3\ArcPy\arcpy\sa\Utils.py", line 53, in swapper
    result = wrapper(*args, **kwargs)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.3\ArcPy\arcpy\sa\Functions.py", line 4177, in Wrapper
    ["Plus", in_raster_or_constant1, in_raster_or_constant2])
RuntimeError: ERROR 000824: The tool is not licensed.

I have tried both Spatial and 3D for my defined license.


texts: Access QGIS API from Anaconda python environment using Conda-Forge

texts:
I am trying to use QGIS in stand-alone python scripts (windows 10, python 3.6.6).

After the following steps...

1) create new conda environment (geospatial)

2) install the qgis conda-forge package through conda install -c conda-forge qgis as per conda-forge page (https://anaconda.org/conda-forge/qgis)

3) run anaconda prompt and activate the geospatial environment

... I run python and try to import qgis.core. However I get the "ModuleNotFoundError" error.

(geospatial) C:\Users\vince>python
Python 3.6.6 | packaged by conda-forge | (default, Jul 26 2018, 11:48:23) [MSC v.1900 64 bit (AMD64)] on win32
Type "help", "copyright", "credits" or "license" for more information.
>>> import qgis.core
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ModuleNotFoundError: No module named 'qgis'

From reading other resources, it seems PYTHONPATH and PATH need to be set, but after setting those (see below), I still get an error on import qgis.core, this time an "ImportError: DLL load failed" error.

See the following anaconda prompt output:

(geospatial) C:\Users\vince>set PYTHONPATH=C:\Anaconda3\envs\geospatial\Library\python

(geospatial) C:\Users\vince>set PATH=C:\Anaconda3\envs\geospatial\Library\bin;%PATH%

(geospatial) C:\Users\vince>python
Python 3.6.6 | packaged by conda-forge | (default, Jul 26 2018, 11:48:23) [MSC v.1900 64 bit (AMD64)] on win32
Type "help", "copyright", "credits" or "license" for more information.
>>> import qgis.core
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "C:\Anaconda3\envs\geospatial\Library\python\qgis\__init__.py", line 72, in <module>
    from qgis.core import QgsFeature, QgsGeometry
  File "C:\Anaconda3\envs\geospatial\Library\python\qgis\core\__init__.py", line 34, in <module>
    from qgis._core import *
ImportError: DLL load failed: The specified module could not be found.

The (Python) paths I set above were "derived" from similar paths I found when one wants to use QGIS API by installing "osgeo4w. So, the main question: When using the Conda-Forge QGIS package, what are the required PYTHONPATHs and PATHs to be set; and/or what other configuration is needed to get the conda-forge QGIS package going?


texts: Calculate area of geometry derived from current feature in QGIS 2.8

texts:
I'm trying to create a new field defined as the area of the convex hull of the current feature, in the Field Calculator in QGIS 2.8.

Calculating the area of the current feature is easy ($area), and calculating the convex hull of the current feature is easy (convexHull($geometry)), but I can't find any way to calculate the area of that resulting geometry (as opposed to the geometry of the current feature).

I think I need to be able to either:

Apply an area calculation to that convex hull geometry in the Field Calculator -- this is what I'd prefer; or
Write a custom function -- there appears to be a free-standing area function in the Function Editor, but I can't figure out how to get a simple custom function defined and called properly; shouldn't need any arguments, as it's just the area of the convex hull of the current geometry; or
Create a new layer from an interim attribute that just contains the geometry of the convex hull -- then I could calculate the areas of those features in the new layer

Any suggestions on any of these approaches?


texts: arcpy.SpatialJoin_analysis giving ERROR 000210?

texts:
I have two polygon feature classes (buildings), one with attributes in green in the example below (building_sample.shp) and one with no attributes in red in the same example (building_sample_NA.shp).

I want the red buildings to take the attribute values of their closest green neighbor using arcpy.

Following the description of the SpatialJoin_analysis function, that should be easy but my code is not working:

  ERROR 000210: Cannot create output

import arcpy
arcpy.env.workspace = "C:/DemoData/outputs"
arcpy.env.OverwriteOutput = True

arcpy.SpatialJoin_analysis (target_features = "building_sample_NA.shp",
                      join_features = "building_sample.shp",
                      out_feature_class = "building_completed",
                      join_operation = "JOIN_ONE_TO_ONE",
                      match_option = "CLOSEST")


texts: Difference between Geographic-3D and Geocentric CRS

texts:
I'll start with the OGP definitions:

Geographic. A coordinate reference system based on a geodetic datum and using an ellipsoidal (including spherical) model of the Earth. This provides an accurate representation of the geometry of geographic features for a large portion of the Earth’s surface. Geographic coordinate reference systems can be two- or three-dimensional. A Geographic 2D CRS is used when positions of features are described on the surface of the ellipsoid through latitude and longitude coordinates; a Geographic 3D CRS is used when positions are described on, above or below the ellipsoid and includes height above the ellipsoid. [...]
Geocentric. A coordinate reference system based on a geodetic datum that deals with the Earth’s curvature by taking the 3D spatial view, which obviates the need to model the curvature. The origin of a geocentric CRS is at the centre of mass of the Earth.

In practice, I would identify the difference between, say, a Geographic-3D and a Geocentric CRS in the fact that the latter has a central point in the 3D space from which angles and elevation are computed (no ellipsoid is needed, Earth is spherical?), whereas the former refers to the datum/ellipsoid, i.e. latitudes, longitudes and heights are computed with respect to the ellipsoid.

Am I saying something wrong or missing something?
Thanks.


texts: Example of simple kriging with varying local means in R

texts:
I'm looking for an example of how to use simple kriging with varying local means in R.

I've tried google, a few books (for example Aplied Spatial Data Analysis with R) and I've found almost nothing. The only useful thing I've found is the Figure 2.1 on page 13 of Gstat user's manual (http://www.gstat.org/gstat.pdf).

Could someone provide an example of simple kriging with varying local means in R?


texts: Can I add OpenCycleMap to ArcGIS Explorer Desktop?

texts:
I'm using the ArcGIS Explorer Desktop (free viewer) - build 1750 (latest).
It comes with several built in base maps, including OpenStreetMap.
Is it possible to add to it OpenCycleMap:
http://www.opencyclemap.org/
as well ?


texts:
Basically I created a MBtiles using Tilemill and it generates a map with six zooms.
I implement it into tilestream.
I added a layer on the openlayers module(Drupal) in XYZ:

site:8888/v2/geography-class/${z}/${x}/${y}.png

The problem is this URL sent by openlayers: site:8888/v2/geography-class/-1/x/x.png

Each time, whatever the zoom I am on the map, the zoom value always stays at -1, the others values are okay.

If I try to reach site:8888/v2/geography-class/1/1/1.png it works.

What did I miss?


texts: Transfer & maintain a relational database in QGIS - backend vs frontend

texts:
I'm going to set up a small db for my own use in QGIS. I can see only a few of the many tables will be spatial, but the ease of use and capacity of the widgets in QGIS are attractive. I won't writing (much) sql on the command line.

How will I;
- transfer the relationships / joins between computers
- share the data and with anyone else (can the joins only be reproduced via a QGIS project)
- use the relationships for this dataset, as a subproject of other analysis

I understand that tables (spatial or not) can exist in postgresql/postgis but important rdb function like foreign keys are not understood by QGIS. Instead joins, default attribute values, some attribute typing, need to be defined in the QGIS widgets. These constraints then become ephemeral with the project file.

Saving the widgets as part of each layer style is I presume possible, but will not provide a robust solution I think because there will be a lot of intermediate tables.

As far as I know QGIS projects don't save to pgsql.

This earlier question explores some of the issues QGIS managing relationships and these posts from Opengis.ch are helpful too. https://www.opengis.ch/author/signedav/ But what workflows are out there?

My use case is for a bunch of data that exists almost entirely in a spreadsheet, it is only once column that needs to be exposed to a spatial join. However to produce a report back the other way, i.e. select a point on the map and reveal the joined attributes, I believe the rdb needs to live in QGIS.    I have spent some time setting up in DBeaver and pgsql, but foreign key joins don't get recognised in QGIS.

Any thoughts appreciated


texts: I am working with elevation data layer. I need to change a few of the pixel values of a raster layer manually. I use ArcGIS, ERDAS, QGIS. Few of the pixel values were deviated widely than its surrounding pixel values as shown in the image. I need to change them to look better.


texts: Is it possible in uDig to use Features from map to create a profile in profile View?

texts:
Its great to have the Profile View in uDig. The user can create a profile withe the profile Tool. Its required to create a line in the Map

I was wondering if there is a possibility to select a feature in the map and use the geometry to calculate the profile

Is there a chance to to so, How can I solve this requirement? is it possible to switch measurement unit and use the same is used for the distance tool?


texts:
I regularly receive files from a service authority in DWF format, which shows the location of their assets in there area.  I note that there is no option to have the data in an alternate format as it is automatically generated.

I can't seem to convert this file into an alternate format for import into QGIS (for example DXF).  I note that it can be opened in Autodesk Design Review but this cannot convert the data (other than printing to a PDF).

How do others incorporate this data format into their work?


texts: Cutting rasters using polygon and obtaining pixel specific data using ArcMap

texts:
I have a georeferenced polygon file that I want to use to cut several different raster files. I am looking to analyse the individual pixels within those cut areas.

One of the raster files has a spatial resolution of 1 arc-second per pixel (raster A). For each of those pixels, I want to obtain a corresponding value from the other raster files (rasters B-E). The other raster files have different spatial resolutions, mostly three arc-seconds.

The need to output as an Excel format where each row represents a pixel and where columns contain values from rasters A-E.

How can this be done using ArcMap?


texts: How to initialize the QGIS from standalone script?

texts:
I tried to initialize QGIS  from a standalone script, but it's giving the error (find in an image).

I reinstalled, but the same issue occurred.
code:
from qgis.core import *

QgsApplication.setPrefixPath(r"C:\Program Files\QGIS 3.12", True)
qgs = QgsApplication([], False)
qgs.initQgis()
vlayer = QgsVectorLayer(path_to_airports_layer, "Airports layer", "ogr")
QgsProject.instance().addMapLayer(vlayer)
qgs.exitQgis()

How can I initialize correctly?


texts:
I think that must be a "simple" way to import several Geotiff files into R instead of copy and paste n times Rfd<-readGDAL() and then save them as individual objects.

I found the following code in a post regarding how to import multiple shapefiles ...

get all files with the .shp extension from working directory
setwd("D:/GIS_DataBase/GIS_Tirol/Tirol_Verbreitungskarten/Verbreitungs_Daten")
shps <- dir(getwd(), ".shp")
the assign function will take the string representing shp and turn it into a variable
which holds the spatial points data
for (shp in shps) assign(shp, readShapePoints(shp))
plot(get(shp[1])) # i.e.
...done*

But it does not translate to raster files with readGDAL {rgdal} function.

Other post suggested the use of for or/and function, when importing several *.csv files into R. However, I have been equally unsuccessful into make it work.


texts: How to create MXD in folders programmatically,

texts:
I have an MXD as a template and about 650 folders.  Each folder corresponds to a species.

What I want to do is to make an MXD (same as the template) in every folder of each species.

How to do this programmatically?


texts: How to display special German characters "ÄÜÖß" in a map?

texts:
I have a shapefile with names for Regions. These don't seem to contain special characters "Ä", "Ü", "Ö", "ß" but ugly hieroglyphs.

What can I do to display them correctly. If I import the shape as UTF-8 it doesn't work.

The names are in a dbase table. So is it possible to save a dbase table in UTF-8 or is there a different way of dealing with it?


texts: Simplify a shapefile/GeoJSON polygon without "uncovering" territory?

texts:
I've been using Mapshaper to simplify some very complex polygons (boundaries of countries). The goal is to create relatively simple polygons which can be used to determine whether a given lat/lon point is inside a given country. It's OK if I (somewhat rarely) wind up with a false positive, but it's not OK to ever wind up with a false negative.

Using the standard Mapshaper simplification tools, I wind up with a polygon that isn't always covering 100% of the territory covered by the original polygon. That's a problem. Is there any tool/technique that exists to create a simplified polygon that will fully enclose all the territory the original polygon covered?

Update:

To clarify: While I'll be processing one country at a time, I will be querying against a collection of countries to find which of them contains a given point. My intention is that the processed polygons will sometimes overlap, such that points right by the border between two countries will often fall in both of those countries according to the polygons. This is OK, and expected.


texts: How to successfully deploy a web map application built with Web AppBuilder for ArcGIS on local server?

texts:
Following some research, I have built a web map application using web app builder for ArcGIS using my ArcGIS Onlince (AGOL) account which runs properly. With aims to extend some functionalities and themes of the application, I downloaded the app and have inserted in my Internet Information Services (IIS) - C:\inetpub\wwwroot of my pc.

when I run the .index file, it does not yield any results. It shows loading sign but never fully loads.

When I use the link; "http://"machinename"/"name of application"/ it equally does not work. It start to load and later shows "this page cannot load" check your proxy settings. So far from my findings it appears it is proxy settings which is the issue but I am unable to resolve it at this level.

Is there another simple or good means by which I can deploy my application on my PC? What are the steps or where is a link to where I can research this further?


texts:
I currently have a table (not a text box) that I'd like to have overlaid on the map in layout view of ArcMap 10.3. It, unfortunately, has a long title in the column. So, when I go to put the table on the map, it looks very elongated and generally not great. Is there a way that I can text-wrap the titles, and potentially also center-justify the data within the table? Neither the Options menu under Customize nor the Properties menu appear to have any helpful fixes.

I guess if there really isn't any other option, I could just make a text box and then align everything by hand, but obviously this isn't preferred.


texts: Having trouble replaceDataSource loop in Python

texts:
In ArcGIS (using 10.4) I have my layout set up, shapefile A symbolized and labeled. What I need to do is print this layout to pdf, then change the source of shapefile A to the next in a
series of several hundred, have it look the same, change the title, and print the pdf. I'm trying to loop through each shapefile in a folder. The shapefiles are polygon. The legend stays the same because each shapefile has the same number of
polygons, symbolized the same way. The labels should stay the same, but move as the polygons they label change shape and position. The extent is constant.
This is easy to do interactively in ArcGIS: go into the Layer Properties, Source tab, Set Data Source...
In trying to do this in Python, I'm having difficulty. My code maybe isn't looping. It prints one map. I think the replaceDataSource {dataset_name} needs to be a variable (the ""), but I can't figure how to do that.

>>> import arcpy
... import os
... mxd = arcpy.mapping.MapDocument("CURRENT")
... df = arcpy.mapping.ListDataFrames(mxd, "*")[0]
... lyr = arcpy.mapping.ListLayers(mxd, "*", df)[3]
... output_dir = r"\\path\PDFs"
... arcpy.env.workspace = r"\\path\Set1Test"
... shapelist = arcpy.ListFeatureClasses()
... for shape in shapelist:
...     lyr.replaceDataSource(r"\\path\Set1Test", "SHAPEFILE_WORKSPACE", "")
...     TextElement = arcpy.mapping.ListLayoutElements(mxd, "TEXT_ELEMENT", "MapTitleText")[0]
...     TextElement.text = lyr.datasetName
...     arcpy.mapping.ExportToPDF(mxd, r"\\path\PDFs\Map " + lyr.datasetName)

I've also tried this, with the same result:

>>> import arcpy
... import os
... mxd = arcpy.mapping.MapDocument("CURRENT")
... df = arcpy.mapping.ListDataFrames(mxd, "*")[0]
... lyr = arcpy.mapping.ListLayers(mxd, "*", df)[3]
... output_dir = r"\\path\PDFs"
... for file in os.listdir(r"\\path\Set1\Set1Test"):
...     if file.startswith ("Set"):
...             lyr.replaceDataSource(r"\\path\Set1\Set1Test", "SHAPEFILE_WORKSPACE", "")
... TextElement = arcpy.mapping.ListLayoutElements(mxd, "TEXT_ELEMENT", "MapTitleText")[0]
... TextElement.text = lyr.datasetName
... arcpy.mapping.ExportToPDF(mxd, r"\\path\PDFs\Map " + lyr.datasetName)


texts: How to buffer a vector shapefile using ogr python?

texts:
I am trying to learn how to use ogr in python using the country and populated places datasets from http://www.naturalearthdata.com/downloads/50m-cultural-vectors/. I am attempting to use filters and buffers to find points (ne_50m_populated_places.shp) within a specified buffer of a named country (filtered from feature class ADMIN in ne_50m_admin_0_countries.shp). The problem appears to be that I don't understand what units to use for buffer(). In the script I have simply used an arbitrary value of 10 to test if the script works. The script runs but returns populated places from around the Carribean region for named country 'Angola'. Ideally, I want to be able to specify a buffer distance, say 500km, but can't work out how to do this as my understanding is buffer() is using the units of countries.shp that will be in wgs84 lat/long format. Advice on method to achieve this would be much appreciated.

# import modules
import ogr, os, sys

## data source
os.chdir('C:/data/naturalearth/50m_cultural')

# get the shapefile driver
driver = ogr.GetDriverByName('ESRI Shapefile')

# open ne_50m_admin_0_countries.shp and get the layer
admin = driver.Open('ne_50m_admin_0_countries.shp')
if admin is None:
  print 'Could not open ne_50m_admin_0_countries.shp'
  sys.exit(1)
adminLayer = admin.GetLayer()

# open ne_50m_populated_places.shp and get the layer
pop = driver.Open('ne_50m_populated_places.shp')
if pop is None:
  print 'could not open ne_50m_populated_places.shp'
  sys.exit(1)
popLayer = pop.GetLayer()

# use an attribute filter to restrict ne_50m_admin_0_countries.shp to "Angola"
adminLayer.SetAttributeFilter("ADMIN = ANGOLA")

# get the Angola geometry and buffer it by 10 units
adminFeature = adminLayer.GetFeature(0)
adminGeom = adminFeature.GetGeometryRef()
bufferGeom = adminGeom.Buffer(10)

# use bufferGeom as a spatial filter on ne_50m_populated_places.shp to get all places
# within 10 units of Angola
popLayer.SetSpatialFilter(bufferGeom)

# loop through the remaining features in ne_50m_populated_places.shp and print their
# id values
popFeature = popLayer.GetNextFeature()
while popFeature:
  print popFeature.GetField('NAME')
  popFeature.Destroy()
  popFeature = popLayer.GetNextFeature()

# close the shapefiles
admin.Destroy()
pop.Destroy()


texts: Deleting multiple features (selected by expression) at once using PyQGIS?

texts:
I'm trying to select multiple features at once to do batch editing. The first selection is supposed to be the expression within the request, but I can't seem to get anything to select. I want the selection to be deleted.

##layer=vector
##zone_index=field layer
##subty_index=field layer
##sfha_index=field layer

from qgis.core import *import qgis.utils
from PyQt4.QtCore import *

layer=processing.getObject(layer)

layer.startEditing()

request = QgsExpression(u"\"zone_index\"='AREA NOT INCLUDED' OR \"zone_index\"='OPEN WATER' OR (\"sfha_index\" = 'F' AND \"zone_index\" = 'D')")

it= layer.getFeatures( QgsFeatureRequest( request))

ids = [i.id() for i in it]
print ids

layer.setSelectedFeatures(ids)

layer.commitChanges()


texts: Breaking single polygon into multiple polygons?

texts:
I created voronoi polygons after which I dissolved them on the basis of one attribute. Now there are cases where in for one single attribute "abc" I have3 polygons at different places with same object id. I would like to break these into 3 polygons as they are at different locations.

Is there any tool in ArcGIS Desktop or QGIS which can help me in doing that?


texts: How to give a table some geometry from a view on QGIS using a join?

texts:
So basically I've got this static table that people update information

Then I need it to give it some geometry from a view table. This geometry is dynamic, that's why  can't just make a new geom column on the static table

Both layers have matching ids, but different info. I'd like to give the geometry from the view to the table, but when doing a join on QGIS, it doesn't offer the option of adding the geometry column.

Is there a way to join the geometry?


texts: PostGIS St_AsMVT latitude incorrect at low zoom levels

texts:
I'm using the St_AsMVT feature in Post GIS to create mapbox tiles.

This works ok at high zoom levels (>9) but at lower levels the latitude appears to be significantly off...which gets worse as the zoom level decreases.

Has anyone else experienced this? Is it a bug in PostGIS? In the GIF the red box should stay in it's original position, but moves north in lower zoom level tiles.


texts: How can I use Planet Labs tiles with CartoDB?

texts: I'd like to use Planet Labs mosaics as a base layer for my CartoDB map. But it seems to require authentication of some sort. What's the right URL to add it as a custom basemap to my CartoDB map?