texts:
Given that Folium builds on the Leaflet JS library, I was wondering if it could provide the same functionality as Leaflet for custom icon markers.

Scanning through the tutorial here, I'm somewhat doubtful that this is possible, but I thought I would ask here to confirm. It looks as if the marker icons are limited to those provided on Bootstrap.

map_1 = folium.Map(location=[45.372, -121.6972], zoom_start=12,tiles='Stamen Terrain')
map_1.simple_marker([45.3300, -121.6823], popup='Some Other Location',marker_color='red',marker_icon='info-sign')
map_1.create_map(path='iconTest.html')

If this is not supported, is anyone familiar with another Python mapping framework that would support the use of markers with custom icons?


texts: Is there a name for deformation-based maps?

texts:
I have already seen maps using deformation to visualize information. For example, this map shows temporal distance between some French cities when traveling by train: http://www.nf-ameublement-pro.com/sc/images/new/newsletter/septembre/carte_france.jpg

Is there a name for this kind of map? Also, is this type of map available in any GIS software?


texts: QGIS - Feature vs Layer fields out of sync. Is this by design?

texts:
I'm scripting in QGIS 3.6 using python. Looking at the API and trying out attribution-editing functions, I found that a vector layer has a QgsVectorLayer::fields() set of fields, and the features inside have a set of QgsFeature::fields().

If I use addAttribute() on the layer (in between the startEditing and commitChanges guards), then immediately query a layer-feature's fields, I find that the fields are out of sync (new field is in the layer's fields, but not in the feature's fields) and calling feature.setAttribute() with the new string fails with a  key error.

I'm wondering:
1. if there is a practical need to have these things out of sync (ie: is it ever desirable / useful)
2. what normal people do to keep layer & feature fields in sync, if it is undesirable.  Is there a simple call where I can flush the layer's fields into all the features fields quickly?


texts: Joining geometry tables in GeoPackage using SQL

texts:
I want to join two geometry tables (table1 containing and table2) inside the same GeoPackage and eventually two tables in two separate GeoPackage. Just to be clear, when I open the GeoPackage in QGIS, I only have table1 which is in my case several lines that have multiples attributes. I want the result to be stored as an updated table1 with more columns. I'm also interested in knowing how to create a new table in the first GeoPackage. Ideally, I wish to join only one column (column1) of table2 and keep all columns of table1. In both table, the common column is named fid.
import sqlite3
conn = sqlite3.connect("C:/temp/bd.gpkg")
cursor = conn.cursor()
req = " ".join([
    "CREATE TABLE table_new AS",
    "SELECT table1.*, table2.column1",
    "FROM table1",
    "LEFT JOIN table2",
    "ON table1.fid = table2.fid"
])

cursor.execute(req)
conn.commit()
cursor.close()
conn.close()

It works but doing so create a table without geometry at all. I tried without the CREATE TABLE line but then nothing is modified. I'm new to SQL and Geopakage so I suppose I'm missing something pretty basic.


texts:
I'd like to define the default value for a field in QGIS using OpenWeather API request.

We use QField for outdoor data collecting, and some of the actual weather conditions are important for us. I've an OpenWeather API key, and I want to insert the temperature, the visibility, the wind speed and direction to the field as human-readable text.

The API request is something similar:

http://api.openweathermap.org/data/2.5/weather?lat=47.6&lon=20.4&appid=my_api_key&units=metric.

In the browser the result is:
{"coord":{"lon":20.4,"lat":47.6},"weather":[{"id":800,"main":"Clear","description":"clear sky","icon":"01n"}],"base":"stations","main":{"temp":0,"pressure":1037,"humidity":73,"temp_min":0,"temp_max":0},"visibility":10000,"wind":{"speed":1.32,"deg":219.501},"clouds":{"all":0},"dt":1551050100,"sys":{"type":1,"id":6676,"message":0.0032,"country":"HU","sunrise":1550986040,"sunset":1551024985},"id":716877,"name":"Nyulasszőllő","cod":200}

I'd like to write this data to the field something like this:
0°C; visibility: 10000 m; wind: 1.32 m/s, 219.501°;

Is it possible to fill my field with this method in QGIS?
Will it work with QField?
How should I define the default value for the form?


texts:
I am trying to access and display some data in a GeoJSON layer.

Can someone take a look at my code below and suggest what I need to do to make it work?

The result of the console.log is:

  Object {latlng: o.LatLng, layerPoint: o.Point, containerPoint:
  o.Point, originalEvent: MouseEvent, type: "click"…}

var wrldbndryData;

function wrldbndryStyle(feature)
{
    return {
        weight: 2,
        opacity: 0.1,
        color: 'black',
        fillOpacity: 0.1
    };
}

function highlightFeature(e)
{
    var layer = e.target;

    layer.setStyle(
    {
        weight: 9,
        color: '#666',
        fillOpacity: 0.5
    });

    if (!L.Browser.ie && !L.Browser.opera)
    {
        layer.bringToFront();
    }
}

function resetHighlight(e)
{
    wrldbndryData.resetStyle(e.target);
}

function showFeature(feature)  <-- *** THIS IS WHERE I AM TRYING TO ACCESS THE GEOJSON DATA *****
{
    console.log(feature);
    if (feature.properties && feature.properties.NAME)
    {
        $('#sidecol').html(' COUNTRY INFO </center' + '
' + feature.properties.NAME + '');
    }
}

function eachFeature(feature, layer)
{
    layer.on(
    {
        click: showFeature,
        mouseover: highlightFeature,
        mouseout: resetHighlight
    });
}

wrldbndryData = L.geoJson(wrldbndry,
{
    style:wrldbndryStyle,
    onEachFeature:eachFeature
});

wrldbndryData.addTo(map);


texts:
I've collected polygon data using ODK's geo_shape. I'm trying to convert it to a geojson for use in QGIS. I'm using python 3, pandas and geojson to do this.

The geo_shape output from ODK looks like this:

51.781223 -1.2891759999999999 0.0 0.0;51.78472...

Coordinates for each vertex contain two extra values (0.00 0.00) that are not needed (z value and accuracy I believe) and the x and y are the wrong way around. Therefore I've put together an ugly method of reshaping each feature's coordinates to something more useable. The processed geo_shape coordinates look like this:

[[-1.2891759999999999, 51.781223], [-1.266689, 51.784727], [-1.2407679999999999, 51.782497]]

and are added to a new column in the dataframe, 'myGeom':

result = [reshape_geom(x) for x in df_poly['habitat_area']]
df_poly['myGeom'] = result

Using the processed coordinates and this question I've got very close.

However, the resulting geojson will not draw in QGIS3 unless the coordinate list is enclosed in 4 square brackets [[[[my_x_1, my_y_1],[my_x_2, my_y_2],[my_x_3,my_x_3]]]].

The code looks like this:

features = df_nona.apply(
lambda row: Feature(geometry=MultiPolygon(row['myGeom'])),axis=1).tolist()

properties = df_nona.drop(['columns to drop'], axis=1).to_dict('records')
feature_collection = FeatureCollection(features=features, properties=properties)

with open('file06.geojson', 'w', encoding='utf-8') as f:
    json.dump(feature_collection, f, ensure_ascii=False)


texts: Same size polygons giving different overall areas

texts:
I am currently working with multiple shape files that I have clipped to the exact same extent, but when I calculate area (acres) within each polygon.
I am getting different results between the different shape files.
Each file represents a different variable, ie: I am working with a slope polygon, aspect, vegetation, etc.
The acreage results (from calculate geometry) are similar, but not the same (9015, 10294, 8999).
There are multiple geometries, each different for each shapefile, and I am comparing the acreage sums across each shapefile that I have clipped to the same size.
The plotted clipped extents are mostly just squares with some irregular shaped polygons.

Fixing the topology errors did not make the acreage's equal.
They are all on the same projection.
Fixing geometry errors did not work either.


texts: Overriding automatic projection definitions for layers in QGIS

texts:
Is there a way to set an override for a specific projection in QGIS?

I want layers that have a specific projection to be set with a custom projection I've defined.

I can define it manually, but wondering if there's a way to do it automatically for all current and new layers in the project.


texts:
This seems like it should be very easy to do, but I am unable to click on a point and have the map recognize that I have clicked on the point. I first tried just creating a simple point, setting the query out fields to "*" and the spatial relationship to "INTERSECT". That didn't work so now I've been trying to create a polygon around where I click and setting the spatial relationship to both "WITHIN" and "CONTAINS" (two separate tries), I've also tried using an envelope instead of a polygon to no avail.

I have tried using  this  example directly and it works with the given map because the features are polygons, however when I change the URL's to my feature server (which is all points) it no longer does anything when I click, the length of the results array is zero.

Here is the code for my current attempt using an envelope (I left my comments in where I tried making a polygon as well) inside my OnSingleTapListener:

clicked = mMapView.toMapPoint(v1, v2);
int dimens = mMapView.getExtent().getDimension();
double pixW = dimens / mMapView.getWidth();
double tolerance = 10 * pixW;
Envelope env = new Envelope(clicked, tolerance, tolerance);
Query query = new Query();
query.setGeometry(env);
query.setOutFields(new String[] {"*"});
//Polygon poly = new Polygon();
//double x = clicked.getX();
//double y = clicked.getY();
//poly.startPath(x-2.0, y+2.0);
//poly.lineTo(x+2.0, y+2.0);
//poly.lineTo(x+2.0, y-2.0);
//poly.lineTo(x-2.0, y-2.0);
//poly.lineTo(x-2.0, y+2.0);
query.setSpatialRelationship(SpatialRelationship.WITHIN);
query.setInSpatialReference(mMapView.getSpatialReference());
featLayers[0].selectFeatures(query, ArcGISFeatureLayer.SELECTION_METHOD.NEW,

Does anyone have any idea why this would not work when the exact same method works for other features that aren't points? The example also has me using an ArcGISDynamicMapServerLayer, is this necessary? I have tried it both with and without this layer although I prefer it without because the layer loads slowly and makes the map look worse.


texts: How to get GPS coordinates of neighborhood polygon from Google Maps?

texts:
Does anyone know how I can find the GPS coordinates for neighborhoods like this:

https://www.google.com/maps/place/SoHo,+New+York,+NY/@40.7236448,-74.0006793,16z/data=!3m1!4b1!4m2!3m1!1s0x89c2598c5e51e2f3:0xce2f3449f490f818

Not just in New York, but ideally anywhere else in the US.


texts: Distribute Hotel Beds According to Population Density in Buffers using ArcGIS for Desktop?

texts:
I am running an ArcGIS for Desktop analysis where I am plotting competing hotel facilities and then distributing their beds in a 5-mile buffer around the site.

I would like to distribute the beds such that they are distributed according to where the population resides within that 5-mile radius, as opposed to evenly across the buffer.

That way, when I run an intersect - if the buffers overlap on top of a densely populated area, it will pick up more beds than if the same area overlapped on top of a sparsely populated area.

Any ideas?


texts:
I want to define no data values when working with gdal.BuildVRT(). According to the gdal.documentation http://gdal.org/python/osgeo.gdal-module.html, by creating a gdal.BuildVRTOptions (srcNodata=-99, VRTNodata=-99), I should be able force the gdal.BuiltVRT() to incorporate nodata values.

Yet, doing so as one can be seen below, does not work.

gdal.BuildVRT(file_open,None,srs_wkt,gdal.GRA_Bilinear,gdal.BuildVRTOptions (srcNodata=-99, VRTNodata=-99))

and I get an error message:

Traceback (most recent call last):
  File "/home/esaiet/Dropbox/Python_programing/Python_scripts/converting_porjection_vrt_II.py", line 45, in <module>
    vrt_ds=gdal.AutoCreateWarpedVRT(file_open,None,srs_wkt,gdal.GRA_Bilinear,gdal.BuildVRTOptions(srcNodata=None))
  File "/usr/local/lib/python3.6/dist-packages/GDAL-2.2.1-py3.6-linux-x86_64.egg/osgeo/gdal.py", line 2872, in AutoCreateWarpedVRT
    return _gdal.AutoCreateWarpedVRT(*args)
TypeError: in method 'AutoCreateWarpedVRT', argument 5 of type 'double'

What am I doing wrong?


texts: Error message when attempting to convert Landuse polygon to Raster in QGIS 2.18.14

texts:
I may be doing something obviously wrong, but I'm trying to convert a landuse polygon to a raster file for later use in QGIS 2.18.14. I've used the Rasterize function within QGIS succesfully previously but for this file I'm getting the following error message

'ERROR 6: One and only one of -3d, -burn or -a is required.'

Any ideas as to what may be going wrong?


texts:
I'm attempting to create ideal motorcycle routes (a tourism map) using network analyst, I have an attribute for road curvature. So far I've been able to do two things: create routes using either popular destinations (hotels, restaurants, etc...) or the midpoints to road segments with low sinuosity as stops. The problem which is clearly illustrated in the map below, is that network analyst doesn't create loops. It will go to a destination and then back tract away from it. Ideally if someone is merely going for a "joyride" they would be able to continue on their route with option of forking off to destinations along the route.
The only solution I've found is to preserve the first and last destination and a new destination on top of first one. Then after solving the route again, it will loop back to the first point. But I would like to improve it so that their are several loops. I'm hoping to create enough so that riders could find a loop around their area, but still be able to get to other loops.
Any additional ideas?

I've added an idealized map above, the green is a closed loop and the red lines are off-shoots to destinations.


texts: overpass-api or osmconvert: which method is more appropiate to use for files from 10 to 400 mb?

texts:
new to openstreetmap and this technology.

update: some mentioned postgis: well - The database schema of postgis allready contains the following tables: nodes, ways, relations, node_tags, way_tags, relation_tags, way_nodes.... well - can i do imports into postgis - as easy as into mysql?

i want to gather data from .o5m & pbf file - several ways are possible:

what is wanted: i look for restaurants - with all of the following tags:  [code]@id @lon @lat amenity name adress [as town, street, housnumber] website and mailadress. [/code]

All interests are based on the regions of Germany see http://download.geofabrik.de/europe/germany.html i have several files i want to work with - ranging form 10 MB (Bremen) to 390 MB (Nordrhein Westfalen) the osm.pbf-files are not too big;

main question: which method is the best and the most appropiate? - to store the  results in a mysql-db or just have big calc-sheets (with csv-data)

To achieve this goal of having only certain features in my csv file, we will have to filter the csv file generated by osmconvert afterwards.
Also, the arguments to the --csv option of osmconvert are the tag keys we want to have in our csv file.
If we want "amenity" and "name" columns in addition to id, lat and lon, we specify "@id @lon @lat amenity name", and then filter the CSV output.

a. use osmconvert to convert osm to .o5m and .pbf file it works fine but am having problem converting osm to .csv or excel file

example:

I can use

 ./osmconvert addis.osm --all-to-nodes --csv="@id @lon @lat amenity name" --out-csv -o=outfile.csv

and subsequently:

grep cafe outfile.csv > cafes.csv
grep restaurant outfile.csv > restaurants.csv
cat cafes.csv restaurants.csv > cafes_and_restaurants.csv

b. on the other hand i can do this way: use overpass-api:

<query type="node">
  <has-kv k="place" v="city"/>
  <has-kv k="name" v="any town "/>
</query>
<query type="node">
  <around radius="1000"/>
    <has-kv k="shop"/>
</query>
<print/>

question: are there more (effficient) ways!? i need to have the most efficient.

b. - the overpass-api-method creates a XML-file (more or less large )
for example a XML file which I created from the overpass api. This i would need to load into a mysql database.
btw: i use this to LOAD XML LOCAL INFILE '/home/jay/Downloads/interpreter1' into table jayDB.xml1 (id);

I'm new to XML so not sure if this is possible.

<osm version="0.6" generator="Overpass API">
<note>
The data included in this document is from www.openstreetmap.org. The data is made available under ODbL.
</note>
<meta osm_base="2013-03-07T14:54:02Z" areas="2013-03-07T10:37:02Z"/>
<node id="240486180" lat="50.9744274" lon="3.0152858">
<tag k="addr:housenumber" v="9"/>
<tag k="addr:street" v="Marktplaats"/>
<tag k="amenity" v="cafe"/>
<tag k="email" v="vandaelekoen67@skynet.be"/>
<tag k="name" v="Paviljoentje"/>
<tag k="opening_hours" v="Mo-Su 09:00+; Tu off; Th 09:00-14:00"/>
<tag k="phone" v="+3251636211"/>
<tag k="website" v="http://www.paviljoentjestaden.be"/>
</node>
<node id="244312208" lat="51.2461401" lon="5.4390455">
<tag k="amenity" v="cafe"/>
<tag k="created_by" v="JOSM"/>
<tag k="name" v="De Club"/>
</node>

however i think due to the format of the file it thinks K and V are the columns. Where as I was amenity to be column and cafe to be the value
i tried to look atthe details https://dev.mysql.com/doc/refman/5.5/en/load-xml.html

as mentioned above: the main question: which method is the best and the most appropiate? - to store the  results in a mysql-db or just have big calc-sheets (with csv-data)

love to hear from you


texts: Uncertainty of a line between two coordinates with uncertainties in QGIS

texts:
I am trying to plot the uncertainty of a line between two coordinates with uncertainties. I have created polygons around 2 coordinates to represent the uncertainty region around each point.  I then create random points in each polygon, and then join them using the points to paths function.  Looks something like this:

I want to plot a region that represents the probability of the lines, and also takes the uncertainty region of the original endpoints into account - should look like this;

The uncertainty region needs to be at the 95% confidence level.  I am new to QGIS, I thought maybe by creating a grid and populating cells within the grid I could somehow come up with a polygon that represents a probability density area?
I am using QGIS Desktop 3.0.1 with GRASS 7.4.0. I believe I will need large numbers of random points/lines to achieve what I want which may make processing any algorithms slow?


texts:
i have field named (depth) in feature class (point) with this spec :

float - 5 digits  >>>  0.7186   or 0.9695

i create new field named (depth_new) with field calculator : round([depth],1)
the result will be >>>  0.7186 = 0.7  or 0.9695 = 1

then i want to use (depth_new) field as label in my point feature class using arcpy, here is some of my statement in my code :

      lyr.labelClasses[0].expression = '"{}" & (left([depth_new],1)) & " " & "<SUB>" & (right( [depth_new],1)) & "</SUB>" & "{}"'.format("<FNT size = '17'><CLR blue = '255'>","</CLR></FNT>")

it is works fine,
but as you can see on my code, i use "right" and "left" for the digit of the number that i want to show as label,

the problem is lets take a look at record with value = 0.7 it will be come 0.7 , but with value = 1 , the label will be 1.1 (so the number is changed, the correct should be 1.0)

i know i can solve this in field properties and check "pad with zeroes" and use 2 digit in properties, but let say i have 100 layers, so should i really do that for every 100 layers? i use arcpy for labelling to skip repitition and save time,

so any workaround or alternative way for this?


texts: List Attribute values selected by location

texts:
I would like to get a list of road features that intersect every county in a state and output it as CSV or Excel.
Ex: county A intersect the Road Names

a
b
c

county B intersect the Road Names

e
f

Any ideas using FME?


texts:
I've updated the code with the suggestions in the code but still not working.

arcpy.env.workspace = r'F:\EGM722 - Customising GIS apps\Assignment\folklore.gdb'
arcpy.env.overwriteOutput = True
arcpy.Delete_management('folklore_shape.shp')
#arcpy.Delete_management(r'F:\EGM722 - Customising GIS apps\Assignment\folklore_shape.shp')
arcpy.MakeXYEventLayer_management('Folklore', 'x_coord', 'y_coord', 'folklore_Layer', "PROJCS['IRENET95_Irish_Transverse_Mercator',GEOGCS['GCS_IRENET95',DATUM['D_IRENET95',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Transverse_Mercator'],PARAMETER['False_Easting',600000.0],PARAMETER['False_Northing',750000.0],PARAMETER['Central_Meridian',-8.0],PARAMETER['Scale_Factor',0.99982],PARAMETER['Latitude_Of_Origin',53.5],UNIT['Meter',1.0]];-5022200 -15179500 10000;-100000 10000;-100000 10000;0.001;0.001;0.001;IsHighPrecision", '#')
arcpy.FeatureClassToFeatureClass_conversion("folklore_Layer", "F:\EGM722 - Customising GIS apps\Assignment", "Folkore_shape.shp")

I've tried arcpy.Delete_management('folklore_shape.shp') and arcpy.Delete_management(r'F:\EGM722 - Customising GIS apps\Assignment\folklore_shape.shp')

Getting error;

Exception in Tkinter callback
Traceback (most recent call last):
  File "C:\Python27\ArcGIS10.5\lib\lib-tk\Tkinter.py", line 1542, in __call__
    return self.func(*args)
  File "F:\EGM722 - Customising GIS apps\Assignment\folklore - 22-04-18 - diff_folder.py", line 132, in new_record
    arcpy.FeatureClassToFeatureClass_conversion("folklore_Layer", "F:\EGM722 - Customising GIS apps\Assignment", "Folkore_shape.shp")
  File "C:\Program Files (x86)\ArcGIS\Desktop10.5\ArcPy\arcpy\conversion.py", line 1891, in FeatureClassToFeatureClass
    raise e
ExecuteError: ERROR 000258: Output F:\EGM722 - Customising GIS apps\Assignment\Folkore_shape.shp already exists
Failed to execute (FeatureClassToFeatureClass).

Creating the event layer works it is the delete part that isn't working.


texts:
I moved from QGIS 2 to QGIS 3 I cannot seem to be able to use the zonal statistics package.
Basically what I want to do is take the world pop GeoTIFF for a country and then calculate the population based on a shapefile/polygon. I was able to do that easily in QGIS 2 but the package seems to be missing in QGIS 3. Any hints?


texts: PostGIS function to move a polygon to centre over new point coordinates

texts:
Background:

(Note: I am trying to do as much of this as possible in PostGIS/PostgreSQL)

I have a PostgreSQL table containing location data of an aircraft in orbit around a city. Every second a new row is added with new position data.

From this table I have created a view which shows a constantly updating point of only the current position. This view therefore contains only one row.

The Task:

I have created a custom polygon in a different table, which acts as a buffer around the current position point. I want to reposition/move this polygon to always be centred over the constantly updating current position coordinates. Does anyone know how to do this? I have already experimented with rotate and scale. I am happy with them. It is moving the polygon that is confusing me.

I have tried creating a view that uses ST_Translate for the polygon, however from my understanding, this function only moves coordinates by a fixed amount from the origin. It doesn’t allow me to specify exactly what lon lat coordinates I want it to move to. Does anyone know a function/how to move a polygon to new coordinates, centred on a point? Again, I am trying to achieve this from PostGIS.


texts: AttributeError: 'NoneType' Connecting GeoRaster

texts:
When I open the option to connect georaster and go again to create the connection, the error below appears. Do you know what it can be?

Show the mesage:

AttributeError: 'NoneType' object has no attribute 'findChild'
Traceback (most recent call last):
File "D:/Profiles/usr/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\oscars-georaster-qgis-plugin\oscars_georaster_dialog.py", line 112, in btnNew_clicked
self.dlgDataSourceManager.findChild(QtWidgets.QPushButton,'btnNew').click()
AttributeError: 'NoneType' object has no attribute 'findChild'

Version Python: 3.7.0 (v3.7.0:1bf9cc5093, Jun 27 2018, 04:59:51) [MSC v.1914 64 bit (AMD64)]
Version QGIS: 3.18.2-Zürich Zürich, 515138c171


texts:
I'm scripting OGR functionality in order to filter some mapinfo tables and convert them into sqlite files.

ogr2ogr -f "SQLITE" newtest.sqlite test.tab -sql "SELECT cols FROM test WHERE .."

where TAB file is about 2GB size (including ID,MAP,DAT & TAB extensions within)

I receive error like this :

ERROR 1: CPLRealloc(-1532233231): Silly size requested.

What does it mean exactly,how to resolve this ?


texts:
This is a follow up to Where to obtain free 15m ASTER image for New Mexico?.

If I want to analyze 10 years or so, do I need to download each ASTER file for each year (e.g. 2003 ASTER, 2004 ASTER,...) and then merge all of them into one image?

An ASTER revisits a specific area twice a month, so do I need to download two ASTER each month for 2003, 2004, 2005, etc. until 2013 and then run ndvi for each ASTER per month?

I do need to perform them crop analysis/ or soil crop health for this from 2003 to 2013.


texts: Visual Basic Error while Creating a new UIButtonControl

texts: I am trying to add a new form that a programmer sent me for a specific UIButtonControl. When I import the form into Visual Basic following the steps they sent me, it imports the form as a module instead of a new form. when I try to use the new UIButtonControl button that I created with the module, it gives me this error: Compile Error. Invalid outside procedure. Both the programmer and myself feel it is because visual basic is bringing the form in as a module instead of a form, can anyone help me fix this please?


texts:
I am investigating using some C/C++ libraries for point-in-polygon lookups and libspatialindex seems popular. I've written some code that seems to do what I want but I'm not entirely sure if what I doing is correct. I have a requirement that the polygons have an arbitrary number of dimensions, they generally range from 5 to sometimes 16.

My problem is that when creating the RTree I have to specify the number of dimensions and so if I try to insert a polygon that has a higher number the insertion fails. Also I'm not quite sure if SpatialIndex::Region is the correct object to use to represent polygons, it seems to actually be a rectangle rather that what I'd call a Region as it only allows 2 dimensions.

A quick fix would be to create individual Regions for each pair of points in the polygon but that is quite inaccurate.

I've scoured the internet for sample code but I can't find anything worthwhile.

Am I doing something wrong?


texts: How to choose proper edit session in arcpy?

texts:
There are 2 databases A and B; Each one has two tables: Table1 and Table2, both were created with the same script, defining for both: fields, versioning and archiving.
It is required to delete some rows in both tables and, of course in both databases.
A script was prepared:
def remove_items(from_table, ids):
    log.add("Start to clean orphan objects from: {0}".format(from_table))
    total_to_delete = 0
    fields = arcpy.ListFields("{0}.{1}".format(user_do, from_table))
    if "INSPECTION_ID" in [f.baseName for f in fields]:
        with arcpy.da.SearchCursor("{0}.{1}".format(user_do, from_table), ["INSPECTION_ID"]) as count_cursor:
            for row in count_cursor:
                if row[0] not in ids:
                    total_to_delete += 1

        row_counter = 0
        with arcpy.da.UpdateCursor("{0}.{1}".format(user_do, from_table), ["INSPECTION_ID"]) as table_cursor:
            for row in table_cursor:
                if row[0] not in ids:
                    table_cursor.deleteRow()
                    row_counter += 1
                    log.add("In table {0} deleted rows {1} / {2}".format(from_table, row_counter, total_to_delete))

            log.add("Total rows deleted from {0} : {1}".format(from_table, row_counter))

with arcpy.da.Editor(arcpy.env.workspace) as edit:
    log.add("This script may take a while to finish.")
    remove_items("Table1", orphan_ids)
    remove_items("Table2", orphan_ids)

The script is run against Database A, and it runs with out any issue.
The script is run against Database B, and it generates the following error:
['RuntimeError: Insufficient permissions [Database_B.Table_2][STATE_ID = 2880]\n', '\nThe above exception was the direct cause of the following exception:\n\n', 'Traceback (most recent call last):\n', '  File "C:/T/Repos/RemoveOrphanInspections.py", line 136, in <module>\n    main(argv[1:])\n', '  File "C:/T/Repos/RemoveOrphanInspections.py",'SystemError: <built-in method __exit__ of Workspace Editor object at 0x0000019EDF973180> returned a result with an error set\n']

Permissions were checked, and the data owner user (the one used to run the script), has the permissions of Select, Insert, Update and Delete.
I make a change to the script:
with arcpy.da.Editor(arcpy.env.workspace) as edit:
    log.add("This script may take a while to finish.")
    remove_items("Table1", orphan_ids)

edit = arcpy.da.Editor(arcpy.env.workspace)
edit.startEditing(False, True)
edit.startOperation()
remove_items("Table2", orphan_ids)
edit.stopOperation()
edit.stopEditing(True)

Now the script runs with out issue.
With this modification I am using two edit sessions, therefore there is two transaction scopes (That was not the original plan).
Versioning, archiving, permissions and fields are the same in both Databases.
How could check the differences between the Database_A.Table2 and Database_B.Table2?
Which property of the table should be checked to choose a proper Edit session management?


texts: Identifying number of facilities within service area in ArcGIS Desktop?

texts:
I am trying to identify the facilities within a Service Area. I have 100 points (facilities) where I created Service Area polygons of 5 km each. Now I need to know how many of these 100 facilities are within these Service Areas. In short, for service area of facility 1, what other facilities are within this area etc? The result needs to be a image of a service area with a table of all other facilities within this area.

I now do an intersect between the facilities and service areas which results in a point dataset of facilities with a column in which service area they fall in, but unfortunatly this also gives back all facilities that fall within their own service area. I have to delete the intersect of facility 1 with its own service area. I have a feeling that there might be an easier way to do this.


texts: Applying calc() on a rasterstack with certain pixels with all NA values

texts:
I have a raststack of 768 NDVI files with 4.1 million pixels each. Now I am performing a pixel based time series method (BFAST) using the following functions:

Function to convert daily ts into monlthy ts

aggregate.daily.to.monthly <- function(daily.ts) {
  s.month <- round(aggregate(as.zoo(daily.ts), as.yearmon, median), 0)
  s.month <- as.ts(s.month)
  return(s.month)
}

xbfast function to convert rasterstack into a ts object and apply BFAST

xbfast <- function(data) {
  s <- bfastts(as.vector(data), dates, type = c("irregular"))
  s.d.linear <- round(na.omit(s, method = "ie", interp = "linear"), 0)
  s.m.linear <- aggregate.daily.to.monthly(s.d.linear)
  result <- bfast(s.m.linear, h = 0.13, season="dummy", max.iter=2)
  return(cbind(result$Magnitude, result$Time))
 }

Performing calculation on rasterstack

Now I use calc() from raster package to apply xbfast on rasterstack ndvi which as dates function attached to it:

 fbfast <- calc(ndvi, fun = function(x){
      t(apply(x, 1, xbfast))
 })

  Error in .calcTest(x[1:5], fun, na.rm, forcefun, forceapply) :
  cannot use this function

The above error appears because certain pixels have all NA values and xbfast function does not work on such pixels. The above code is tested on few other areas and it produces results. Apparently, the error is given for those areas which have certain pixels with all NA values.

In order to skip calculations on those pixels which have all NAs I used na.rm argument but it doesn't run.

fbfast <- calc(ndvi, fun = function(x){
  t(apply(x, 1, xbfast)) }, na.rm = TRUE)

  Error in .calcTest(x[1:5], fun, na.rm, forcefun, forceapply) : cannot
  use this function. Perhaps add '...' or 'na.rm' to the function
  arguments?

Modified code using is.na

fbfast <- calc(ndvi, fun = function(x){
  if (all(is.na(x))) {
    NA
  } else {
    t(apply(x, 1, xbfast))
  }
})

  Error in apply(x, 1, xbfast) : dim(X) must have a positive length

I modified the code:

fbfast <- calc(ndvi, fun = function(x) {
  if (all(!is.na(x))){
    t(apply(x, 1, xbfast))
    } else {
      rep(NA, 2)
  }})

It runs and produces a rasterbrick of 2 layers but with only NAs! What could be wrong here?

I intend to produce three rasters (magnitude of breakpoints, timing of breakpoints and number of breakpoints). First 2 are part of xbfast but I couldn't find a way to produce number of breakpoints.


texts: Getting progressive distance polyline and points in QGIS

texts:
In QGIS, how do I get the distance between a layer of points and a polyline layer? I have tried with a distance matrix but it does not seem to be working,
I have a layer of points and a polyline; the perpendicular to the polyline must be found and then from the point of intersection the distance to the previous vertex must be calculated not the next. Is possible to have distance automatically?
I have tried with a v.net.connect but not working as it serves me. any ideas? I post this photo. I have to find d1, d2 and d3, the distance of the vertex.

This is the example I have export the distance of all points from the left vertex. I have tried with v.net.connect but the distance are from point to point that intersect. Any ideas?


texts: Split a single raster into multiple rasters based on the attributes in R

texts:
I have a land cover raster which has 15 land cover classes (which previously had 70 land cover classes and I aggregated the same). However, I realized that I need 15 different rasters based on each of the land cover types. Is there any way I can achieve this in R?

lulc_2015 <- raster("E:\\LandCover_Data\\2015_lulc")

#Load matrix to be used for reclassification
rec <- read.csv("E:\\LandCover_ReclassifyMatrix.csv")

#Reclassifying the Raster to convert 70 LC to 15 LC
rc <- as.matrix(data.frame(from=rec$V2, to=rec$To))

New_LULC_2015 <- reclassify(lulc_2015,rc)

New_LULC_2015 <- ratify(New_LULC_2015)
rat <- data.frame(
  ID= 1:15,
  LandCover = c("Evergreen","Deciduous","Mixed Forest",
            "grass","Degraded Forest","Scrubland",
            "Dry Grassland","Wet Grassland","Plantation",
            "Settlement","High Elevation Plantation","Coastal",
            "Barren Land","Cropland","Wetlands")
  )

  levels(New_LULC_2015)<- rat

 > New_LULC_2015
 class       : RasterLayer
 dimensions  : 642, 382, 245244  (nrow, ncol, ncell)
 resolution  : 1000, 1000  (x, y)
 extent      : 461951, 843951, 892583.3, 1534583  (xmin, xmax, ymin, ymax)
 coord. ref. : +proj=utm +zone=43 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
 data source : in memory
 names       : layer
 values      : 1, 15  (min, max)
 attributes  :
       ID LandCover
 from:  1 Evergreen
 to  : 15  Wetlands

If I try creating a new raster, I get this error.

> raster(New_LULC_2015@data@attributes[[1]]$ID==1)
 Error in (function (classes, fdef, mtable)  :
  unable to find an inherited method for function ‘raster’ for signature
 ‘"logical"’

I simply need 15 different rasters, for Evergreen, Deciduous etc. I need these as I realized I need to use the predict function from the raster package for a random forest model (Unless I don't really need to create separate rasters).


texts: Is there a way to share a QGIS plugin for some users, not every QGIS users?

texts:
I need to make a PyQGIS plugin for a client, but we have some specific connection to a database. I did a plugin before only for me so I used QGIS plugin reloader. Now I need it to be on other desktop.
My questions are:

If I put my plugin on the qgis repo, do I have any way to make it only usable for the client
What are the other options I have if not?


texts: Syntax for evaluating strings with wildcards

texts:
I'm trying to do a Field Calculation in ArcGIS 10.3 (Python 2.7) that uses the Codeblock (Pre-Logic Script Code) to evaluate the attributes in a different field to inform the calculation).

Here's an example of what I'm trying to do (but getting syntax errors):

def evalStreetType(TYPEFIELD):
    if (TYPEFIELD LIKE "St*"):
        return "Street"
    elif (TYPFIELD LIKE "Av*"):
        return "Avenue"
    elif (TYPEFIELD LIKE "P*k*wy*):
        return "Parkway"
    else:
        return TYPEFIELD

And then the field calculator will say [StreetType] =

evalStreetType(!AbrvStTyp!)

Hopefully I didn't introduce new syntax errors as I simplify things to post in this question. I believe the crux of my problem is that I'm phrasing my "LIKE" statement and using wildards in a SQL sort of way, and it's not how Python does it. But I'm not sure now to evaluate strings with wildcards in python (....obviously).

Can someone help me out!??


texts: I'm using QGIS to do georeference scan of the map. I have to do it in UTM zone 36N. The problem is I want to do it by creating control points and use  background maps available in Open Map Layers (Bings Aerial, Google Map etc.). They are in PseudoMercator- not possible to change coordinates. What is the best way to resolve this problem? I've already received adjacent area shp and georeferenced TIFFs so it's just continuation of the work. There is no exceptions, it has to be in 36 N.


texts: Projecting NLCD data as WGS84 using QGIS?

texts:
I am using QGIS 3.03 Girona and trying to process the CONUS 2016 NLCD (Continental US 2016 National Land Cover Database) 30x30 raster data provided by the MRLC viewer.
I downloaded an area with the boundaries W078 N40 E076 S38 and would like to use this landcover data with Open Street Map vector data and SRTM-1 elevation height maps. That data comes in the WGS84 projection while the NLCD data comes in its own projection:

The CRS given by QGIS is:

Generated CRS (+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs)

However, I have been having trouble reprojecting the data.
If I set WGS84 directly by Set Layer CRS, the raster layer disappears entirely. Zoom to layer shows its now projected fine but massively zoomed out: notice the coordinates at the bottom - thats in degrees!

If I save as with WGS84 set, I get the same (except in black and white).

If I use the warp tool to reproject as WGS84, I get the very same that I started with (the image on top of the OSM data correctly scaled, but slanted over).

Does anyone know how to project it properly?

It should cover the exact same area as my four SRTM-1 elevation tiles.


texts: Display local language for OSM data in TileMill

texts:
I'm following the instructions here to create a minimalist style map for Japan:

https://github.com/aaronlidman/Toner-for-Tilemill

Now, I've downloaded the OSM data, used osm2pgsql, created the tilemill project and it now properly displays in TileMill.

However, all the displayed city name labels are in the Romanized (english) characters.  Is there a way that I can get the Japanese (Kanji, hiragana, etc) displayed in place of thier romainzed dipictions? (Street names do appear in the original Japanese, sometimes with Romanized letters in parens)

I'm assuming that the original Japanese names are contained in the osm.pbf, but somewhere along the pipeline they are being ignored.

Any ideas?


texts: Automated geodatabase replication with offline editors

texts:
I would like to design an application that would use an SDE database with simple features and a geometric network. That geodatabase would then be replicated out to any number of personal geodatabases that would be edited using an ArcEngine application. These engine users could be offline at any point and ideally the changes that were made would be sync'ed/merged in both directions upon reconnection to the network without any versioning or checking in/out data.

Can SDE replication be configured to accomplish this scenario? What other alternative is there to make this happen?


texts: Debugging error 000278 from Append and Merge tools in ArcGIS Desktop?

texts:
Append and Merge tools issues in Arc 10. keep getting error 000278 - failed on input OID ...... (could not write value to one field or the other).

Does anyone know how to fix this?

Never had this problem in Arc 9.3, where if the value did not fit, it would get chopped off.


texts: Extending size (area) of raster using ArcGIS Desktop

texts:
My problem is that I would like to multiply two rasters using raster calculation tool. Both rasters have values 0 and 1, but the covered area of first raster is bigger (100 x 100 miles) then the area of second (80 x 80 miles).

How can I extend the second raster to the same size?

I would just like to add columns and rows with 0. At first I thought that this would be simple, but after searching on internet I didn't found any solution.


texts: How to change delimiter in attribute table of QGIS?

texts:
I have imported a .kml file into QGIS, and the data is all shoved together into a single field like this:

Fine, it's pipe delimited - easy to change the delimiter like in LibreOffice calc imports I thought. But could not find a solution within QGIS.

Next attempt: change the delimiter in the kml file to the common default comma delimited using Sed: sed -e 's/|/,/g' in.kml > out.kml

Still no joy - it looks exactly the same as above, only with commas instead of pipes. So someone please help me, and I'm sure it will be of use to others.

(P.s. a colleague and I have tried the same thing in ArcGIS and the same problem is encountered).


texts: GeoJSON vector layer labels are not fixed in respect to WMS layers

texts:
I have three layers on my map. Two are from a WMS service (orthophoto and green linestrings) and one is vector layer (points - yellow labels) with GeoJSON source:

OL3Source = new ol.source.GeoJSON({
  //SRID: 102067
  'url': mapp['maprouter']+'?act=get_layer_data&layer='+layer_id+url_params
});

The GeoJSON looks like this (there are 3 times 2 same points - it's an issue in our database):

{"type":"FeatureCollection","features":[{"type":"Feature","geometry":{"type":"Point","coordinates":[-642906.776502149,-1040547.16261086]},"properties":{"feature_data_getter_route":"get_cf_data_overview","feature_data_container":"popup","layer_id":"celni_fotky.overview","id":2}},{"type":"Feature","geometry":{"type":"Point","coordinates":[-647061.5440096,-1036457.97874901]},"properties":{"feature_data_getter_route":"get_cf_data_overview","feature_data_container":"popup","layer_id":"celni_fotky.overview","id":1}},{"type":"Feature","geometry":{"type":"Point","coordinates":[-642906.776502149,-1040547.16261086]},"properties":{"feature_data_getter_route":"get_cf_data_overview","feature_data_container":"popup","layer_id":"celni_fotky.overview","id":7}},{"type":"Feature","geometry":{"type":"Point","coordinates":[-647061.5440096,-1036457.97874901]},"properties":{"feature_data_getter_route":"get_cf_data_overview","feature_data_container":"popup","layer_id":"celni_fotky.overview","id":6}},{"type":"Feature","geometry":{"type":"Point","coordinates":[-642906.776502149,-1040547.16261086]},"properties":{"feature_data_getter_route":"get_cf_data_overview","feature_data_container":"popup","layer_id":"celni_fotky.overview","id":12}},{"type":"Feature","geometry":{"type":"Point","coordinates":[-647061.5440096,-1036457.97874901]},"properties":{"feature_data_getter_route":"get_cf_data_overview","feature_data_container":"popup","layer_id":"celni_fotky.overview","id":11}}]}

The label is set as follows:

OL3Layer = new ol.layer.Vector({
            source: OL3Source,
            visible: true,
            style: styles.get(layer_id)
        });

where

new ol.style.Style({
 image: new ol.style.Icon(({
   anchor: [0.5, 1.5],
   anchorXUnits: 'fraction',
   anchorYUnits: 'fraction',
   src: 'images/map/photo-marker.png'
 }))
});

and the view is set like this:

view: new ol.View({
 center: ol.extent.getCenter(mapp['CR_feataure'].getGeometry().getExtent()),
 zoom: 8,
 projection: ol.proj.get('EPSG:102067')
})

The projection is enabled like this:

proj4.defs("EPSG:102067", "+proj=krovak +lat_0=49.5 +lon_0=24.83333333333333 +alpha=0 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel +towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56 +units=m +no_defs");

The projections are all EPSG::102067. WMS is provided by MapServer. Coordinates for the vector layer are converted on geometry from lat/lon using the following query:

UPDATE cf_metadata SET geom = ST_Transform(ST_SetSRID(ST_MakePoint(degrees(lon),degrees(lat)), 4326),102067);

The problem is that when I zoom in or out of the map. The WMS layers stays fixed in relation to each other but the vector labels are not fixed at the point and are gently sliding away from their position.

Zoomed out:

Zoomed in:

Do you have any idea how to fix this?


texts: Dynamically viewing Python script results in ArcMap?

texts:
I want to view the output of ArcPy scripts (run in Python IDLE) in ArcMap.

I tried refreshing the ArcMap data frame after I run the script.

This does not work. Please advise me on this?


texts: st_read option does not work if parameter is separated from equal (=) sign

texts:
This is a small pet peeve. I'm reading a latin1 encoded shapefile in R using st_read and the encoding option.

I just found out that if I write encoding=latin1 everything works as expected but, if I write encoding = latin1 I get a warning message and the encoding is not set correctly.

Is this the expected behavior? Has anyone else noticed this?


texts: How do you symbolize "Unique Values, Many Fields" with more than 3 fields?

texts: I have data that has 5 fields with unique values that I need to symbolize at one time. Is there a work around to allow me to do this or a python code I can run to use more than three fields?


texts: Using shapefiles as a reference for accuracy assessment in ecognition

texts:
I classified my project and got an ESRI shapefile with points in it. Those points are hand made samples with the right class assigned to them. Now I need to use those samples to check my own classification via AA - but how do I use this shapefile?

Loaded it as a thematic layer but I can not interact with it in the AA section of eCognition.


texts: How to find out what projection my data is using?

texts:
I have a set of data that I created from gps locations of images. The locations are in decimal. I don't know what projection was used when the data was captured.
I used QGIS to display the data with a forced projection (EPSG:28355 - GDA94 / MGA zone 55) but the data looks wrong.
The resulting map is off by nearly 90 degrees and the circular streets are skewed. So I am assuming that the projection is wrong and I can't work out how to change it.
Can someone help?


texts:
Recently I had ArcGIS Pro crashes because of an add in I installed. I followed the instructions for a clean uninstall of Pro and reinstalled it. But after re-install the add in still caused pro to crash.
Where does Pro keep the config settings for addins?
I was attempting to tweak some settings and now every time the add in is loaded pro will crash. I've downloaded a new add in file, I've removed and reinstalled the add in, and I've done a clean uninstall of pro but nothing works.
The add in:
https://github.com/USEPA/EPA-Metadata-Editor-Pro/wiki


texts: Retrieving Elevation from Color in Raster DEM?

texts:
I'm fairly new to GIS work.

I have a DEM in raster format, that is unfortunately not properly meta-referenced. In other words it has a legend where each color is assigned with the respective elevation, but no program I have could work with it.
So what I think I need to do is, assign each color the respective elevation value and then turn it into a proper DEM to work with (the terminal result should be a *.XYZ file).
But I don't really know how...

Programs I use(d):

QGIS
Global Mapper
JOSM
L3DT

The DEM file with and without legend:

With Legend: Download (2 MB)
Without Legend: Download (7 MB)


texts: Calculate intersection area in square meters of two polygons that use EPSG:4326 in Python


texts: Which version of Spatialite does QGIS support?

texts:
At work, we utilize QGIS 2.8.3. I am interested in some features of SpatiaLite, as styling support (since version 4.3.0) and topology support (since 4.4.0, a Release Candidate version).

Where can I found which version of SpatiaLite a specific version of QGIS supports?


texts: How to control what kind of tags / flags are imported with osm2po?

texts:
With the help of osm2po, I have successfully inserted the contents of the italian files (italy.osm.pbf) in a pgrouting-enabled database (postgresql 8.4)

It displays fine in Qgis, using rule based labeling (thanks to underdark for the .qml file).

Now I wonder: how can I find out what attributes are included in the .pbf and how to control what tags / flags are imported in the database? Any help very much appreciated.


texts: PyQGIS Standalone QVectorLayer works and QRasterLayer not work without QgsApplication / initQgis()

texts:
I'm trying to make a standalone python app and I want to use in some parts of app functions of PyQGIS.

I'm testing things to understand how it works.

Doing tests I have found a situation that I do not understand.

I can load a QgsVectorLayer without QgsApplication.initQgis(), but with QgsRasterLayer this not works.

Can someone explain why this happends? (I know it is recommended to start always with QgsApplication.initQgis().

Here is my test code.

Works QgsVectorLayer but not QgsRasterLayer.

from qgis.core import QgsVectorLayer,QgsRasterLayer

file='MyTrack.gpx'+"|layername=track_points"
vlayer = QgsVectorLayer(file, "Termes", "ogr")
if vlayer.isValid():
    print('Load Correct QgsVectorLayer')
else:
    print('Not Load QgsVectorLayer')

for feature in vlayer.getFeatures():
    print(feature.geometry().asWkt())
    break

uri="dpiMode=7&identifier=mdt:Elevacion4258_25&url=http://www.ign.es/wcs/mdt"

DEM  = QgsRasterLayer(uri, 'my wcs layer', 'wcs')
if not DEM.isValid():
    print("QgsRasterLayer failed to load! ")
else:
    print("QgsRasterLayer Load Correct!")

Works QgsVectorLayer and QgsRasterLayer (take about 1 minute to load)

from qgis.core import QgsApplication
qgs = QgsApplication([], False)
qgs.initQgis()

from qgis.core import QgsVectorLayer,QgsRasterLayer

file='MyTrack.gpx'+"|layername=track_points"
vlayer = QgsVectorLayer(file, "Termes", "ogr")
if vlayer.isValid():
    print('Load Correct QgsVectorLayer')
else:
    print('Not Load QgsVectorLayer')

for feature in vlayer.getFeatures():
    print(feature.geometry().asWkt())
    break

uri="dpiMode=7&identifier=mdt:Elevacion4258_25&url=http://www.ign.es/wcs/mdt"

DEM  = QgsRasterLayer(uri, 'my wcs layer', 'wcs')
if not DEM.isValid():
    print("QgsRasterLayer failed to load! ")
else:
    print("QgsRasterLayer Load Correct!")


texts:
I  would like to determine the surface ﬂow network, this based on the spatial analysis of DEM data and the creation of separate layers for surface ponds and preferential surface ﬂow pathways.
More accurate I would like delineation surface pond and surface flow pathway over urban areas based on DEM.
I want to know what software I can automatically calculate?


texts: How to use SEXTANTE in a script outside QGIS?

texts:
I´m trying to do some geoprocessing using the QGIS Sextante plugin outside of the QGIS application (and integrate it with different other libraries like arcpy, etc.).

import sys, os
from PyQt4.QtGui import QApplication
QApplication(sys.argv)
from qgis.core import *
import qgis.utils

sys.path.append("C:/Users/.../.qgis/python/plugins")

import sextante
sextante.core.Sextante.Sextante.initialize()
input = "D:/Python_Test/a.shp"
output = "D:/Python_Test/b.shp"

sextante.runalg("qgis:convexhull",input,None,None,output)
print "Finished!"

Unfortunately the process fails with the following error:

Traceback (most recent call last):
  File "<interactive input>", line 1, in <module>
  File "C:/Users/.../.qgis/python/plugins\sextante\core\Sextante.py", line 376, in runalg
    alg = Sextante.runAlgorithm(algOrName, None, *args)
  File "C:/Users/.../.qgis/python/plugins\sextante\core\Sextante.py", line 276, in runAlgorithm
    if not param.setValue(args[i]):
  File "C:/Users/.../.qgis/python/plugins\sextante\parameters\ParameterVector.py", line 59, in setValue
    layers = QGisLayers.getVectorLayers(self.shapetype)
  File "C:/Users/.../.qgis/python/plugins\sextante\core\QGisLayers.py", line 79, in getVectorLayers
    layers = QGisLayers.iface.legendInterface().layers()
AttributeError: 'NoneType' object has no attribute 'legendInterface'

What am I doing wrong? Sextante seems to be initialized correctly because

sextante.alghelp("qgis:convexhull")

gives me

ALGORITHM: Convex hull
    INPUT <ParameterVector>
    FIELD <ParameterTableField from INPUT>
    METHOD <ParameterSelection>
    OUTPUT <OutputVector>

I´m using Windows 7 Professional 64Bit and the OsGeo4W installation with QGIS 1.8.0 and python 2.7.4.


texts: ArcMap Zonal Statistics as Table 2 output does not contain all zones

texts:
I'm encountering an issue with zonal statistics in which my output table contains less zones than my input (the output contains roughly 1/3 of the input zones). It is a large data set with about 65,000 zones (vector) and a raster grid of population density with a cell size of 5 meters. There is a lot of overlapping among the zones, but I am using Zonal Statistics as Table 2, which should be able to handle that. From other users' accounts of this issue, I've come to understand that zones are excluded in the output table if they do not overlap any of the cell centers in the raster grid. However, in my data set, all zones are >800 meters in diameter and the cell size of the raster grid is 5 meters, so I don't think cell size can be the root of the issue here. I also have made sure in the environment settings to set the cell size to 5 meters, set the snap raster as input raster grid, and to set an appropriate boundary and mask. I'm using ArcMap Desktop Version 10.5.

Does anyone know what could be causing this issue?


texts: Qgis problem clipping a raster .asc file, the output is unreadable

texts:
I am currently using QGIS (1.8.0) to prepare files for use in MaxEnt to predict species distribution in a new territory in Argentina. I need to

a) cut the raster files (bioclim, and others) to a single province, and

b) create a .asc file from a .shp file (not sure if this is possible in QGIS)

My problem is with point 'a'. I have read the various questions on this site regarding raster clipping which were very useful. I have set the CRS for each file (WGS84 EPSG4326), and the Raster/Extractions/Clipper tool using a mask does create a file.

However, when I open this file using TextPad (with which I can view the original larger file) the text is in an unknown language/script and does not open in MaxEnt.

I was worried the problem my lay in the fact that my .shp mask did not fully fall within my original raster, so I tried clipping by extent, so I clipped only within the extent of the raster, but got the same results.
Has anyone any suggestions as to what might be the problem?


texts: How can I convert from EUREF-FIN coords to Lat/Long?

texts: See question title, how can I translate EUREF-FIN coords to Lat Long?


texts: Multiple users overwrite each other's changes in One Drive shared file geodatabase

texts:
We are using file geodatabases along with Micorsoft teams one drive.  The beauty of this is that users all have the need files on their C drives and it is super fast.  When a user edits a file it automatically gets updated on everyone elses PC.  When we try going off the network it is painfully slow.  So hence why we are not using Enterprise GDB along with the additional cost.

Where I am running into an issue is that if two users are trying to create points at the sametime in a feature class with the file geodatabase and then save the point it saves on their version on the C drive.  However, because they both are using the same object ID which everyone makes it to the server gets posted.  So if both users create a point because objectID goes in incremental order it would be objectID #7 for example.   Hence the overwriting.

As I understand it you cannot have two users writing to the same feature class at the same time even if they both were using the same geodatabase.

So my question:
Is it possible to change the objectID from being incremental to something else like maybe a globalID or allow us to programmatically assign the ObjectID. Can I make it so the ObjectID is not the unique identifier but something else would be.  How can I work around this issue?


texts: Automatic sheet size in text label in QGIS

texts:
I'm in Print Composer and want to develop a template. In that template, I would like to add a text Label which automatically represents the size of the sheet in the format "Sheet size: 841x1189".

There is got to be an expression for that but I can't find it. Can someone recap it for me?


texts: Help with Workflow for an Enterprise GeoDatabase Management

texts:
[This is a cross post from ESRI's forum--not even getting many views there. Apologies for the duplicate.]

Please see this for reference: http://img42.com/YA1Qo

I have been tasked with creating an ESRI Javascript API application; but I could also have access to ASP.NET and/or php environments if needed. The attached file show a data model I have been working on. It shows only a few fields. It is for a university. Briefly, here is what the model describes: The 'Buildings' table will have all the buildings on campus. The 'RoomData' table will have all Rooms from all Buildings. The 'FloorPlans' table will have all Floor Plans for every Floor of every Building. Let's forget the Insurance table for now.

The database will be SQL Server 2012. The Server OS will be Windows 2012. ESRI's ArcGIS for Server 10.3 will be the GIS Server.

The challenge before me is to enable non-technical staff to Add a building to the enterprise geodatabase and then add individual Rooms for that Building, and Floor Plans for that Buildings for each Floor of that Building.

Solution 1: I foresee a web-based solution as better: A 'wizard' which would guide the users through adding data through online forms. Not a problem. But the problem is Spatial data. Specifically, how can each Building's polygon be made into the Building's table in a web interface? The users could have an ArcMap document where they could draw the Polygon but how could that be sent to the data's correct field? Maybe some python script?

Solution 2: ArcMap based solution: There can be two MXDs/Building. Each MXD can have all the Rooms (or Floor Plans) data for a specific building. So these two MXDs would cover all the required data for a Building. There can also be an additional field called 'Building_ID' with hard-coded value for a Building (which would exist in another table in the geodatabase). This is the current solution in a 'Pilot' project and kind of works: Each of the two MXDs get Published to the GIS Server--through an interface with SQL Server (ArcSDE). However, this is not a true relational database.

Let me add that I may be completely off or maybe over-complicating matters. I thought a wizard based approach would be better but am more than willing to listen to other ideas.


texts: Add column to show intersection with another layer

texts:
I have a table containing a list of properties.
I have a successful query to return only those that intersect another layer.
However i want to return all the properties now and instead return a new column that shows InBuffer = True or False.

This query works to get only the intersections but dont know how to change it to what i need.  Any ideas?

SELECT p.[FEATURE_ID] AS [FEATURE_ID],
  p.[Transactionuniqueidentifier] AS [Transactionuniqueidentifier],
  p.[geogtab] AS [geogtab],
  p.[Price] AS [Price],
  r.[linenum] AS [linenum]
FROM [dbo].[GB_PP_HS2_1617] p
inner join [dbo].[GB_HS2_ROUTE_C] r
on r.railroute.STBuffer(500).STIntersects(p.geogtab) = 1

In this example geogtab is a geography column.

If this is easier to perform in POSTGIS im happy to use that instead.

I have tried the following without success.... query runs forever without finishing

    SELECT p.[FEATURE_ID] AS [FEATURE_ID],
 p.[Transactionuniqueidentifier] AS [Transactionuniqueidentifier],
 p.[geogtab] AS [geogtab],
 r.geog.STBuffer(500).STIntersects(p.geogtab) AS [InBuffer],
 r.[linenum] AS [linenum]
FROM [dbo].[GB_PP_HS2_1617] p
left join [dbo].[GB_HS2_ROUTE_C] r
on r.geog.STBuffer(500).STIntersects(p.geogtab) = 1


texts: How to create raster from data (x,y,z in Cartesian coordinate System)?

texts:
I have simulated many trajectories of an moving object. From simulations I have obtained data files where the trajectory is saved in terms of x,y,z coordinates recorded at equal time intervals. Each recorded position also indicated the kinetic energy of my object. Sample data for first trajectory is below. All my data has been recorded in Cartesian coordinate System

      t        x      y        z     Etot
1  0.00 696621.4 167730 1680.960 1192.526
2  0.01 696621.4 167730 1680.959 1192.526
3  0.02 696621.4 167730 1680.958 1192.526
4  0.04 696621.4 167730 1680.952 1192.526
5  0.06 696621.4 167730 1680.942 1192.526
6  0.08 696621.4 167730 1680.929 1192.526
7  0.11 696621.4 167730 1680.901 1192.526
8  0.13 696621.4 167730 1680.877 1192.526
9  0.15 696621.4 167730 1680.850 1192.526
10 0.16 696621.4 167730 1680.834 1192.526

Now, I would like to create a raster file from the simulated data (2m DEM). The desired output is the object spatial frequency, aka assign a probability to each pixel of the object passing though it?

I have been trying to work by way through the Applied Spatial Data Analysis in R book, but i get stuck, i have never worked with spatial data before.

QUESTIONS:

whats the easiest way to achieve my objective?
which R functions/packages should I use?
do I need to assign a CRS ?

EDIT (adding info) :
x,y,z ale coordinates from a 3-D cartesian coordinate system and describe the trajectory of a moving object. I have discarded a z dimension (height) and only plotted the positions of the rock within a 2-D system :

I have 90 more trajectories (within this area). I want to visualize the spatial deposition patter -> see which pixels are more likely to be hit by the object or not.

I was thinking of proceding:

Create an empty raster cell with 2m resolution
Overlay the 90 trajectories obtained from simulations (similar to Figure 1) and record how many times (and ideally with what energy) does a pixel get hit or trespassed.
Interpolate to empty cells
Resultant raster would express rockfall spatial frequency


texts:
I have this statement working in arcpy but I want this to implement in gdal, especially with gdal_calc, if possible.

Rastercalc = Con(((C-B) > 2) & (B > 100), A, Con(((A-D)>15) & ((A-E)>15),0,A))

I already tried writing this in gdal_calc as below but the output is not the same

gdal_calc.py --format GTiff --type Byte -A %s --A_band %s -B %s -C %s -D %s --D_band %s -E %s --E_band %s --outfile %s --calc="A*((C-B) > 2)+ (A*((A-E) < 30)* ((A-D) < 30))* ((C-B) <= 2)" --overwrite

I even tried this using numpy.where() with no luck

outcalc = numpy.where(numpy.logical_and((C - B) > 2, B > 100), A, numpy.where(numpy.logical_and((A - D) > 30, (A - E) > 30), 0, A))

Would somebody knows if this possible in gdal_calc (this is preferred way) but any other option in gdal would be great including numpy.where()?

I am pasting the part of my script and sampled results. The input data is actually a multiband (14 bands) dataset.

script:

--arcpy implementation----

            A = Raster(infile + os.sep + 'Layer_'+str(b))
            D = Raster(infile + os.sep + 'Layer_'+str(b-1))
            E = Raster(infile + os.sep + 'Layer_'+str(b+1))
            B = Raster(tempMed)
            C = Raster(tempMean)
            Rastercalc = Con(((C-B) > 2) & (B > 100), A, Con(((A-D)>15) & ((A-E)>15),0,A))
            Rastercalc.save(out_band)

            ##------ gdal_calc implementation------

            runCalc     = 'gdal_calc.py --format GTiff --type Byte -A %s --A_band %s -B %s -C %s -D %s --D_band %s \
                        -E %s --E_band %s --outfile %s --calc="(((C.astype(int16)-B > 2) * (B.astype(int16) > 100)) *\
                        ((A.astype(int16)-D <= 15) * (A.astype(int16)-E <= 15))) * A" --overwrite'\
                        % (infile, b, tempMed,tempMean,infile, b-1,infile, b+1, out_band)

spectral profiles of original layer, layer produced in arcpy, and gdal.

original

arcpy

gdal


texts: How to create a barrier feature class to use in a geometric network?

texts:
I am working with the NHD and need to isolate stream segments upstream of barriers.  I have too many barriers (thousands) to do this manually using the Geometric Network barrier tool on the toolbar. How can I create a geometric network feature class (or something) of these barriers?
Thanks!
katie


texts: Getting layer name of selected features with ArcPy?

texts:
I'm working on Python add-in that changes geometry of an object (polygon orthogonalization). The tool itself works properly when I give it an input line of layer name (or just drag layer from table of contents to script window).

Now I need to deal with single selected feature from unknown layer. For example I have 10 layers in the current mxd, edit session started for 5 of them and I need the one which feature is selected from.

How to get a layer name of current selected feature?

UPD. So what i have now. When tool works from toolbox it have an input parameter layer which is further converted to path:

layer = arcpy.GetParameterAsText(0)
dsc = arcpy.Describe(layer)
layerName = dsc.baseName
fd = dsc.path
layerPath = os.path.join(fd, layerName)

This works good when I execute this tool from toolbox for all the layer contents. This tool can also work with only selected features and I already tested it by replacing all the input with a strict line layer = "building" and adding this layer to TOC.

Now I'm trying to edit this tool to work as add-in with button.

What I need now is to know how to get a layer name when I simply select one polygone from some layer.


texts: Why are my circles plotted on maps, using buffer data from GeoPandas, not the same size as if I plot on Google Earth?

texts:
I am using GeoPandas and matplotlib to plot some circles of desired radius around a set of points.
I find that the circles are not the right size on the maps, but the equivalent sized circles drawn on Google Earth are correct.
Here is the general code sequencing.
points_wgs84 = gpd.GeoDataFrame(
    pdf_frame,
    geometry=gpd.points_from_xy(pdf_frame.longitude, pdf_frame.latitude),
    crs={'proj': 'latlong', 'ellps': 'WGS84', 'datum': 'WGS84', 'no_defs': True}
).rename(
    columns={'geometry': 'measurement_coordinates'}
).set_geometry('measurement_coordinates')

points_3857 = points_wgs84.to_crs(epsg=3857)

Here are the radius values for the circles to be plotted:
143     3120.545547
1430    3120.545547
2421    3120.545547
4789    3126.559006
2011    3126.559006
3592    3126.559006
59      3133.507389
4835    3133.507389
3549    3133.507389
1174    3133.507779
2225    3133.507779
2639    3133.507779
4481    3154.907523
2359    3154.907523
1901    3154.907523
2144    3176.645992
4073    3176.645992
5010    3176.645992

points_3857['disc'] = points_3857.buffer(points_3857['radius'].astype(float))
points_disc_3857 = points_3857.set_geometry('disc', drop=False, inplace=False, crs='epsg=3857')

points_disc_3857['circle'] = points_disc_3857.boundary
points_circles_3857 = points_disc_3857.set_geometry('circle', drop=False, inplace=False, crs='epsg=3857')

Then plotting the result with the code:
plt.rcParams["figure.figsize"] = (20, 20)
ax = points_circles_3857.plot(alpha=.3)
ctx.add_basemap(ax)
points_3857.plot(ax=ax, c='r',alpha=0.3)

I see a map that looks like this:

The circle are supposed to be in the range of 3120 - 3180 m in radius.
Now comparing, I find that the equivalent sized circle on Google Earth measure at 2600 m.
I used the intersection of the roads at about 170 degrees as the common reference to make the measure.
I believe that Google Earth uses true Mercator, and that 3857 is Web Mercator, but I cannot figure out how to get the measurements in GE and GeoPandas plotting to match.
I am looking to understand if I have the wrong coordinate transformations for the CRS, and how I can fix the relationship match.


texts: Creating report and inserting comma at end of it in ArcGIS Desktop?

texts:
When I run a create report and want to use the comma when you want to sum the acres ? Is there a way around to do that ?

Here is the attachment where I use the yellow highlight.

It is so much easier if you just add the comma !

I am using ArcGIS ArcMap 10.0 SP4.


texts: Filter data with PostgreSQL or geoserver?

texts:
I'm working on a PostgreSQL database with PostGIS. I need to filter an 8 million records table in order to publish a view through geoserver.

I know two ways to do that, one is that I create view in PostgreSQL to filter the data and publish on geoserver, the other one is let the geoserver making the filtering using sld.

My question is that which one is the better solution regarding to gain better performance.


texts: I have a text file with consecutive coordinates of gps data  for ongoing and incoming route(on Same road). i want to convert them in line format. how can i normalise those extra coordinate so that i get only single line for a single route.


texts: Incorporating several Python script tools into one tool and have user input determine section of script run?

texts:
I have a series of script tools created and work well in ArcMap 10.2.2 just fine.  These tools generate a series of buffer features from an input feature class.  I have these as separate tools each tool with a different set of buffer distances.

I want to be able to make one tool where the input feature class is chosen, then the set of buffer distances chosen out of a list from the script tool UI.

Example below of what is needed

import arcpy
arcpy.env.workspace = "C:\Buffers.gdb"
InPut = arcpy.GetParameterAsText(0)

#"select buffer situation B or A (not optional must be A or B)"

# Situation A Buffer 1
arcpy.Buffer_analysis(InPut,"first_Buff","16 METERS","FULL","ROUND","ALL")
# Situation A Buffer 2
arcpy.Buffer_analysis("first_Buff","second_Buff","5 METERS","FULL","ROUND","ALL")

# Situation B Buffer 1
arcpy.Buffer_analysis(InPut,"first_Buff","20 METERS","FULL","ROUND","ALL")
# Situation B Buffer 2
arcpy.Buffer_analysis("first_Buff","second_Buff","10 METERS","FULL","ROUND","ALL")


texts: Retaining Older Features with OpenLayers BBox Strategy

texts: To limit the number of features requested, I have used OpenLayers BBox Strategy with ratio of 0.5. The problem is, when I pan the map a new set of features are displayed and the older features are removed. Is there any way to keep those older features when new features are requested?


texts:
I'm using a tile service that has different imagery sources at different zoom levels e.g. when zoomed out clouds and older imagery are shown, but as you zoom in to higher resolution imagery the clouds are lost and imagery is more recent.
Ideally I'd like to provide a tile service where the imagery is consistent through all zoom levels. Are there any tools I could use to generate or proxy a tile collection based on specific zoom levels from the original tile source, interpolating the other zoom levels from the desired levels?
I've had a good look around MapProxy as a potential solution, but I couldn't see any useful options there for this usage case.
If I can't find anything I intend to code a solution up to generate the tiles, but I imagine this is a problem that has already been solved somewhere by one or more tools that I've not encountered yet.


texts: How to correct topology area boundary errors (slivers/gaps) on arcs in arcgis / arcsde

texts:
I use arcgis 10.1 (arcinfo licence) with arcsde to connect to a SQL Server database.  The initial dataset has around 100k polygons (used for land registration, having a resolution that is around 1cm and using an MTM projection).  Many geoprocessing operations were made on the dataset to build several related layers (copy / paste / split / merge / join,  etc.).  Note:  I never did any polygon simplification/generalization that would change the shape or amount of vertices in the polygon layers.

At the end of these operations I made a topology (1mm tolerance) between my initial dataset (rank #1) and the final dataset (rank #5) using the "Area Boundary must be covered by boundary of" rule to see if I landed back on my feet.  I expected errors only on the polygons that I had modified.  However, I got many (several thousand) errors on parts that I had not modified.  The symptoms are always the same:  they always affect polygon arcs.  The width of the sliver/gaps is always between 0,01mm and 10mm. Vertices between my layers are always identical (looking at 1:1 zoom with X,Y coordinates) but the path between vertices is changed for some unknown reason.  I have enclosed a few pictures to show the problem (the zoom pictures are at 1:1 zoom).  Any idea about what could have happened and how to correct this?  The main constraint is that the initial layer must remain unchanged (for legal reasons).

I have tried to change the tolerance of the topology to 1cm - it does solve 50% of my problem, but at the expense of the initial layer being slightly modified (I made an intersection between the original layer and my "initial" layer and they revealed many small changes).

I have also tried the automated "snap" tool using "Edges" and "Vertex" options, but almost nothing was changed (probably because the vertices are already in the right place.  Only the arcs between vertices have problems).


texts: QGIS 3 doesn't recognize Description from Value Map in attribute table

texts:
I have some .csv files that defines Value Maps. When I set them in the Attributes Form panel, it shows the right values and description names according to the .csv files.

But as can be seen on screenshot, in attribute table are not shown description values, but only raw values. I searched at almost every settings panels, but I can't figure it out.

In QGIS 2.x it shows description values by default. Does anybody know whats the point of this problem?


texts: How to add graphics to map centre in OpenLayers?

texts:
I am using the following code to add point graphics to the map center whenever there is an extent change, how ever the point remains at [0,0] and doesnt change with change in extent. How can I fix this?

map.on('moveend',function(evt){
    var map=evt.map;
      var extent1=map.getView().calculateExtent(map.getSize());
    var b= getCenter(extent1);
    var c=transform(b, 'EPSG:3857', 'EPSG:4326');
    var lon=c[0].toFixed(1);
    var lat=c[1].toFixed(1);
    var point_feature = new ol.Feature({ });
    var point_geom = new ol.geom.Point(
      [lon,lat]
    );
    point_feature.setGeometry(point_geom);
    var vector_layer = new ol.layer.Vector({
      source: new ol.source.Vector({
        features: [point_feature]
      })
    })
    var style = new ol.style.Style({
      image: new ol.style.Circle({
        fill: new Fill({
          color: '#3399CC'
        }),
        stroke: new Stroke({
          color: '#fff',
          width: 4
        }),
        radius: 8
      }),
    });
    map.addLayer(vector_layer);
    vector_layer.setStyle(style);
  })


texts: Preventing user from selecting different tool when custom tool is selected?

texts: I am using VB.Net to write an add-in. I'd like to keep the tool that I created selected even when the user attempts to select a different tool, and only deselect when my tool is clicked on again.


texts:
I run into errors with distance-based spatial weights in pysal. When running a spatial lag model (e.g. twosls_sp.GM_Lag), I get the following exception:

Exception: All entries on diagonal must equal 0.

Trying the example here, the diagonal is 1 (i.e., an observation is also its own neighbor), and the output is different from that given in the docs:

import pysal
points=[(10, 10), (20, 10), (40, 10), (15, 20), (30, 20), (30, 30)]
w = pysal.DistanceBand(points,threshold=11.2,binary=False)
w.weights[0]
    [0.10000000000000001, inf, 0.089442719099991588]
w.neighbors[0]
    [1, 0, 3]

After row standardization, there are only zeros and nans in the weights matrix:

 w.transform = 'r'
 w.weights[0]
     [0.0, nan, 0.0]

Is this a bug in pysal, or likely a problem with my installation?

I am running the Enthought Canopy Python installation, with pysal installed separately, on a Mac OSX 10.8. Pysal is v1.11.0.


texts:
I'm making the gui for my QGIS3 plugin using PyQT5 and QT designer. I have a very long list of items in a QComboBox and I need the user to be able to select one or more items. I know that QComboBox allows only single selection, so I thought that I can use QListWidget. The problem is that the list of selectable items is very long and I'd need a dropdown menu or something similar.

Is there a way to style a QListWidget like a QComboBox or to allow a multiple selection using QComboBox?

EDIT: I added an immage of my GUI. This is a long list and, as you can see, it's the reason why I can't use QListWidget or QTableWidget. A dropdown checkbox group should be perfect and I'm surprised that QT doesn't make this available.


texts:
I need to convert a .xls file to a point shapefile. The Excel file contains points in Lat, Lon with many attribute fields. I wrote this script to convert from .xls to .dbf and then to shapefile but it creates a .dbf of strings instead of numbers.

How can I fix this?

import os
import glob
import arcpy
from arcpy import env
arcpy.env.overwriteOutput = True

env.workspace = r"E:\test"
os.chdir(r"E:\test")

xlslist = arcpy.ListFiles("*.xls")
for xlsfile in xlslist:
    InFile = xlsfile
    outFile = os.path.splitext(InFile)[0]+".dbf"
    arcpy.ExcelToTable_conversion(InFile, outFile, Sheet="Sheet")

    outlayer = outFile
    outlayer = os.path.splitext(InFile)[0]+"_Layer"
    spatialreference = "GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]];-400 -400 1000000000;-100000 10000;-100000 10000;8.98315284119522E-09;0.001;0.001;IsHighPrecision"
         arcpy.MakeXYEventLayer_management(outFile,"Longitude","Latitude",outlayer,spatialreference,"#")

shpfile = os.path.splitext(xlsfile)[0]
arcpy.CopyFeatures_management(outlayer,shpfile)
del outlayer


texts: Using field calculator to find the next nearest (next minimum value)

texts:
I'm new to Python, but I'm attempting to find the next smallest value in a field, using another field as a case field (within field calculator in ModelBuilder). I've drafted the following, but it does not appear to be functional.

Does anyone have any ideas?

Expression:

NNDist( !ToBreak!, !FID_LicTypesZIP!, !MIN_ToBreak!, !NNDist!)

Code Block:

def NNDist(ToBreak, Customer, Minimum, NNDist):
    for Customer:
        for ToBreak:
            if ToBreak>Minimum and ToBreak<NNDist
                return ToBreak
            else:
                return -1


texts: Linearly comparing one point to all points within 20' foot area, for every point in dataset using ArcPy?

texts:
I have a very large dataset of GPS point that I am trying to write a python script for that will linearly compare one points to the ones around it. Original I was using the ID to determine which point to compare to the next one but do to how the points was collected the ID are not always in order, results in a lot of extra points being selected.

I believe the best solution would be to compare one point to all points within a 20’ foot buffer area of it. I am not sure how to best write a code that does this, I was thinking of creating an in-memory buffer of all points and then using that to select the points that I want to compare to each other.

Here is part of my original script that I am using now (this not using a buffer but an ID# to compare points):

JONum = arcpy.GetParameterAsText(0)
JO = "{}" "{}" "{}" "{}".format('JONumber = ',"'",JONum,"'")
OutPut = arcpy.GetParameterAsText(1)
arcpy.MakeFeatureLayer_management('GPSPipePoint', 'Temp', JO, OutPut)
arcpy.AddMessage('MN Feature Layer')

layer = 'GPS_GasMain_Temp'
fields =  ['ID', 'NAME', 'SIGNAL', 'DEPTH', 'Z_COORD','OPERATORNAME','COLLECTIONDATE']
all_records = [i for i in arcpy.da.SearchCursor(layer,fields)]
all_records = sorted(all_records, key=lambda x: (x[6],x[7],x[1]))

ids_to_select2 = []
for index in range(1,len(all_records)-1): #Skip first and last record
    if (abs(all_records[index][2] >= (1.2 * all_records[index-1][2]))):
        ids_to_select2.append(all_records[index][0]) #If record >+-20%

sql2 = """{0} IN({1})""".format(arcpy.AddFieldDelimiters(datasource=layer, field=fields[0]), ','.join([str(i) for i in ids_to_select2])) #Build a query using ids_to_select list. For example 'GISID' IN(2023,2024)


texts: How to fetch polygons from PostGIS into a PHP array?

texts:
Currently I fetch polygons from a database table using ST_AsText(). The results are in text format, like this:

POLYGON((84 104,212 48,321 122,243 179,275 238,131 240,84 104))

I am not quite sure if ST_AsText() and string manipulation is the way to go for good performance. Is there another way to fetch the points?

How can I get all points of a polygon as an array in PHP?


texts:
Check validity in QGIS and Check geometry in ArcGIS Desktop 10.0 give me different results on the following shapefile

ArcGIS does not return an error at the following geometry:

but QGIS does when using check validity:

Is there an official name for geometry errors like this? If they are bad, how should I fix it? I tried QGIS v.clean, manually editing the polygon and ArcGIS repair geometry without success.

BTW the polygon depicted is GU 14917 (Vertex 397)


texts: How to keep same extent for my 10 raster files?

texts:
I am not able to keep same extent for my raster files and stack them together.
library(raster)
# for reproducible example
r1 <- raster(xmn=72.23081, xmx=77.79537, ymn=34.26277, ymx=37.35659,res=1)
values(r1) <- 1:ncell(r1)
r2 <- raster(xmn=72.51037, xmx=77.68996, ymn=34.51186, ymx=37.098,res=1)
values(r2) <- 1:ncell(r2)
r3 <- raster(xmn=72.2514, xmx=77.9430, ymn=34.20765, ymx=37.40778,res=1)
values(r3) <- 1:ncell(r3)

rasters <- list(r1, r2, r3)
s <- stack(rasters)
# Error in compareRaster(x) : different extent


texts: Saving symbology in GeoPackage with ArcMap?

texts:
I'm working with a GeoPackage file containing polygons in ArcMap 10.6.2
I have hundreds of different classes symbolized one by one that I have stored in a LYR file.

It could be possible to save all this symbology inside the GeoPackage file as I can do with QGIS 3.4 (in style_layer table)?

It is not possible to read the GeoPackage symbology saved with QGIS from ArcMap. So, I'm trying to save this symbology in the GeoPackage with ArcMap.

Any other ideas or solutions?


texts: How to draw polygon by using different point in different UTM zone

texts: I am going to calculate the area of agricultural land the lays in different UTM zones (36 and 37) by using different GPS collected data.  I have GPS XY coordinates and I want to use these to draw a polygon. How?


texts: Specify a GDB feature class using Python Tkinter dialogue box

texts:
I tried to specify a feature class from a GDB dataset through Python Tkinter tkFileDialog.askopenfilename. Below is the script. It worked for me in the past when I was dealing with .txt or .las files. However, this time it is not giving me the GDB feature class, which has no file extension here. Please advise.

import arcpy
from arcpy import env

import Tkinter,tkFileDialog
from Tkinter import *
import tkMessageBox
import ttk
from tkFileDialog import askopenfilename
from tkFileDialog import askopenfile

env.workspace = "C:\\temp\\miniPilot\\miniPilot.gdb"

def browseFor_inFC():
    inFC_name = tkFileDialog.askopenfilename(parent=root, initialdir=r'C:\temp\miniPilot\miniPilot.gdb\data', \
                                                   title='Choose a feature class')
    inFC.set(inFC_name)

def cancel():
    root.destroy()

def remove():
    fields = arcpy.ListFields(inFC)

    for field in fields:
        print field.name
        if field.name <> "OBJECTID" and field.name <> "SHAPE" and field.name <> "SHAPE_Length":

            arcpy.DeleteField_management(inFC, field.name)

root=Tk()
root.title("XXXXXX")

# The first frame named 'mainframe'
mainframe = Frame(root, relief='flat', borderwidth=6)
mainframe.grid(column=0, row=0, sticky=(N, W, E, S))
mainframe.columnconfigure(0, weight=1)
mainframe.rowconfigure(0, weight=1)

# set up Tkinter variables, which are global variables, for the first gui -- the root window
inFC = StringVar()

# widgets on the mainframe
Label(mainframe, text="Specify the feature class: ").grid(column=1, row=1, sticky=E)
inFC_entry = Entry(mainframe, width=36, textvariable=inFC)
inFC_entry.grid(column=2, row=1, sticky=(W, E))
inFC_browser_Button=Button(mainframe, text="...", height=1, width=3, command=browseFor_inFC).grid(column=3, row=1, sticky=W)

# two button widgets
refract_Button=Button(mainframe, text="OK", command=remove).grid(column=0, row=2, sticky=E)
cancel_Button=Button(mainframe, text="CANCEL", command=cancel).grid(column=4, row=3, sticky=W)

for child in mainframe.winfo_children(): child.grid_configure(padx=5, pady=5)

root.mainloop()   # create an event loop


texts: QGIS plugin that adds a point along a line at a specified distance?

texts:
I'd like to have a polyline layer that has distance from the line origin inserted as marker points along the lines. The distances are stored in the attribute table as a field. The lines have directions.

Is there a plugin that can read the distances from the field, and create points along the lines for each record?

The context is a roads layer that will have maintenance distance markers from certain intersections. There are thousands of records and the road shapes are complex, not straight, so an automated process would be ideal.


texts:
I am using ArcMap 10.8.1.
I have shapefiles that are a series of overlapping polygons of drone imagery footprints. I want to loop through each shapefile, and for every shapefile loop through and rasterize every polygon/row in each of those shapefiles. I wrote some code that I adapted from here (Polygon to raster: Creating 1 raster for each polygon in feature class) where I:

find a list of all the shapefiles that end in "*footprints.shp" in an intermediate STEPS folder
create an output subfolder in STEPS\rasterfootprints for each file I process to contain the rasterized footprints
SearchCursor through and create rasters of each polygon in the shapefile.

The following code executes:
import arcpy
import math
import os

root = r"C:\users\myusername\Desktop\ORTHO_FOOTPRINT_AUTOMATION" #change the filepath here
inputs = root + r"\INPUT_CSVS" #input directory of geotag CSVs
steps = root + r"\STEPS" #location of temporary file storage of intermediate processing steps
rfoots = root + r"\STEPS\rasterfootprints" #subfolder in steps to store the 100's of rasterized individual photo footprints
outputs = root + r"\OUTPUT_SHAPEFILES" #location of final

arcpy.env.workspace = steps
arcpy.env.overwriteOutput = True
shplist = arcpy.ListFiles("*footprints.shp")
for shp in shplist:
    prefix = shp.rsplit('footprints.shp')[0]
    os.makedirs(rfoots + '\\' + prefix)
    outfolder= (rfoots + '\\' + prefix)
    PIDname = "PID" #Sequential PhotoID field
    InFeatures=shp
    with arcpy.da.SearchCursor(shp, PIDname) as sCur:
        for row in sCur:
            out_Raster = "PID{}.tif".format(row[0])
            arcpy.MakeFeatureLayer_management(shp, "Layer")  # make a new layer
            arcpy.PolygonToRaster_conversion(shp, "COUNT", os.path.join(outfolder, out_Raster))
            arcpy.Delete_management("Layer")

But it seems to ignore the search cursor through the field "PID" (type:Long) and produces the same raster of the entire area (i.e. a raster of ALL photo footprints) for every file that should be of a single photo.

When it should be giving me something like this (created from running manually in QGIS):

This is what my attribute table looks like:

What do I need to change in my script?


texts: Creating image for google map tiles from KML file

texts: I am totally new to this. I have a KML file of zip code boundaries. I want to create image from this KML, so that can use to overlay google maps. Can someone please guide me here?


texts: Counting overlapping polygons in PostGIS using ST_Union very slow?

texts:
I have a set of 1000 polygons that overlap and I want to count the areas of overlap. I am using this blog post's instructions to create multilines from the 1000 polygons then use this line file to generate polygons for the overlapping areas and then count the overlaps.
http://boundlessgeo.com/2014/10/postgis-training-creating-overlays/

This works fine for datasets of 100 or fewer polygons but hangs on anything more than that on the first step. It seems like the ST_Union is what's so slow, but ST_Collect seems not to work here because it generates a multiline with far fewer nodes and doesn't create non-overlapping polygons when used in the second step. It's too bad because ST_Collect is so fast - but doesn't give the same result. Anyone have ideas on how I can modify this process to work faster on a dataset of 1000 polygons?

Step 1

CREATE TABLE boundaries_polygons1000 AS
SELECT ST_Union(ST_ExteriorRing(wkb_geometry)) AS geom
FROM polygons1000;

Step 2

CREATE SEQUENCE polyseq_polygons1000;
CREATE TABLE polys_polygons1000 AS
SELECT nextval('polyseq_polygons1000') AS id, (ST_Dump(ST_Polygonize(geom))).geom AS geom
FROM boundaries_polygons1000;

Step 3

ALTER TABLE polys_polygons1000 ADD COLUMN count INTEGER DEFAULT 0;
UPDATE polys_polygons1000 set count = p.count
FROM (
  SELECT count(*) AS count, p.id AS id
  FROM polys_polygons1000 p
  JOIN polygons1000 c
  ON ST_Contains(c.wkb_geometry, ST_PointOnSurface(p.geom))
  GROUP BY p.id
) AS p
WHERE p.id = polys_polygons1000.id;


texts: Combining different geometry types in ArcGIS for Desktop?

texts:
Two geodatabases have been created using split by attribute tool from CAD drawing file. One gdb is from polyline layer and another from polygon layer of input drawing file.

I want to combine contours which are polygons in one feature class and contours as polylines in another feature class.
Or
Is there a tool to import similar layers based on 'Layer' filed of CAD drawing irrespective of geometry type into a new output gdb ?


texts: How to add a polygon in Google Earth on a Website

texts:
As it is said in the title i have to add a polygon in Google Earth which is integrated in a website.
I managed to integrate Google Earth in the Website, so i can view the world but i can not use any functions like searching specific places or add paths. Does anybody know what i could do to fix this? Or did i used a wrong method to integrate Google Earth?
I hope this is the right question for this forum.

Ps


texts: Given a collection of IPoints is there a quick way to get the extent that contains them?

texts: I have a bunch of IPoints.  I need an extent which contains all of them.  Is there a quick way to do it?   I had the idea of converting them to screen coordinates and then figuring it out manually but that seems like the wrong way to do it.


texts: Draw polygon from sdo.geometry with OpenLayers

texts:
I have some polygons in my oracle database, the column type is sdo.geometry, and I want to draw the polygon using OpenLayers.

The problem is OpenLayers reads WKT type and I can't convert to this type, because my version of Oracle is "Oracle express Edition", and function like: SDO_UTIL.TO_WKTGEOMETRY(geometry) or geometry.get_WKT() returns "(null)"

.

I already have something on my database cuss:

I'm using C#, and I haven't found any library for my problem.


texts: Spatial join including attributes - like SQL right join?

texts:
I've run in to a little problem. I have two sets of polygons with important attributes. I need to join over the attributes from the first set to the other, but I can't seem to figure out how. Theoretically it seems like an SQL right join, I think. Can you help?

The only way I can get this to work, I think, is the somehow copy the blue polygon. The first problem is that not all of the red polygons are completely inside the blue, but I have selected by centroid (I need to do this). The second problem is that there could be 1, 2, 3 or more red polygons per blue.

The example below would be Blue1+Red1 and Blue1+Red2 and so on.

I have ArcEditor 10.1, QGis and FME.


texts: Using Google Maps static tiles with Leaflet?

texts:
I am trying to use tiles generated with the Google static API with Leaflet and I am facing issues because the generated tile seems to be displayed "off-centred".

My original question has been posted on Stack Overflow.

I am asking the question here in case this site is more suitable (I will also move the question here if the answer is given in this site).


texts:
I am working on Shreve stream order. I made raster where is basin with rivers. I created Shreve stream order but problem is that I should put labels to map, but I don't know how (only putting text manually).

So I tried create vector from that raster and then use NVSVectorToolsSet, but I can't find this toolbox on my ArcMap 10.3 neither 10.1. It is any other way how put that labels from raster to map. It is quite big basin so I don't want put labels there manually if there is a better and faster way.


texts: Using an SQL Query to skip or run a code in ArcPy/ Python

texts:
Using ArcMap 10.4 and Python I am trying to run a code (arcpy.ExportReport) on a bunch of different frequency tables. Which I have completed the only problem is that if I have a table that does not match my definition query the code does not run and will not export the report.

Basically, how do I skip the code from running when the definition query parameters are not met.

Code Block:

            if:
               report_definition_query =SQL_NAME2arcpy.mapping.ExportReport(lyr,Template,PDF,"DEFINITION_QUERY",report_definitonn_query =SQL_NAME2)

           else:print "Blank Report"


texts: Add File/personal Geodatabse to Geoserver

texts:
Is it possible to add file or personal geodatabase to Geoserver?
If not, how we can get rid of shapefile limitations in Geoserver (without using a database like PostGIS)?


texts: Looking for a data set of geo tracking for a group of users

texts:
I run a study on user commuting patterns and I need a data set where each line would be something like:

userID | Timestamp | lat | long

Then, I would like to analyze each user separately and derive insights on his/her commuting patterns. I understand that this kind of data set is rather tricky in terms of privacy, but maybe someone here has an idea?

So far I looked into the following sources:

OSM but here the data points do not contain user IDs (for clear reasons)
Uber in this case, the "user" is a taxi driver and his/her commuting patterns are not interesting for my purposes.
GeoLife looks promising but I fail to download the data set

Any ideas? Sources?


texts: Issue with Sextante r.map calculator and modeler

texts:
I am trying to build a model with Sextante modeler and the r.mapcalculator algorithm but when I run the model I get the error message:

‘Could not load layer: ‘….’Check sextant log to look for errors in algorithm execution.’

The Sextant Log reads:

sextante.runalg("modeler:ndsi.","H:/NDSI/MODIS/20120607_SurfaceReflectance
   /0607_Processed Images/0607_b4_warp","H:/NDSI/MODIS/20120607_SurfaceReflectance
      /0607_Processed Images/0607_band6_warp","H:/NDSI/MODIS/test/NDSI")

The algorithm I have put into r.mapcalculator is a basic algorithm to calculate NDSI for snowcover analysis:

e.g. (amap-bmap)/(amap+bmap)

amap = band 4

bmap= band 6

I have tried on different computers as well as only using the r.mapcalculator through the Sextante Toolbox. Is there a different calculator I can use in my model? Is there a different way I should write out my algorithm?


texts: I have omi (aura satellite) level 2 data sets for the aerosol product in he5 file format. Now I want to read that he5 file and save it as geotiff format in python. Could you help me anyone regarding this issue?


texts: I wonder if I can do simple arithmetics with coordinates. Can I do substraction of coordinates like 2.570801 (lng) 46.89023 (lat) / 17.424316 (lng) 55.15377 (lat)?! Or will that get messy without a CRS?!


texts: Calculating many shortest paths with pgRouting

texts:
I 145013 Danish addresses, and 10 big cities in Denmark as points in QGIS. Furthermore, I have four vector layers that connect all the points. They are small roads, main roads, freeways and ferry routes. See the screen dump from QGIS below for an illustration:

Now I would like to calculate the shortest distance from each address to each of the 10 big cities. I will use pgRouting for this, and I have followed the nice tutorial on finding 1 shortest path here http://anitagraser.com/2013/07/06/pgrouting-2-0-for-windows-quick-guide/

However, I can't figure out how to extend the problem from 1 path to (145013 * 10 =) 1450130 paths.

What method should I use in pgRouting, and how do I incorporate all the 1450130 paths?

I use:

QGIS 2.1.0
postgres 9.2
postgis 2.0
pgRouting 2.0
OS: Win 7 64bit

I have shared my data here: https://www.dropbox.com/sh/rn8gdxr1bki4q4x/fFRd0PerkP

It's in in a QGIS project file.

I appreciate all comments and answers!


texts: How do I make OpenLayers.Class.create() work in Openlayers 2.12?

texts:
At version 2.12 some functionality has been changed and moved to the "deprecated.js"
I have the following code which I want to be able to run without "deprecated.js". how can I convert it ?

OpenLayers.Layer.MapInfoRESTHandler = OpenLayers.Class.create();
OpenLayers.Layer.MapInfoRESTHandler.prototype =
    OpenLayers.Class.inherit( OpenLayers.Layer.Grid, {...

I know OpenLayers.Class() is funtion now , but I'm not sure how to use it in current syntax so it will create the same class as the code above (and same with inherit)