texts: SRID of geoemtry exported in PostgreSQL with ogr2ogr

texts:
I use ogr2ogr to export geometric data from EDIGEO files to some PostgreSQL table in append mode, but the output geometry SRID is always set to zero.

Can I do something to specify the output geometry SRID ?

the command I use is :

ogr2ogr -append -dialect SQLITE -f Postgresql PG:"dbname=test
port='5434' user='postgres' active_schema=napo"  thefile.thf
-sql "SELECT cast(commune_id.idu as smallint) AS insee ,section_id.tex AS section,
parcelle_id.tex AS numero,setsrid(parcelle_id.geometry,3942) AS geom,
setsrid(parcelle_id_label.geometry,3942) AS point
FROM commune_id
LEFT JOIN section_id ON section_id.idu = section_id.idu
LEFT JOIN parcelle_id ON parcelle_id.idu = parcelle_id.idu left join
LEFT JOIN parcelle_id_label ON parcelle_id.object_rid = parcelle_id_label.ogr_obj_lnk"
-nln napo -a_srs EPSG:3942 -s_srs EPSG:3942

The table is created by the first call of the command.


texts: Delineate watershed using manual points in RSAGA

texts:
I'm currently trying for an automated solution to delineate watersheds within a large basin. I have a collection of locations (X,Y coords) for various tributaries reaching a larger river, and the end-game is to have each sub-basin delineated so I can calculate areas etc. I have been working in R, using SAGA and the RSAGA packages, (R 3.2.2, SAGA 2.1.2, OSX Yosemite 10.10.5).

So far, I have been able to fill the sinks:

rsaga.fill.sinks('dem_saga.sgrd', 'dem_sinkfilled.sgrd', method = "xxl.wang.liu.2006")

And have been able to make a raster of the same size and resolution as the dem, where all cells are 'NA' except for a single cell at the desired location to delineate the watershed (target.sgrd). This was achieved using rasterize command in R (e.g. rasterize(target_point, dem, 1)).

Unfortunately, although I am pretty sure the command I want is the 'Upslope Area' module, when it runs, I end up with a raster where every cell = 100.

I'm not sure exactly how to make a reproducible example, but my code:

rsaga.geoprocessor(lib = 'ta_hydrology', 4,
                    param = list(TARGET = 'target.sgrd',
                                 ELEVATION = 'dem_sinkfilled.sgrd',
                                 AREA = 'target_catch.sgrd',
                                 METHOD = 0))

Finally, the tutorial I have been working through is here. Jumping down to "Creating a Watershed map".


texts: Using one map, multiple layouts, each with different layers in ArcGIS Pro?

texts:
I was wondering whether it's possible in ArcGIS Pro to create multiple layouts for one map, which contains different layers (e.g. aerial photo, buildings and soiltype).

After prepping the map in 'map view', I'd add a new layout for e.g. a map of the building, and add the map frame. In the contents view I'm toggling layers on and off for what I need in that particular layout (e.g. aerial photo + buildings). Then, I'd add another new layout for the soil type, toggling layers on and off in contents (e.g. aerial photo + soiltype). Unfortunately, it seems that toggling layers on and off in contents view affects all layouts.

Is there any way to toggle layers on and off within a single layout view of the same map, without affecting the other layouts?

I just can't seem to figure this out but there must be something I'm not grasping about ArcGIS Pro, because this seems like pretty basic functionality.

I'm using ArcGIS Pro 1.4.


texts: Is OpenLayers appropriate for creating maps of a fantasy world?

texts:
Please excuse my ignorance, I've been directed here from Stack Overflow where I've been looking at a variety of HTML5/JQuery packages to create a map of a fantasy world that exists in a game I maintain a player fansite for.

The features I desire are

Free to use
Requires no server side installations
Allows for markers/annotations
Supports layers (for example roads)
Able to be zoomed/scaled
Able to switch from a "global" view to a more "local" map (eg from continents down to cities and streets)
Requires no knowledge of GIS by the end user beyond the sort of experience Google Maps offers
Supports some sort of overlay when a region is moused over/into

I did not consider open layers because I thought it might be too slow or overkill for my project, but as I'm looking at a variety of jquery solutions which only work together questionably well I'm now willing to consider a more substantial learning investment if it will help me achieve my goal.


texts:
I'm sure this question is really really basic but I appeal to you because I hate the idea that I have spent 20 minutes trying to do such a simple thing...

All I want to do is style a polygon layer by increasing complexity of a pattern.

ex. low value polygons will have very spaced out lines, high values will have lots of lines very close to each other...

Again, I know this sounds incredibly low level, but I need some help.


texts:
Does anyone know if grid control used within the batch mode of ArcToolbox available via ArcObjects (see image 1)?

I have used the IGPValueTable control (as used in the Union tool, image 2), however, this does not appear to allow you to add domains such as the field type drop down comboBox shown in the attached image.


texts: Selecting features by location in QGIS (boundaries_touches)

texts:
Is there a tool or function in QGIS which can select features of one layer which are located on the boundaries of another layers features? In ArcGIS it's the "select features by location" tool and the relationship set to "boundaries_touches", but I need to work it into a QGIS Model for accessibility.
Edit:
The shapefile for the selection is a grid and looks like this:

The selection reference is the grid dissolved into Basins and looks like this:

The operation should select every feature/grid cell which touches the borders of the features of the basin shapefile.


texts: Height and X Y coordinates of points on Google Earth

texts: How can I get information on height and X,Y coordinates using Google Earth data for a given terrain?


texts: How to reproject data in multiple UTM zones? (DRCongo, Rwanda and Uganda)

texts:
I am working in the area near to Virunga National Park, in DRCongo, Uganda and Rwanda.

I have established my working area in the zone 35 south. However, I also have to work with data from Uganda and Rwanda, which is not only in that zone. The data from Uganda and Rwanda is in GCS WGS84.

Should I first project my data to UTM 35S, even though not all the area falls into zone 35S?
and then,

Should I reproject all of my data to South Pole Lambert Azimuthal Equal Area or Albers equal-area? Which is better to my working area?

Thanks for your help!


texts:
Does someone have an example of a use of the leaflet-geotiff plugin? I have difficulties. Even the demo in the files don't work for me :
https://github.com/danwild/leaflet-geotiff-2
(the raster doesn't appear)
It works on this website but I don't know why : https://danwild.github.io/leaflet-geotiff-2/


texts: arcpy.da.SearchCursor feature orders mixed up in Oracle enterprise geodatabase?

texts:
Something very weird happen when I use a SearchCursor over a feature class in an Oracle Enterprise Geodatabase

The order of the features in the SeachCursors changes depending of the field_names parameter that is used.

For instance:

source = "\\\\cnatrtd8\\geo\\ArcGIS Server\\Connections\\GEODEV011.sde\\GEO09E01_CS_FRA_GEN"
with arcpy.da.SearchCursor(source, ["OBJECTID"]) as cursor:
    ids = [row[0] for row in cursor]
    print(ids)

with arcpy.da.SearchCursor(source, ["OBJECTID", "CD_CS_FRA"]) as cursor:
    ids = [row[0] for row in cursor]
    print(ids)

>>>
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60]
[39, 40, 41, 42, 43, 44, 45, 46, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 47, 48, 49, 50, 51, 52, 53, 55, 56, 54, 57, 58, 59, 60]
>>>

This happen with many different features class, but only in the Enterprise Geodabase. The same features class in a File geodatabase dont do this.

What could be the cause of that?


texts: Name for meridians rotated 90 degrees about a meridian

texts:
Let's say you drew lines similar to meridians, in that they all intersect at two opposing points--but instead of being at the poles, those points were instead somewhere along the equator.

Is there a term for this?


texts: What is the method used by wikipedia to solve Trilateration problem, is it the Linear Least Square? or is it just a direct solution of a set of linear equation systeme with exact distance? but in this case how this method can find solution even when the three circles does not intersect at one point?


texts: coordinate system for calculating distance at the nation scale

texts:
I want to use ArcGIS 10.1 to calculate the straight (euclidean) distances from zip codes to their nearest city centers in about 300 U.S. metropolitan areas. I have two point shpfiles: one includes the centroids of all zip codes and the other includes the points that represent city centers. Both of the two shpfiles are in the GCS_North_American_1983, which prohibits me from doing the distance calculation.

I have converted the GCS into North America Lambert Conformal Conic to do the calculation. It worked. But I am not sure whether this is the best way to do it. Any advice is appreciated!


texts:
I'm trying to export vector files from QGIS to KML.  When I open the KML in Google Earth I get erratic lines all around the globe.  All layers set to WGS84.  I have read the other topics on KML export and tried the suggestions with no success.  Seem so simple but something is off.

Here is the boundary I'm trying to export.  It is a line not poly.

This is what i get when I open the Exported KML to GE


texts: How to create a map showing NDVI result in erdas?

texts:
I'm using Spot 5 image to run an NDVI analysis. The band that I used in processing is 4,3,2. The area was 80% covered by mangrove.
My problem here how can I create a map to show or present the result?
Basically I run the NDVI using the function automatically under unsupervised > ndvi in erdas.


texts: UnsafeLoadingWarning when updating anaconda environment Arcgis Pro 2.5

texts:
I updated arcgisPro to 2.5.1 and it had a version mismatch with my arcpy environment, so I did the following to create a clone.

To create a new environment that is up to date with arcgispro-py3:
 1. Generate a list of additional packages installed in your current environment,
    with the conda command:

  conda env export > env.yaml
   2. (Optional) If you have additional dependencies installed through pip,
      find those with pip freeze:
  pip freeze > requirements.txt
   2. Create a new environment by cloning arcgispro-py3, and activate it:
  conda create --clone arcgispro-py3 --name my-env --pinned
  activate my-env
   3. Add back missing conda packages from your current environment:
  conda env update -n my-env -f env.yaml
   4. (Optional) Add back missing packages from pip:
  pip install -r requirements.txt
   5. (Optional) Make it the default for the ArcGIS Pro application and the
       "Python Command Prompt":
  proswap my-env

However, at step 3,  I received the message:

C:\Users\USERNAME\AppData\Local\Temp\_MEI79962\conda_env\env.py:68: UnsafeLoaderWarning:
The default 'Loader' for 'load(stream)' without further arguments can be unsafe.
Use 'load(stream, Loader=ruamel_yaml.Loader)' explicitly if that is OK.
Alternatively include the following in your code:

The file env.py (and the folder _MEI79962) appear to be temporary, so I am not sure what script I should edit to make this work. Has anyone solved this problem?


texts: Creating new feature classes from rows in polygon feature class?

texts:
I'm trying really hard here to get an arcpy equivalent for ARCGIS' ModelBuilder's Iterate Feature Selection. I need to export the individual polygon rows in a feature class as new feature classes.

Any help?


texts: ArcGIS Pro Symbology: Polygon boundary is displayed at the wrong size

texts:
This feels like a silly question, but ArcGIS Pro is displaying the boundaries of my polygons as very wide. Any suggestions on what settings to change to put it back to normal?

Image 1 below is the current display, with a 0.7 outline boundary. Image 2 is no boundary.


texts: How do I change the snapping sessions for a layer in QGIS 3

texts:
In QGIS 2 I would go the the settings menu and select the snapping options but this option is not there in version 3.  I have found the snapping options in digitising options but these don't seem to allow one to turn snapping off and on on a layer by layer basis.

I have done standard searches but can't find anything for GIS 3 there may be stuff there but it must be buried in all the hits for older versions.

The problem I am having is that I can't move vertexes small distances and changing snapping tolerance does not seem to make any difference.  This was very confusing when trying to come to grips with the new vertex tool.

Update: I have also found Project->Snapping options which again does not allow control on a layer by layer basis but did allow me to set the snapping radius in meters which helped which has relived the immediate pain but I still need to be able to turn snapping off for this layer.


texts:
I am using ArcGIS 10.5 and i am trying to get metadata information of a featureclass by calling following statement from arcpy within the ArcGIS Environment.

import arcpy
from arcpy import metadata as md

But i am getting following error

Runtime error
Traceback (most recent call last):
  File "<string>", line 1, in <module>
ImportError: cannot import name metadata

Do i need to enable some extension or import this library some other way to make this work ?


texts: Why the WFS point is not shown to the correct location after a WFS Update request

texts:
I've been working on the implementation of an edit the position of a point functionality.
I use leaflet and geoserver and give the possibility to the user to change the location of a point.
It works as expected but its a bit slow. As you can see in the image, although the WMS point has moved correctly, the WFS point is still in the old position. It takes some time (approx. one minute) in order to be displayed correctly.

Any ideas why this is happening?

EDIT

This is how I implement the edit. I actually send a POST request using Python Requests library. The XML of the POST request looks like this:

<wfs:Transaction service="WFS" version="1.1.0"
  xmlns:ogc="http://www.opengis.net/ogc"
  xmlns:wfs="http://www.opengis.net/wfs"
  xmlns:gml="http://www.opengis.net/gml">
 <wfs:Update typeName="geonode:{{ layer_name }}">
  <wfs:Property>
   <wfs:Name>the_geom</wfs:Name>
    <wfs:Value>
     <gml:Point>
       <gml:pos>{{ coords }}</gml:pos>
     </gml:Point>
    </wfs:Value>
  </wfs:Property>
 <ogc:Filter>
   <ogc:FeatureId fid="{{ coords }}"/>
  </ogc:Filter>
 </wfs:Update>
</wfs:Transaction>


texts: Creating water surface from vector line in QGIS

texts:
I want to create a raster surface representing the flood level in a creek based on output from a 1D hydraulic model.

I receive a 3D profile of the water surface in a creek line as a 3D Polyline (AutoCad format) from our hydraulic modelers.  Each vertex of the polyline represents the coordinates and level of the water surface based on cross-sections used in the modelling.

I want to import the string into QGIS and generate a water surface as a raster for analysis - for example depths, extents, comparison of flood level to building floor levels, etc).

Basically I want to extend a surface horizontally and perpendicular from the string until it intersects a nominated raster surface (existing ground), along the length of the string.  I would assume that the elevation varies linearly between each vertex of the polyline (which would be accurate enough for my work).

Is there a way to generate this surface in QGIS?  I would appreciate any guidance.

Thanks

PS I have edited this question to hopefully make it clearer.


texts: Adding all the layers in a KML file in PyQGIS

texts:
I've a KML file that contains several layers, if I use the usual method of adding vector layers from a file it only adds the first layer (the Point layer)

vlayer = QgsVectorLayer("C:\Users\HOME\Desktop\Test.kml", "Test", "ogr")

QgsMapLayerRegistry.instance().addMapLayer(vlayer)

How to add all the layers in this file?


texts: Is there an online source of raw contour data?

texts:
My goal is to get a vector image (SVG preferably) of a local mountain. The data clearly exists, but as an armchair cartographer without access to raw layer data I've been reduced to combing the web for sources of just contours.

I usually find it mixed with relief shading and land ownership boundaries which make the image recognition (live tracing in Illustrator, for example) quite difficult.

Has anyone found a mapping site with a web interface that provides such a thing?


texts: cartodb: query for the points within a radio from another point - postgis

texts:
I've georreferenced a couple of points in a CartoDb table, so now I have a the_geom column of type point

I'd like to get all the points that are x meters away from a certain x,y point.

So far now, following this article http://unserializableone.blogspot.com.ar/2007/02/using-postgis-to-find-points-of.html, I tried the following:

SELECT * FROM my_contacts
WHERE
distance(
  transform(PointFromText('POINT(-34.6043183 -58.380924)', 4269),32661),
  the_geom
) < 1000

(in this case my x,y point would be -34.6043183 -58.380924 and I'm looking for locations within 1000 meters from it)

and I get this error:
    ERROR: Operation on two GEOMETRIES with different SRIDs

I tried without the transform, and I get the same error

How can I find out the SRIDs of the points in my db?

And how can I translate a lat, lon google map point to that desired SRID?

update:

with this query I found out the SRID of my points

SELECT ST_SRID(the_geom) FROM my_contacts

they are all 4326, so I tried with

and this seems to work rather ok:

SELECT *
FROM my_contacts
where
distance(
  transform(PointFromText(
    'POINT(-34.6675645 -58.3712721)', 4326),4326),
  the_geom
) < 33.62
order by distance

but I just got thiw 33.62 by trial and error, I don't know how to translate it to meters, and the results when trying with a different point don't seem to be very consistent...

thanks a lot


texts: Parsing a raw string to select by expression

texts:
This works as expected:

lyr = iface.activeLayer()
lyr.selectByExpression("\"fieldname\" > 0")

According to this this article about parsing raw string in python the following lines of my code should work just the same but they don't.

query_dict = {"query1":r"\"fieldname\" > 0"}
lyr.selectByExpression(query_dict["query1"])

What am I missing out on?


texts: Changing the line width/weight in Leafletr addPolylines function depending on an attribute

texts:
I would like to display flow lines (OD-lines) using leaflet in Shiny and sofar I can successfully map all the lines I want. But I also want to change the line width so that each line will have a different linewidth depending on the attribute they have. I am trying to get the values for weight with the following code but it does not work:

  values <- reactive({
    test[quel==input$zelle&hr==input$std&vmitt==input$vmitt&tag==input$tag]
  })

  observe({
    anzhl <-isolate(values()[,10])
    print(anzhl[1]*2)
    str(anzhl)

    leafletProxy("zones", data = values()) %>%
    removeShape(layerId="liID") %>%
    addPolylines(lng=~x_coor, lat=~y_coor, stroke=TRUE, layerId="liID", color="black", weight=sum(anzhl())/100)})

I get the following error message:

Warning: Error in observerFunc: could not find function "anzhl"
Stack trace (innermost first):
    56: observerFunc [E:\.../server.R#118]
     1: runApp
ERROR: [on_request_read] connection reset by peer

What am I doing wrong and is it ever possible to change the line widths with the values called within an interactive session.


texts: Problem loading from MySQL in QGIS 2.10 on Mac

texts:
I cannot load a vector layer from MySQL in QGIS 2.10 (Pisa) on a Mac. If I goto Layer -> Add Layer -> Add Vector Layer -> Database I do not see MySQL in the drop down.

I found this link. I did not go through all the steps in that answer. Rather I jumped to step 6 where it says to run ogrinfo --formats. This is what I got:

$ ogrinfo
dyld: Library not loaded: /usr/local/opt/mysql/lib/libmysqlclient.18.dylib
  Referenced from: /usr/local/bin/ogrinfo
  Reason: image not found
Trace/BPT trap: 5

After some mucking around I found that I already have a similar file at:

$ ls /usr/local/mysql/lib
libmysqlclient.20.dylib libmysqld-debug.a   mecab
libmysqlclient.a    libmysqld.a     pkgconfig
libmysqlclient.dylib    libmysqlservices.a  plugin

I had installed MySQL from a DMG file if that helps. I copied all these files recursively to /usr/local/opt/mysql/lib and also renamed libmysqlclient.20.dylib to linmysqlclient.18.lib to end up with following:

$ ls /usr/local/opt/mysql/lib
libmysqlclient.18.dylib libmysqld-debug.a   mecab
libmysqlclient.a    libmysqld.a     pkgconfig
libmysqlclient.dylib    libmysqlservices.a  plugin

Now I get:

$ ogrinfo --formats | grep 'MySQL'
  -> "MySQL" (read/write)

But Step 7 of the answer does not work as expected:

>>> from osgeo import ogr
>>> sorted([ogr.GetDriver(drv).name for drv in range(0,ogr.GetDriverCount())])
['ARCGEN', 'AVCBin', 'AVCE00', 'AeronavFAA', 'BNA', 'CSV', 'CartoDB', 'CouchDB', 'DGN', 'DODS', 'DXF', 'EDIGEO', 'ESRI Shapefile', 'ElasticSearch', 'GFT', 'GME', 'GML', 'GMT', 'GPKG', 'GPSBabel', 'GPSTrackMaker', 'GPX', 'GeoJSON', 'GeoRSS', 'Geoconcept', 'Geomedia', 'HTF', 'Idrisi', 'Interlis 1', 'Interlis 2', 'KML', 'LIBKML', 'MSSQLSpatial', 'MapInfo File', 'Memory', 'NAS', 'ODBC', 'ODS', 'OGDI', 'OSM', 'OpenAir', 'OpenFileGDB', 'PCIDSK', 'PDS', 'PGDump', 'PGeo', 'PostgreSQL', 'REC', 'S57', 'SDTS', 'SEGUKOOA', 'SEGY', 'SOSI', 'SQLite', 'SUA', 'SVG', 'SXF', 'TIGER', 'UK .NTF', 'VFK', 'VRT', 'WAsP', 'WFS', 'Walk', 'XLS', 'XLSX', 'XPlane']

MySQL is not listed in above. What am I missing here? Tried running QGIS again and don't see MySQL in the drop down. Any help?

I have also tried copying mysqlclient.20.lib to /Library/Application\ Support/GDAL/1.11/PlugIns/ with no success.


texts: How to use 931 GP service with v10 Map Services?

texts:
Pinched the code from this profiler sample.
Having some trouble setting it up with my custom Esri Flex application.

My Flex application is using the Web Mercator Aux Sphere WKID and is using v10 map services, and therefore is 3857 WKID.

When using the tool within my application, I get a 500 error.  This is the request is sends:

http://sampleserver2.arcgisonline.com/ArcGIS/rest/services/Elevation/ESRI_Elevation_World/GPServer/ProfileService/execute?f=json&Display_Segments=true&Image_Width=600&Input_Polylines={%22features%22%3A[{%22geometry%22%3A{%22paths%22%3A[[[16616585.72030409%2C-4240972.45511179]%2C[16390943.612806262%2C-4351041.775842438]]]%2C%22spatialReference%22%3A{%22wkid%22%3A3857}}}]}&Image_Height=250

sampleserver2 is on a 9.3.1 setup.

If I swap out the 3857 at the end of that link for the older 102100, I get a correct JSON response.

So how do I get around this?

I can't see an equivalent service on the v10 sample servers.
My environment is all v10, and I would prefer to keep my app within the newer 102100 WKID.

Further info on the change in WKIDs.


texts: Is there an arXiv-like electronic archive used by the remote sensing community?

texts:
Mathematicians, physicists, and (theoretical) computer scientists use arXiv extensively to share their preprints. Indeed, one could say that it's the established way to share one's research.

Is there a similar preprint service that's used by remote sensing, GIS, and related communities?

A related question: Preprint services other than arXiv (for other fields) on Academia.SE.


texts: The ee.Reducer.toList() reducer does not return the same number of features as the feature collection

texts:
I have a feature collection (samplePlots) with 2755 features. When I use the code below, I expect a list with 2755 elements, but the result is a list with only 2515 elements. What would be the reason? Unfortunately I can't share the feature collection.
image.reduceRegion({
  reducer: ee.Reducer.toList(),
  geometry: samplePlots,
  scale: 30
})


texts: Retrieve selected feature on ArcGIS Online as string?

texts: For now I am trying to create a webmap with the ArcGIS API (4.x) for JavaScript. I have a bit understanding of the workings of JavaScript, however I couldn't find code in which it's possible to select a feature (let's say a building) on the web map and return a string when selecting on this feature (such as the name of selected feature. I want to use the string of the name (e.g.) for another action.


texts: Failed to plot lat/long coordinates in decimal degree format using my QGIS 2.14.2

texts:
I am using QGIS 2.14.2 but whenever I try to plot my coordinates(Lat/Long; decimal degrees) from a CSV file, no coordinates are plotted. I even tried converting the coordinates from lat/long to UTM before re-uploading the CSV file but it still failed to plot the points.
Interestingly, the layer's attribute table shows that everything is well captured. BUT there are no points plotted on the display. What do I do?


texts: Get Raster Catalog Names from a File Geodatabase Programmatically

texts:
I'm trying to get only raster catalog entries from a file geodatabase using the arcgisscripting geoprocessing API. So far I am able to do this using the ListDatasets function using the ALL argument. The problem is that it also returns raster datasets as well and I suspect there are other types it would return if they were present in the file geodatabase. I've tried the Raster argument but it returns the raster datasets and omits the raster catalogs. I guess I could filter the items returned from gp.ListDatasets("", "Raster") from the items returned from gp.ListDatasets("", "ALL") leaving me with what I'm looking for, however there might be some other type I don't know of that would be left in the list if it were present.

Here is the code I'm starting with:

import arcgisscripting
gp = arcgisscripting.create(9.3)
gp.workspace = r"C:\FileGeodatabase"
datasets = gp.ListDatasets("", "ALL")
for dataset in dataset:
    #Only want raster catalogs here...

I'm looking for a way to get only raster catalogs. What approach would ensure I'm working only with raster catalogs?


texts: Adding a time/datestamp when creating a video thumbnail URL in Google Earth Engine?

texts:
How can I add a text label showing the date of capture of each image in the animation?
var s2 = ee.ImageCollection("COPERNICUS/S2_SR")
    .filterDate('2020-03-01', '2021-03-20')
    .filterBounds(point)
    .filter(ee.Filter.lt('CLOUD_COVERAGE_ASSESSMENT', 5));

var gifParams = {
  'region': polygon,
  min:0,
  max:2000,
  bands: ['B4', 'B3', 'B2'],
  'dimensions': 600,
  'crs': 'EPSG:3857',
  'framesPerSecond': 2
};

print(s2.getVideoThumbURL(gifParams));


texts: Setting a relative path to load a QML style in QGIS plugin

texts:
I would like my plugin to apply a QML style to an output layer. This consists of three parts:

Packaging the QML file into the plugin.
Referring to the QML file using a relative path.
Setting the style on the output layer.

For number 1, I have saved my style as default.qml in a /styles/ subdirectory of the plugin, and I have added it like so into the resources.qrc file.

<RCC>
  <qresource prefix="/plugins/phylo_tree" >
     <file>styles/default.qml</file>
  </qresource>
</RCC>

I have then used resources.qrc to generate resources.py with pyrcc5

For number 2 - referring to the QML file by a relative path - my understanding is that the path should be :plugins/phylo_tree/styles/default.qml, as mentioned in this answer.

For number 3 - I have used the code in this answer to add a style to the output layer, using the postProcessAlgorithm function. This works if I use an absolute path to the QML file, but not when I use the relative path :plugins/phylo_tree/styles/default.qml: When using the relative path, the operation completes without any error message, but no style is applied.


texts: OSM layer does not print in GeoServer printing module

texts:
I am new with the GeoServer printing module, in GeoServer version 2.6.0, freshly installed a few days ago.

I can now print WMS layers and I think I configured the OSM data correctly. But it still does not work for printing OSM, so I must have mistakes somewhere.

I am using OpenLayers and the default url for OSM seems to be "tile.openstreetmap.org"

I have set up the OSM url in the hosts section in the config.yaml file:

- !dnsMatch
   host: tile.openstreetmap.org

My printing url is constructed as follows:

var url = 'http://xxx.xxx.x.xxx:8080/geoserver/pdf/print.pdf?spec=';

var options = {
    units: 'm',
    srs: 'EPSG:900913',
    layout: 'A4 portrait',
    dpi: '300',
    mapTitle: 'This is the map title',
    comment: 'This is the map comment',
    layers: [
        {
            baseURL: 'http://tile.openstreetmap.org',
            singleTile: false,
            type: 'OSM',
            maxExtent: [-20037508.3392,-20037508.3392,20037508.3392,20037508.3392],
            tileSize: [256, 256],
            extension: 'png',
            resolutions: [156543.0339, 78271.51695, 39135.758475, 19567.8792375, 9783.93961875, 4891.969809375, 2445.9849046875, 1222.99245234375, 611.496226171875, 305.7481130859375, 152.87405654296876, 76.43702827148438, 38.21851413574219, 19.109257067871095, 9.554628533935547, 4.777314266967774, 2.388657133483887, 1.1943285667419434, 0.5971642833709717]
        },
        {
            baseURL: 'http://xxx.xxx.x.xxx:8080/geoserver/abc/wms',
            opacity: 1,
            singleTile: true,
            type: 'WMS',
            layers: [layers],
            format: 'image/jpeg'
        }
    ],
    pages: [
        {
            center: [center.lon, center.lat],
            scale: map.getScale()
        }
    ]
};

var printUrl = url + JSON.stringify(options);

As you can see below, the WMS layers are correctly printed, but the OSM layer does not appear, even though there was no error message during the printing.

http://we.tl/nI6y7lNNhV


texts:
Is it possible that arcpy converts all inputs of the a value registered as GPlinearunit (kilometers, meters, feet, miles, etc) into units of XY coordinates system?
I need to perform calculations with it, since ArcPy GPlinearunit recognizes the parameter as a string.

I have:

import arcpy, math

source = arcpy.GetParameterAsText(0) # < feeature point
w = arcpy.GetParameterAsText(1) # <- Double
d = arcpy.GetParameterAsText(2) # <- GPlinearunit

for row in arcpy.da.SearchCursor(source, ["SHAPE@XY"]):
    x, y = row[0]
    arcpy.AddMessage("{}, {}".format(x, y))

rad = float(w) / 57.2957795
arcpy.AddMessage("{}".format(rad))
dX = x + float(d) * math.cos(rad)
dY = y + float(d) * math.sin(rad)
arcpy.AddMessage("{}, {}".format(dX, dY))

Error message:

  Traceback (most recent call last): File
  "C:\AIR-GIS\toolbox\scripts\Perfil.py", line 26, in  dX = x +
  float(d) * math.cos(rad) ValueError: invalid literal for float(): 5
  Kilometers Failed to execute (Perfil).

I found a solution: convert all units in meters.

digit = int(filter(str.isdigit, str(d)))
arcpy.AddMessage(str(digit))
letters = str(filter(str.isalpha, str(d)))
arcpy.AddMessage(str(letters))

if letters == "Meters":
    dfloat = float(digit)
elif letters == "Kilometers":
    dfloat = float(digit) * 1000
elif letters == "Unknown":
    dfloat = float(digit)
elif letters == "Feet":
    dfloat = float(digit) * 3.28084

rad = float(w) / 57.2957795
dX = x + dfloat * math.cos(rad)
dY = y + dfloat * math.sin(rad)


texts: How I can reverse numbers in a string without reversing any text?

texts:
I have a string, I need to reverse just the numbers in the string without the text:

String: [321 test]

The result will be: [123 test]


texts: Seeking neighborhood boundaries for Saudi Arabia?

texts:
I currently have set of coordinates (lat,lon) that I would like to associate with neighborhoods or districts in Saudi Arabia. Therefore,I might be looking for neighborhood boundaries (polygons coordinates). What would be the best way to achieve this? I do not mind paid services, but of course I would prefer open/free services if possible.

Edit: searching in Google Maps selects the district boundaries, but I couldn't download this polygon using their API.


texts: Does Processing plugin v2.10.1 not recognise SAGA?

texts:
Version 2.9.3 of the Processing plugin works fine, I can run the SAGA tools successfully. Updating the Processing plugin to version 2.10.1 results in SAGA disappearing entirely from the toolbox. The path is fixed at:

C:\Users\Me\.qgis2\python\plugins\processing\algs\saga

Has anyone else come across this issue?

I'm using QGIS 2.8.2 on Windows 7 64-bit.


texts:
I am working with ArcGIS 10.2 and using python to use the function for zoom to selected feature.  I need to have this zoom to selected feature work for any project and I also need drop down menus for all of the SQL's (Quarter section, township, range, and meridian) that I currently have.

#Zoom to Selected Feature in SearchMap Dataframe

import arcpy
mxd = arcpy.mapping.MapDocument('CURRENT')
df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
df.zoomToSelectedFeatures()
df.scale = df.scale * 1.1

# Clear selections on the Southern Reference Grid layer DataFrame

clearLyrs = ["Southern Reference Grid"]
for clearLyr in clearLyrs:
    lyr = arcpy.mapping.ListLayers(mxd, clearLyr, df)[0]
    arcpy.SelectLayerByAttribute_management(lyr, "CLEAR_SELECTION")

arcpy.RefreshActiveView()


texts: How can I make a PostGIS ST_Distance_sphere query more efficient?

texts:
Within PostgreSQL/PostGIS I have written a SQL query to get the closest distance for each element within a point coordinate table from a polyline table.  The query comes back and works correctly, but it takes about 1 minute per record (see SQL statement below note: I have limited it to the first 10 records so that I can derive how long it would take per record otherwise it would take days to process all 366,681 records). I have created a GIST index on both the "the_geom" table fields.

The other table which is polyline has 369,962 records.

How can I make a PostgreSQL/PostGIS ST_Distance_sphere query more efficient?

SELECT "ParAdd", "ParZip", min(distance) FROM
  (SELECT a."ParAdd", a."ParZip", ST_Distance_sphere(a.the_geom, b.the_geom) distance FROM
    (SELECT "ParAdd", "ParZip", the_geom
      FROM reflux_blue_06_21_2012_15147_va_md_premise_res_prosp_filtered_t limit 10
    ) a,
    ( SELECT the_geom FROM reflux_blue_06_21_2012_15147_va_md_premise_gis_main_pipe_4326
    ) b
  ) as testing group by "ParAdd", "ParZip" ORDER BY "ParAdd", "ParZip"


texts: Disable QGis Tips Panel when application start

texts:
I am trying to disable the QGIS Tips Panel with config files. I use QGIS Lyon.
- QGIS2.ini
- And a custom file define with the parameter --customizationfile

I know the user can click on "I've had enough tips, don't show this on start up any more!" and he never see this window again. But I would like this Panel never show (even the first time).

I start QGIS with this command

qgis.bat --configpath 'C:\gis' --optionspath 'C:\gis' --customizationfile 'C:\gis\QGISCONFIG.ini' --lang en_US

Maybe with a plugin ?


texts: What is the most convenient coordinates format to store data for spatial queries: WGS84 or OSGB36?

texts:
I have datasets that come in different coordinate formats as WGS84 (4326) and OSGB36 (27700).
Clients are consuming the "published" data as format WGS84.
So at the moment, all the data since the earliest stages are saved as WGS84.
However all the intermediate transformations are based on ST_Intersects, ST_DWithin, etc. So can I consider as rule of thumb: working with data as OSGB36 for faster calculations (cartesian coordinates) and converting to WGS84 only for the published data?


texts:
I'm currently developing an indoor positioning system for a college project. I searched for a good mapping platform and I found that openstreetmaps would suit me best, but now I'm stuck in generating an accurate indoor map with real long and lat values. How could I create an accurate indoor map with long and lat values without using GPS?

Any idea how to create such a map or any suggestions about accurate indoor mapping?


texts: Find out what Expression has been used to populate an attribute column in QGIS

texts: I need to find out what expression has been used in the field calculator to populate a number of columns in the the attribute table in QGIS. Does anyone know how to get it to display this info? Like in excel when you select a cell with a formula it will display the formula in the box at the top of the page.


texts:
How can I detect polygons / loops which are not topologically closed?  I'm not wanting to fix any open polygon / loop.  I merely just wish to identify them.

I do know that for a polygon to be "valid" it must be closed when using st_IsValid.  But an invalid geometry doesn't mean it's an open polygon.

I've been able to use st_IsValidDetail to help me find things like self-intersections, but now I just want to simply identify features which are not topologically closed.

I believe st_IsValidDetail should be able to do this for me quite easily if I know the correct reason to filter for.

The problem I run into is I do not know where online I can find a list of all the reasons for an invalid geometry as reported by st_IsValidDetail.  Is there a comprehensive list somewhere?


texts: Editing many-to-many relationship classes with ArcGIS API for JavaScript?

texts:
I am seeking some guidance on working with many-to-many relationship classes thru the ArcGIS Server Javascript API. I am designing a traffic count data management application. Count data locations are stored in a points layer, and there is a non-spatial table of count requests. My plan was to relate the two with a many-to-many relationship class. Count points can be requested any number of times, and each count request can request new data any any number of pre-existing count locations. A many-to-many relationship class seems like the obvious choice here because it enforces data integrity.

There is a complicated process to relate records thru a many-to-many relationship class in ArcGIS Desktop. I had first tried doing this thru adding rows to the relationship class table, but later found out this is not supported. This leads to two questions as I look to implement this functionality in my web application:

Will I be able to access a many-to-many relationship class table thru Javascript and edit/delete records using .applyEdits()?
Is a many-to-many relationship class really necessary in this case? Or is there a more straightforward way to accomplish the same result?


texts:
Can parametric SQL View be used for postgis raster map in geoserver?
Or how can I query postgis raster map with the dynamic request?


texts: create a fishnet using lat/lon table text file

texts:
I have a text file with lat/lon table values and want to use this information to generate a fishnet using python. Any heads up would be deeply appreciative. Also please point to if this has been already addressed.

Thanks,
Ibe


texts: Rotating annotation by angle using center point of annotations (with baseline alignment) in ArcMap

texts:
Is it possible to rotate annotation (in ESRI Annotation FeatureClass)  by angle so that new position would be calculated based on annotation center point (not alignment, as it is set to baseline, center)?

OR alternatively- would it be possible to change annotation alignment to center/center without moving annotation?

I have thousands of annotations and I would like to update their rotation based on updated angle field in database. All annotations have VERTICALALIGNMENT set to “baseline”.


texts:
How would I import a .csv file with polygon coordinates like the sample below into QGIS?

I tried using the tool to convert geotrace to WKT and then importing the WKT via importing a delimited file but the information came in as a dbf/table.

  {"geodesic":false,"type":"Polygon","coordinates":[[[-77.89444542198672,40.73437605298494],[-77.89417592740148,40.73437605298494],[-77.89417592740148,40.734915042155414],[-77.89444542198672,40.734915042155414],[-77.89444542198672,40.73437605298494]]]}


texts: Add layer to map in qgis3 via a python script

texts:
I'm trying convert my old plugin to QGIS3 and one part is to create a shapefile, and afterwards load it to the Layer menu like this:

from qgis.core import QgsVectorLayer, QgsProject
layer = QgsVectorLayer('c:/dev/test.shp', 'test', "ogr")
print(layer.isValid())
QgsProject.instance().addMapLayer(layer)

The python console returns True and the Layer appears in the layer list, hence I can't add features to it when I have toggled editing:

However, if I remove the layer from layer list and add the layer with the same code in the Python console I have no problem with adding a feature.

Any suggestions for what I can try?


texts: No way to De-authorise ArcGIS Server and regain license?

texts:
At v10 of ArcGIS Server, from what I understand of the new license model, it is not possible to de-authorise ArcGIS Server (like you can with Desktop)

Clients have visibility over their licences with the handy Customer Care Portal tool.
- From here you can see how many licenses you have authorised and how many you still have available.

Scenario:
- Purchase 1 license for ArcGIS Server (license level irrelevant) CCP: 0 Auth, 1 Avail)
- Install & Authorise in a dev environment   CCP: 1 Auth, 0 Avail)
- Roll over to production, install and attempt to authorise = Fails as no license available
= No ability to release license, so have to get in touch with ESRI Inc Customer Care to reset the license

This might not be too tricky if your over in States, but what about the rest of the world?  Get in touch with distributor --> to get in touch with Redlands/Customer Care --> and then back again = Takes lot of time.

Just wondering if I am either missing something or if there is a better workaround to the above?

At 9.3.1 and prior, it was possible to install using the same ECP several times, and through the online registration process, if Redlands saw that an ECP with one server license was getting authorised lots of times, then the client would get a phone call.
= I imagine the new model tries to take the work away from ESRIs end, but can be inefficient when trying to roll over.

I also assume that if my machine died, and I had to reinstall again, still would have to get in touch with Customer Care to reset license?


texts: Combining overviews for a large multiband raster image?

texts: I have several large rasters, complete with overviews. I have combined these as separate bands of a multiband virtual raster. However, QGIS performs very slow viewing this vrt. Is there a fast way to produce an overview for this vrt?


texts: Extracting building footprints from Sentinel-2 imagery

texts:
My attempt to extract building footprints from Sentinel-2 images using machine learning algorithm trained on Sentinel-2 images produced a lot of false positives and there is no sign that the algorithm actually learnt anything. When I tried the same architecture on another kind of dataset (MNIST, CIFAR-10), it worked perfectly.

My question is: Does it really make sense to use Sentinel-2 imagery with resolutions 10m for building extraction application?

Red = False Positives; Blue = Labels; Green = True Positives


texts: I want to have multiple maps in layout view, based on selected spatial data. Is this possible?

texts: I have a point data set with meta data. Based on this meta data I want to display the selected data in two or more different layout view windows in ArcMap 10.3 at the same time. Is this possible?


texts: Increase resolution of generated hill shade on QGIS

texts:
I want to export a hill shade (GeoTIFF) created from a DEM as an image for printing, but the resolution of the hill shade is way too low (900 columns, 1450 lines).

Any way to upscale it?

I tried raster > projection > warp, checked 'change size' and put higher numbers but it didn't make any difference.


texts: Clipping Basemap using shapefile as clip feature in ArcGIS Pro

texts:
I am trying to isolate the United States from the rest of the world using a provided basemap in ArcGIS Pro.
I am trying to use a shapefile of the United States I obtained from the living atlas as the Clip Feature. When I open the clip geoprocessing tool the basemap does not appear in any of the drop downs. When I drag the basemap to the input I am given the error 000840 which states, "The value is not a Feature Layer."
Is there anyway around this or are basemaps not editable?


texts: Changing DEM resolution when creating a slope raster?

texts:
I want to create a slope map using a LiDAR DEM. I think that I should use the slope tool; but I need a slope map with a 2m resolution and I'm not sure how to do that.

How can I have a slope map with a 2m resolution using a DEM with 1m resolution. See DEM's characteristics in picture below.


texts: Count the number of segments of a polygon based on the angles

texts:
I wish to count the number of segments of a polygon that has a length less than 25 meters when the angle between two vertices is not close to 180 degrees i.e. if the angle is between 178 to 182 degrees I want to ignore this case even though the segment length is less than 25 meters. I note the basic part of the code below is in VBA. Could anybody help me please to find why I get the wrong results? Thanks.

     Dim pFCursor As IFeatureCursor
Set pFCursor = pFeatureClass.Update(Nothing, True)

Dim pFeature As IFeature
Set pFeature = pFCursor.NextFeature

Dim pGeomColl As IGeometryCollection
Dim pSegColl As ISegmentCollection
Dim pCurve As ICurve
Dim pLine1 As ILine
Dim pLine2 As ILine
Dim dVtxAng As Double
Dim i As Long, j As Long
Dim pi As Double
pi = Atn(1) * 4

Dim counter As Integer
counter = 0
Dim angleDeg As Double

Do Until pFeature Is Nothing
FID = pFeature.Value(intPosFID)

pFilter.WhereClause = "FID=" & FID
pActiveView.PartialRefresh esriViewGeoSelection, Nothing, Nothing

Set pSelection = pLayer
pSelection.SelectFeatures pFilter, esriSelectionResultNew, True
'******************Code for measuring angles and minimum lenghts

    ' Get all disjoint paths
    Set pGeomColl = pFeature.Shape

    ' Loop thru the paths
    For i = 0 To pGeomColl.GeometryCount - 1
        ' Get all segments that make up this path
        Set pSegColl = pGeomColl.Geometry(i)
        Set pCurve = pSegColl

        ' Loop thru the segments
        For j = 0 To pSegColl.SegmentCount - 2

            ' Get the next two lines
            Set pLine1 = pSegColl.Segment(j)
            Set pLine2 = pSegColl.Segment(j + 1)

            ' Calculate the left side angle (in degrees)
            dVtxAng = CalcAngleBetweenLines(pLine1, pLine2)
             angleDeg = dVtxAng * 57.2957795

            If angleDeg >= 178 And angleDeg <= 182 Then
            counter = counter
            'if angle is close to 180 degrees then ignore the points and its vertices
    GoTo Nextsegment
            End If

            'if angle is not close to 180 degrees and lenght of segment is less than 25m then count the segment
    If angleDeg < 178 Or angleDeg > 182 And pLine1.Length < 25 Then
             counter = counter + 1

            End If

Nextsegment: Next j

Next i

'Store number of points in the field
pFeature.Value(intPosminLenght) = counter
pFCursor.UpdateFeature pFeature

counter = 0
Set pFeature = pFCursor.NextFeature
Loop

End Sub


texts: Avoid backtracking while creating loop in Google Map API 3

texts:
I am trying to build a solution where i would create loop by google map. My first and only business rule is there should not be any repeated route. But with my demo http://www.mindcrew.in/road%20runner/example.html it does show.

Please help me on this .


texts: Is it possible to use UTF-8 characters for shapefile attributes?

texts:
I heard that according to a specification shapefile attributes stored in dbf-files cannot have UTF-8 as encoding.
And then I found out that an additional file with the ending CPG can be used to define the code page to identify the character set used in the dbf file.

For me this sounds contradictory.

So my question is now - is it possible to use UTF-8 characters for shapefile attributes?

and - additional question, since encoding and character sets are always an unpopular topic, do You know a page where the topic is described in detail?


texts: Selecting attributes by coordinates from DataFrame in Python

texts:
I wonder how to return the right information needed from this df. Currently I get an empty return. I think the problem is my selection because the coordinates are correct and all variables are as well. If I take away just one & from the selection it seems to work,  The code is running like this, you can try your best. I really think it's just a little thing missing here.
Who can help?
import os, requests
import pandas as pd
pd.set_option('display.float_format', lambda x: '%.11f' % x)

#example bounds
xmax, ymax, xmin, ymin = 3845447.52241509, 8560465.238354621, -2057434.3252381992, 4031530.0276551154

direc = os.getcwd().replace('\\', '/')
#wd
url = 'https://modis-land.gsfc.nasa.gov/pdf/sn_bound_10deg.txt'
tiles_txt = requests.get(url, allow_redirects=True)

open(direc+'/tiles_bounds.txt', 'wb').write(tiles_txt.content)
text = open(direc+'/tiles_bounds.txt', 'r').read()
text = text.splitlines()
del text[0:7]
del text[-2:]

v,h,hv = [],[],[]
for iv in range(len(text)):
    v.append('v' + text[iv][1:3].replace(' ', '0'))
for ih in range(len(text)):
    h.append('h' + text[ih][5:7].replace(' ', '0'))
for VH in range(len(text)):
    hv.append( h[VH] + v[VH] )

df_bounds = pd.DataFrame(columns=['tile','horizontal', 'vertical','xmin_m','xmax_m', 'ymin_m', 'ymax_m'])
df_bounds['tile'] = hv

hor, ver, xmin_m, xmax_m, ymin_m, ymax_m=[], [], [], [],[], []
for h in range(len(df_bounds['tile'])):
    hor.append(df_bounds['tile'][h][1:3])
    ver.append(df_bounds['tile'][h][4:6])
df_bounds['horizontal'] = hor
df_bounds['vertical'] = ver
df_bounds['horizontal'] = df_bounds['horizontal'].astype(int)
df_bounds['vertical'] = df_bounds['vertical'].astype(int)

tile_width = 1111950.5196666666
tile_height = 1111950.5196666666

for i in range(len(df_bounds['tile'])):
    xmin_m.append(-20015109.354 + df_bounds['horizontal'][i] * tile_width)
    ymin_m.append(-10007554.677 + (17 - df_bounds['vertical'][i]) * tile_height)
df_bounds['xmin_m'] = xmin_m
df_bounds['ymin_m'] = ymin_m

for i in range(len(df_bounds['tile'])):
    xmax_m.append(df_bounds['xmin_m'][i] + tile_width)
    ymax_m.append(df_bounds['ymin_m'][i] + tile_height)

df_bounds['xmax_m'] = xmax_m
df_bounds['ymax_m'] = ymax_m
#####here its's broken HERE
################################################################################################
tile_names = df_bounds['tile'][(df_bounds['xmin_m'] >= xmin) & (df_bounds['xmax_m'] <= xmax) & (df_bounds['ymax_m'] <= ymax) & df_bounds['ymin_m'] >= ymin]
####################################################################################################
tile_names = [i + ',' for i  in tile_names]
tile_names = list(tile_names)
tile_names = ''.join(tile_names)
tile_names = tile_names[:-1]
os.remove(direc + '/tiles_bounds.txt')

print(df_bounds, '\n', '\n')
print(tile_names)


texts:
I"m trying to nest expressions in this raster calculator conditional statement, but I get an error.

This statement works, so I tried to follow it as a template:

Con (IsNull("ras1"),"ras2","ras3")

Separately, this statement works

"rasB"-(5+(1/(1+"rasC"))*5)*("RasD"<=20)

And I'd love to merge them into one calulation liek this

Con(IsNull("rasA"),
("rasB"-(5+(1/(1+"rasC"))*5)*("RasD"<=20),
("rasB"-(5+(1/(1+"rasC"))*5)*("RasD"<=500))))

Where "rasA" is null, return "rasB"-(5+(1/(1+"rasC"))5)("RasD"<=20). Otherwise return "rasB"-(5+(1/(1+"rasC"))5)("RasD"<=500)

I get the error shown below so I'm trying to figure out if Con() requires a layer like the error says or if can I use an expression and I just entered it wrong.


texts: Tool (Simplify Polygon) parameter is not visible anymore

texts:
I am having a strange problem in ArcGIS Toolbox.

I am trying to use the Simplify Polygon (Cartography) tool from ArcToolbox. When I open the tool, there is no possibility to input Simplification Tolerance. Please see the attached image.

I can run the tool as a batch and from the Python window without any problem.

I am using ArcGIS 10(Editor) with Service Pack 5 on Windows 7 64Bit (English) and Windows Server 2008 64Bit (German). On both machines, I have the same problem.
Any ideas?

I have reinstalled the software but the problem remains.


texts: How to map districts without showing features in other districts?

texts:
We want to have a map of each individual district, with the shops on it, so that the person responsible for each district can have a map of their territory. The reason I want to merge them is so that I don't end up with a map of one district with lots of dots around it from other shops all over the state.

I'm using QGIS 1.7.0.  I have a shapefile that splits my state into districts, and a geoJSON that has a list of shop locations.  I want to make a file of each district, with a dot at each shop.  Is it possible to have QGIS put both pieces of information into one file (I haven't had any luck with merge shapefiles).  After that, it would be wonderful if I could get it to split each section into it's own file, with the dots that it contains.  If there is not a way to get it to do that automatically, would the most efficient thing to do be putting a mask layer over the map and just cutting out a hole for each district, one at a time?

I'm mostly self-taught with QGIS, so I don't necessarily know some of the non-basic stuff.  But I'd appreciate any help y'all could give me.


texts: Select points situated on zero degree slope

texts:
I need to choose points from a shapefile, which are positioned on raster cells with 0 degree slope. Creating a slope raster from DME is not a problem, but I don't know how to write a script that chooses points positioned only on zero slope.

I don't know if use GDAL or it's possible over ArcPy or anything.


texts:
I would like to create a grid (say 10mx10m) with unique values (i.e. grid 1,2,3,4) and then overlay a point file and assign a grid value to each point. All points within the first grid have field grid ID = 1 all points within the second grid have field grid ID = 2, etc.

I have created grids with Xtools Pro> Create Fishnet which have unique FID but haven't figured out how to get that associated with the point data.


texts: Split polygon manually- geometry validity error problem

texts:
I want to split a polygon layer manually using the Split tool in QGIS however I'm getting an error

  no feature slit done: The geometry is invalid. Please repair before
  trying to split.

I have checked the geometry using Vector-Geometry Tools-Check Geometry Validity and it reports 0 errors encountered.

I have tried creating a new simple multi-polygon layer and split works as suggested.

I can split the other polygons in this same layer, just not the one I'm interested in.

Any ideas for what else I can try?


texts: Dynamic legend with multiple layers in leaflet

texts:
I have a choropleth map I am putting together with two different layers. each layer has its own symbolization and legend. I have so far been successful in having the legend change when a new layer is turned on, however if a layer is turned off, the legend does not change to match the remaining layer. Additionally, if all layers are turned off, the legend will remain for whatever layer last activated a legend change.

my code that is working for the layer change so far is:

map.on('overlayadd', function (eventLayer)
{
    if (eventLayer.name === 'Housing'){
        map.removeControl(currentLegend );
        currentLegend = houselegend;
        houselegend.addTo(map);
    }
    else if (eventLayer.name === 'Families') {
        map.removeControl(currentLegend );
        currentLegend = childlegend;
        childlegend.addTo(map);
    }
})

The code I have tried so far to make the legend change when a layer is removed or all layers are removed is:

function legendChange(e) {
    if (map.hasLayer(Families) && !map.hasLayer(Housing)){
        map.removeControl(currentLegend );
        currentLegend = childlegend;
        childlegend.addTo(map);
    }
    else if (map.hasLayer(Housing) && !map.hasLayer(Families)){
        map.removeControl(currentLegend );
        currentLegend = houselegend;
        houselegend.addTo(map);
    }
    else if (!map.hasLayer(Housing) && !map.hasLayer(Families)){
        map.removeControl(currentLegend );
    }
}

The map loads right now, so I'm confident my syntax is close and that the problem is likely with my logic.


texts: Deploying services, data, & code to an ArcGIS production server?

texts:
We have a production GIS viewer that runs on a production server. Testing new development/data is done on a test server and we develop on a different server also. At the moment, a lot of the deployment is manual work, with risk for errors. I would like to automate whatever is possible to reduce that risk.

I was want to maintain and deploy new data (Almost all data for our viewer is kept in File geodatabase format).

My main questions are these:

What kind of scripts/tools do you use to deploy new maps,services,data to your production server,

do you deploy them manually? In ArcCatalog or the web manager?
do you copy and paste the cfg files, do you change the settings of those files manually or with scripts?

Do you keep stuff in version control? What do you keep there?
How do you manage different settings between servers (min and max number of SOC's per service for example)
What else do you recommend? Are there good docs to read about it?


texts: Why is the pdf output of my QGIS map different from the one I see in the Print Composer?

texts:
I'm working on a map of migrant deaths and missings in and around the Medditerranean Sea. It includes a legend. When I work in the print composer mode, the legend seems to be looking fine. Here it is:

When I look at the pdf output of my map, however, it looks different:

The key issues here are the extra pink moon-like shape on the left, and the extra "1" at the bottom of the graduated circles.
I suppose this has something to do with the fact that I'm trying to hide part of the "wrong" graduated circles, which I'm overlapping with two images: correct version and a white image to cover the remaining bits. One can see the outlines of these images over here:

Question: why is the output of the legend in the pdf of the map different from the one in the print layout, and how can I make sure QGIS generates the correct output?


texts: Average value of polygons within a set distance from a point

texts:
I have two vector layers: one polygon and one point layer. The polygon layer has a value from 1 to 10. I am trying to find the average value for all polygons within or partly within a 1-mile radius of each of the points in the points layer. What would be the best way to go about this?

I am using ArcGIS 10.2 and I have access to the Spatial Analyst extension.


texts:
I try for too ling to read the envi file in python. Normally this is very easy.
First, I extract necessary values from the .hdr:

for line in open("../../local_storage_only/aviris/ang20180626t145910_rfl_v2q2/ang20180626t145910_corr_v2q2_img.hdr"):

  if "data type" in line:
    _data_type = int(line.split("= ")[-1])
  if "interleave" in line:
    _interleave = line.split("= ")[-1]
  if "byte order " in line:
    _byte_order = int(line.split("= ")[-1])
  if "wavelength units" in line:
    _wavelength_units = line.split("= ")[-1]
  if "correction factors" in line:
    _correction_factors = line.split(" { ")[-1][:-2]
    _correction_factors = [float(element) for element in _correction_factors.split(" , ")]
  if "wavelength = " in line:
    _wavelength = line.split(" { ")[-1][:-2]
    _wavelength = [float(element) for element in _wavelength.split(" , ")]
    _wavelength_dict = dict(zip( _wavelength, range(len( _wavelength))))
  if "bands = " in line:
    _bands = int(line.split(" = ")[-1])
  if "lines =" in line:
    _rows = int(line.split(" = ")[-1])
  if "samples =" in line:
    _columns = int(line.split(" = ")[-1])
  if "data ignore value " in line:
    _ignore_value = int(line.split(" = ")[-1])
print(line)

Secondly, I open the file and convert it in to the right format:

tmp_file = open("../../local_storage_only/aviris/ang20180626t145910_rfl_v2q2/ang20180626t145910_corr_v2q2_img",
               encoding="utf-8")
tmp_img = np.fromfile(tmp_file,
                      np.dtype(np.float32).char)

tmp_img[tmp_img == _ignore_value] = np.nan

However, if I plot the results it looks like this:

I also tried to change the order of height, width and bands. I altered the sorting mechanism of np.reshape("A", "C", "F"). Does someone know what I do wrong?


texts: Add new spatial reference to ArcGIS with WKT

texts:
I'm working with a layer on ArcMap with a spatial reference that is not included in ArcMap. How can I create and add the spatial reference because I already know its WKT?
I use ArcMap 10.3


texts: Run Shortest Path Dijkstra with pgrouting

texts:
I use PostgreSQL 9.5.3, PostGIS 2.2 and pgRouting 2.3.2.

I have  the same problem as Shortest Path Dijkstra because I use the same data with the same columns but when I put

 SELECT seq, id1 AS node, id2 AS edge, cost, b.the_geom, b.osm_name FROM pgr_dijkstra('

            SELECT id,
                     source::integer,
                     target::integer,
                     length::double precision AS cost
                    FROM at_2po_4pgr',
            30, 60, false, false) a LEFT JOIN at_2po_4pgr b ON (a.id2 = b.id);

into SQL window I got

  ERROR: column link "cost" is ambiguous

How do I change cost atrribute?

I made it using this code:

SELECT seq, id1 AS node, id2 AS edge, di.cost, geom_way, osm_name  FROM pgr_dijkstra(
'SELECT id, source, target, st_length(geom_way) as cost FROM aa_2po_4pgr',
23941, 24159, false, false  ) as di  JOIN aa_2po_4pgr  ON di.id2 = aa_2po_4pgr.id


texts: Separating a shapefile into separate shapefile based upon their features using a python script

texts: I have a shapefile of a US. I need to make an individual shapefile for all the states in the US using a python in ArcGIS. How can I do it?


texts: Convert CSV file data points to line feature in QGIS

texts:
I have added csv file data points into QGIS. Then I opened properties and tried to convert points inside the plot as single line. (used option ---> properties ---> symbology --->  marker--->geometry generator --->Line string) But, I can't see any data in my plot after selecting this option.
Is there any other way to convert points to line?


texts: Getting Country from GeoFeatureLayer using ArcObjects?

texts:
I'm currently using ArcGIS and I'm trying to create a basic application where I can click on a country and then get back information such as the shape, size, etc.

I've looked at the code snippets and the documentation and up to date this is what I've come up with.

IPoint pt = new PointClass();
pt.X = e.X;
pt.Y = e.Y;

IGeoFeatureLayer layer = myMapControl.get_Layer(0) as IGeoFeatureLayer;
ESRI.ArcGIS.Geodatabase.IFeatureCursor cursor = GetAllFeaturesFromPointSearchInGeoFeatureLayer(1.0 , IPoint , layer , myMapControl.ActiveView)

So at this point I can click on the map, get the X,Y coords from the mouse  - get the base layer from the map control and then get the FeatureCursor.  Looking at the cursor info I do see things like CNTRY_NAME but nothing else of use.  Am I going about this in the totally wrong direction?  I've looked at http://help.arcgis.com/en/sdk/10.0/arcobjects_net/conceptualhelp/index.html quite a bit lately in addition to this site, but I'm coming up blank.


texts: Display a fixed title (img or txt) on a Mapbox map

texts: I don't find anything on the Mapbox documentation about this ; so I guess I've missed the logic : is there a simple way to add a title on the map that is not affected by the zoom ? (like the Mapbox logo at the bottom left)


texts: How to identify correct key_on parameter from JSON file for Folium Choropeth

texts:
I'm trying to make a Folium Choropleth map. But I'm stuck for long time on completing it. My map shows boundaries correctly, even shows legend correctly but the map is not filled with choropleth colors.
I'm running my code on Jupyter Notebook. I have combined the geopandas dataframe and pandas dataframe as df_mergeplz
My simple code is
plzchoromap = folium.Map(location = [50.935173,6.953101],tiles='cartodbpositron',zoom_start = 11)

plzchoromap.choropleth(geo_data=df_mergeplz,
                       data = df_mergeplz,
                       columns=["Postal Code","Population"],
                       key_on="feature.NUMMER",
                       fill_color='YlGnBu',
                       fill_opacity = 0.3,
                       line_weight=2,)
plzchoromap

my JSON file looks like
{
 "displayFieldName": "NUMMER",
 "fieldAliases": {"OBJECTID": "OBJECTID","NUMMER": "Postzustellbezirk (Nr.)"},
 "geometryType": "esriGeometryPolygon",
 "spatialReference": {"wkid": 4326,"latestWkid": 4326},
 "fields": [
      {"name": "OBJECTID","type": "esriFieldTypeOID","alias": "OBJECTID"},
      {"name": "NUMMER","type": "esriFieldTypeString","alias": "Postzustellbezirk (Nr.)","length": 5}
            ],
 "features": [
      {"attributes": {"OBJECTID": 16,"NUMMER": "51065"},
      "geometry": {"rings": [[[7.0088577711265421,50.964379468107985],,,....,]]}
      },
      {"attributes": {"OBJECTID": 18,"NUMMER": "50823"},
      "geometry": {"rings": [[[6.9339208198918785,50.96181483933605],,,....,]]}
      },
 .
 .
 .
 .
 .

   }

I want to use the parameter NUMMER inside features.attributes.NUMMER.
But passing features.attributes.NUMMER throws error that
~/opt/anaconda3/envs/ML1env/lib/python3.8/site-packages/folium/folium.py in get_by_key(obj, key)
    285             def get_by_key(obj, key):
    286                 return (obj.get(key, None) if len(key.split('.')) <= 1 else
--> 287                         get_by_key(obj.get(key.split('.')[0], None),
    288                                    '.'.join(key.split('.')[1:])))
    289

~/opt/anaconda3/envs/ML1env/lib/python3.8/site-packages/folium/folium.py in get_by_key(obj, key)
    284
    285             def get_by_key(obj, key):
--> 286                 return (obj.get(key, None) if len(key.split('.')) <= 1 else
    287                         get_by_key(obj.get(key.split('.')[0], None),
    288                                    '.'.join(key.split('.')[1:])))

AttributeError: 'NoneType' object has no attribute 'get'

I have referred several other posts which primarily showed correction on key_on parameter. I doubt I am missing something obvious here.
Summary:
key_on = feature.NUMMER #outputs unfilled folium map
key_on = feature.attributes.NUMMER #throws error.

Can someone please let me know what is the right key_on parameter for me or a guide on how to play with JSON object to find the correct key_on parameter.
Just to add on my df_mergeplz looks like below:
OBJECTID    NUMMER  geometry    minx    miny    maxx    maxy    Unnamed: 0  Postal Code City    Administrative Region   Population  Area    Area Codes  Neighborhoods   Latitude    Longitude
0   1   50769   POLYGON ((6.84866 51.08444, 6.84873 51.08420, ...   6.772530    51.007784   6.973939    51.084964   15  50769   Cologne North Rhine-Westphalia  63970   44.800  2203, 221   Chorweiler  51.046949   6.877612


texts: Reverse parameters of a projection map / guess which geographic project is used

texts:
I'm user of a soaring simulator named "Condor Soaring".

This simulator is at the center of a community of enthusiastic people who build sceneries
(landscapes) which represents quite accurately real world.

Glider pilots can train to achieve some tasks (such as flying a triangular circuit or other shapes).

Track distance can be quite long (more than 200 km). Time to achieve task is also important.
Tactical glide computer can be used (like XCSoar) (both for real gliding
but also for simulated gliding)

Condor Soaring tasks are stored in a text file (like a INI file).
XCSoar circuits are stored in a text file (XML file).

I wrote Python script to convert Condor task file to XCSoar
https://github.com/scls19fr/pycondor

XCSoar is using geographical coordinates but Condor soaring is using X, Y coordinates.
Maximal values for X and Y depend of landscape considered.

A DLL named NaviCon.dll provides some functions to get theses informations about a given landscape
(GetMaxX, GetMaxY)

This DLL also provides functions XYToLat, XYToLon to convert respectively XY coordinates
to latitude and to longitude.

This part of my work development is working fine but I would like to be able to convert
XY coordinates to Lat, Lon coordinates without the use of this DLL.

I wrote a script which output using DLL calls, for each installed landscape, max_x and max_y
and have coordinates (lat, lon) of the 4 corners.

    PosX    PosY    Lat         Lon
0   0       0       43,22330093 7,807540417
1   207360  0       43,1682663  5,257740498
2   0       345600  46,33418655 7,741717339
3   207360  345600  46,272892   5,0517416

I also did a script which output this kind of plot (for 20x20 regularly xy spaced points)

a scatter plot with 40x40 points

CSV data can be found here

I wonder how I can know what kind of projection was done, what geodetic datum was "probably" used. I also would like to know what parameters for this projection should be set to be able to calculate myself Lat, Lon without DLL calls but with the use of a Python projection library such as pyproj or GDAL... )

PS : some other similar questions
Identifying Coordinate System of Shapefile when Unknown?
https://stackoverflow.com/questions/3203723/how-do-i-guess-which-geographic-projection-system-is-used-for-a-dataset


texts: Stop custom attribute form from resizing in Qgis

texts:
When I open a custom attribute form for a feature in Qgis it always gets resized to the default size of the original Qgis attribute form. I've tried to set the sizePolicy in QT Designer to Fixed instead of Preferred to prohibit this, but it didn't help.

If I resize the form when I open it by clicking on a feature in Qgis and then resize it, it still opens with Qgis default height/width when I click the next feature.

Is there any kind of setting to tell Qgis to use the size that the form is designed to?

(Running Qgis 1.8.0, QT Designer 4.8.4 on Windows 7)


texts: Compute difference in lon / lat when moving at angle

texts:
I have the following problem: I have been given an origin as a lon / lat pair, a distance in kilometers and an angle on the ground, relative to the cardinal directions. My goal is to compute the destination (as a lon / lat pair) which is reached after moving in the direction given by the angle for the given distance.

How might I go about computing the destination? Is there a closed formula?


texts:
If you were to make viewshed analysis for a point in the landscape and whether you can see a smoke in the distance, what would be the total area that you would calculate? 20 km from the observer point, 30 km? Or less?

Elevation of observer point is 400m, the source fire is on a mountain, below it's peak, on 950m. Distance is 25km, between them there is flat ground and clear atmosphere. Obviously, the viewshed analysis shows that it's visible. But how clearly visible? That's what I'm wondering.


texts: QGIS Dissolve tool re-projects and rotates layer

texts:
I'm trying to use QGIS (2.6.1-Brighton) to aggregate polygons according to an ID attribute.

I'm using the shapefile of school district polygons in Michigan found here:

http://www.mcgi.state.mi.us/mgdl/?rel=ext&action=sext

MI Geographic Framework School Districts (v14a) is the name of the shapefile I'm working with.

I'm trying to create a shapefile of Intermediate School Districts, which in Michigan are a level of aggregation above school districts.

I should be able to do this using the dissolve tool, as there is an ID variable in the v14a shapefile (called ISD).

However, upon doing so, the resulting polygon layer appears to have been re-projected and rotated.

Here's the input polygons:

The CRS of this layer is strange, but here's the definition:

+proj=omerc +lat_0=45.30916666666666 +lonc=-86 +alpha=337.25556
+k=0.9996 +x_0=2546731.496 +y_0=-4354009.816 +no_uoff +ellps=GRS80
+towgs84=0,0,0,0,0,0,0 +units=m +no_defs

And here's the result of the dissolve:

With CRS:

+proj=omerc +lat_0=45.30916666666666 +lonc=-86 +alpha=337.25556
+k=0.9996 +x_0=2546731.496 +y_0=-4354009.816 +gamma=0 +ellps=GRS80
+towgs84=0,0,0,0,0,0,0 +units=m +no_defs

Any idea what's going on here??

I would be particularly grateful if there was a way to reproject the rotated layer so as to undo the rotation, as I had already done a fair amount of stuff by hand to add new attributes to the merged layer (so that starting from scratch would require a lot of repeated manual labor).


texts: Understanding why ArcSDE FeatureLayer is read only and FeatureTypes property is null?

texts:
I have a WPF FeatureLayer backed by an ArcSDE 10 feature class that I hope to use for editing. After the layer is done initializing its readOnly property is true, I think that this points to some problem with the ArcSDE setup but I can't figure out what it could be.

I've verified that the connection is working, inspecting the FeatureLayer shows that it has downloaded the remote table name. The post installation setup for ArcSDE had no trouble, same with creating some empty feature classes, adding them to an mxd file and using it to host a feature service. Is there any other setup I may have missed?


texts: Tabulate Area tool "doubling up" for Overlapping Buffer Zones in ArcGIS Desktop?

texts:
This question is related to a previous question I asked as Raster Version of Spatial Join (one to many)

In ArcGIS Desktop 10.1 I used the Tabulate Area tool to crosstabulate values between an overlapping "buffrast" raster (with identities of 1917 different buffer zones) and "RR1.tif" (with land cover codes for each pixel within the 1917 buffer zones). Detailed methods:

Region Group was carried out for RR1.tif with the Add Link option chosen, and the Snap Raster environment enabled for buffrast. Tabulate Area (Spatial Analyst) was used with the following inputs: input raster or feature zone data = buffrast (zone field = ident), input raster or feature class data = RR1.tif (class field = LINK), processing size = 25.

This worked perfectly for most of the points, with total areas adding up to ~70000 for each of the 1917 buffer zones; however I run into problems where two or more buffer zones overlap each other - when this happens ArcMap "doubles up" on the total area of each buffer zone resulting in Areas of ~90000 and up.

How can I run the Tabulate Area tool for such overlapping buffer zones?


texts: Adjusting multiple text elements in multiple mxd files?

texts:
I am trying to make changes to titles on multiple mxds instead of having to open 1 by 1 but it is only making changes to 1 mxd file, I have even tried to make duplicate copies of the mxd to test but again it only makes 1 change. The code runs with no issues, am I missing something ?

import arcpy, os
from os import listdir

mxdPath = r"C:\Users\Admin\Desktop\Final Maps\New folder"

mxdFiles = [ f for f in listdir(mxdPath) if os.path.isfile(os.path.join(mxdPath,f)) and f[-3:] == 'mxd' ]

try:

    map = arcpy.mapping
    for mxdFile in mxdFiles:
        mxd = map.MapDocument(mxdPath + "\\" + mxdFile)
        for elm in map.ListLayoutElements(mxd,"TEXT_ELEMENT"):
            if elm.text == "CountyProject130":
                elm.text = "Testing709"

    mxd.save()
    del mxd
    Print "Complete!"

except Exception as e:
    print(e.message)


texts: SSRS compatible shapefile for post codes?

texts:
I am trying to make some data visual for a report to show how we are breaking up postcode regions.

I work at a university and we have a list of postcodes that we assign to campuses. I want to color the postcode regions on a map to a color representing the campus to make what area they cover visible.

I have found that SSRS has a map component, and that the Australian Bureau of Statistics releases post code shapefiles, but can't get the files to load in SSRS.

Would anyone know of a way to convert the shapefiles to get them working, or another way to get visualize this data?


texts: QGIS --- Create and move point using external python in an offline environment

texts:
I'm pretty new to QGIS and GIS.
I have an python application that get periodically position [lat,long,height] updates of several moving objects.

The idea is that the existing application sends the data to QGIs and QGIS prints the objects as a plain circle on the map. In the next loop, QGIS receives a updated set of the objects and moves/reprints them at the given location.
My Questions are:

How do I access the API from an external python script?
Do I have to create a layer for the objects?
How do I create a layer from an external script?
How do I send an object to QGIS or print it?
How do I move/update these obejcts or is it faster to delete and reprint all objects?
Could QGIS get the object information from a SQLite-DB?

Best Regards,
Sven


texts:
I'm pretty new to GIS matter and have no idea how to get the right projection of a point.

What's the projection for point (582865.200436, 3541848.067073)?

I have a .prj and corresponding .shp file. The following, is the content of '.prj' file: PROJCS["WGS_1984_UTM_Zone_39N",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",51.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]

I have found PROJECTION["Transverse_Mercator"] in the above content and based on real location of the file (approximately around Iran/Asia) i guess the corresponding projection should be 2059.

How can I find out the correct projection and unit of a '.shp' file out of it's corresponding '.prj' file or any other evidences?


texts: I need to extract a specific layer from a Geospatial PDF using GDAL/OGR

texts:
I have GDAL/OGR 1.10 installed on a 64bit Windows 7 machine. I have some Geospatial PDFs that I want to crop to the outline of the map and convert to Rasters.  However I cannot clip to the neatline.  They have been created using Esri Production Mapping, so the neatline is offset from the actual edge of the map, and clipping by the neatline includes the Axis labels and a lot of whitespace.

I can interrogate the Layers of the PDF by using gdalinfo -mdd LAYERS and I know which layer represents the actual neatline/the edge of the map.

What I need to know is how to use GDAL or OGR to convert this layer to a shapefile.


texts: shapely split function does not work on short lines

texts:
I am trying to split a line with multiple points. I adopt the solution in here, but sometimes it does not work.
from shapely.geometry import LineString, MultiPoint
from shapely.ops import split

line = LineString([(0.1, 0.01), (0.001, 0.001)])
splitter = MultiPoint([line.interpolate((i/6), normalized=True) for i in range(1, 6)])
gcline = split(line, splitter)
str(gcline)

The output is
'GEOMETRYCOLLECTION (LINESTRING (0.1 0.01, 0.001 0.001))'

I guess it is because the line is too short. The method works fine for LineString([(0.1, 0.1), (0.001, 0.001)]).


texts:
I have river networks (polyline) of several adjacent countries as shapefiles that almost match up at the borders, but not quite.

Is there any way to automatically connect them to a consecutive river network and define a tolerance distance. I was thinking along the lines of ArcGIS's Integrate feature but it didn't do the trick. Open to suggestions in ArcGIS 10 sp5 (or 10.1) and Qgis. Exact preservation of the geometry is not necessary.

Update: Sometimes there are only small gaps between the networks, but mostly they overlap in very weird ways, so edgematch and extend line are not an option. The following picture should clarify my problem:


texts: Exact values not possible in graduated colors symbology in Arcmap?

texts:
I have a shapefile with a field containing floating point numbers ranging from about -200 to 20. Many records have a value of exactly 0. I would like to symbolize the layer as follows:

< 0 (green)
= 0 (yellow)
`> 0 (red)

It looks as though in Symbology > Quantities > Graduated colors only ranges are possible. I first selected three classes and defined the bottom 2 cuttoffs as zero:

However when I try to set the color ramp in the symbology I see that an exact value of 0 is not possible:

Also checking the layer I see that all of my zero values have a red color, which isn't what I want.

How do I force the middle category to be exactly 0, and the top and bottom categories to exclude zero?


texts:
I'm trying to join up some data so I can make some rasters. I am doing this in Python but when I create the join in Python in results in an empty attribute table for my shapefile. The columns from both the shapefile and the view I am using in the join are in the shapefile after the join and my view shows up as a joined table in Joins and Relates but there are no data at all.

My Python is:

arcpy.AddJoin_management("Area1_wellspiezometers", "LocId_Calc", "_107062_1_1999", "LocationName")

Where

Area1_wellspiezometers is a shapefile
LocId_Calc is an field in Area1_wellspiezometers
_107062_1_1999 is a view created with arcpy.MakeTableView_management from a table in an MDB.
LocationName is a field in _107062_1_1999

Alternatively when I do this through the GUI I right click on Area1_wellspiezometers > Joins and Relates > Join... and enter my info as such:

This makes the join as expected and everything is dandy.

What is the difference between these two methods and how can I make the Python do what the GUI is doing?


texts:
I am trying to reclassify my vector layer but I seem to have a mistake in my Code Block. The following is an example of what i want;

Old  values:

10000
20000, 163000
230000 t/m 270000, 30000, 290000

New values:

10000
20000
30000

The code I have in my Code block is the following;

    def Reclass( Indexvalue ):
  if (Indexvalue == 10000):
    return 10000
  elif (Indexvalue == 20000) or (Indexvalue==163000 ):
    return 20000
  elif (Indexvalue>= 230000 and Indexvalue <= 270000) or (Indexvalue == 30000) or (Indexvalue == 290000):
    return 30000
  elif (Indexvalue>= 40000 and Indexvalue <= 44000 ):
    return 40000
  elif (Indexvalue>= 50000 and Indexvalue <= 73000 ):
    return 50000
  elif (Indexvalue>= 80000 and Indexvalue <= 97000 ):
    return 60000
  elif (Indexvalue>= 110000 and Indexvalue <= 115000) or (Indexvalue>= 130000 and Indexvalue <= 140000 ):
    return 70000
  elif (Indexvalue>= 150000 and Indexvalue <= 161000) or (Indexvalue ==220000) or (Indexvalue == 221000) or (Indexvalue>= 170000 and Indexvalue <= 194000 ):
    return 80000
  elif (Indexvalue>= 200000 and Indexvalue <= 202000) or (Indexvalue == 280000):
    return 90000
  elif (Indexvalue== 210000):
    return 100000

I found this code in a similar question like I had, but it does not work with me. Both VB and Python expressions give a syntax error so I think that the def and the if functions are not properly used by me.

Can someone point me in the right direction?


texts: Creating a contiguous raster stack from different sources (vectors and DTM bands)

texts:
This is a conceptual question I've been trying to implement for a while, and not being very good at raster/data processing or python I'd like to hear best practices in dealing with an extent such as the ConUSA.

Problem: I would like to create a tool which would allow me to append csv files with latitude and longitude with values from raster layers and polygons. Basically a spatial join tool but for standard layers...

Ideas: I've been struggling with whether to deal with shapes in postgis or flat files. The problem I have with postgis is I suck at the whole SQLAlchemy engine and such and never know how to properly apply it in a loop when pulling out chunks of shapes by county or state or things like that

My current thought process is this: Take the NLCD2014 in EPSG:5071, and create band with NED 1arc second after doing a gdal_translate of the that gives me 2 bands. Few concerns:

1 - Would it be best to keep the NED in it's 1 degree tiles, project them to 5071 and layer everything into 1 degree x 1 degree tiles. I know it would be easier to process than a merged tif... Are there any good resources or examples of this type of work?

2 - I believe it wouldn't be difficult to import a country wide layer of DFIRM polygons and export as 1 degree by 1 degree tiles... right? Is there a better way?

UPDATE WORKFLOW EXAMPLE:

Setup: I get a set of Latitude / Longitudes that I need to sift through and assign a height, floodplain, and our own proprietary flood depths from raster files we create. I'm trying to find the best way to loop through, I believe if I group into counties or HU10 or something and loop that way it would be best, the following is implemented in crappy python:

1) NED tiles - I created a vrt and then created a compressed .tif (lzw) and clipped to the ConUSA, this becomes one raster band I use as lookup for DTM on a country wide basis through gdal and geotransformations and pixels and stuff.

2) DFIRM - I use the polygon shapes as for point in polygon , been experimenting with trying to rasterize but I believe a PostGIS DB with a county or census block index would be quickest, I just struggle with PostGIS engine connection and such

3) Flood Raster - Depth of flood water on HU2 basis, right now we run 2D models, create return periods and I go through and stack the 10yr 50yr 100yr 250yr 500yr and 1000yr rasters into a 6 band stack by HU2 and pull out the depths

So after writing it all out it seems probably best to group everything into HU2 stacks right? Rasterize the DFIRM polys with integers corresponding to floodplains and do the same?


texts:
I am trying to export a map using Python/ArcPy script. The problem that the title of the exported map containts somthing like "A(c)" instead of the accent "é" I wrote. The date (dt) is printed (in french) with no problem.
This part of my scripts looks like:

    txt1="Densité  \n"+ dt
    elmts[1].text=txt1

I am using Python 2.7, ArcGIS 10.2, Notepad++ as text editor and pythonWin IDE.


texts: Moving point A to point B with attributes using QGIS

texts:
I have an A layer containing points with a XX attribute.
I would like to move (on points A) the points of a layer B with the same XX attribute as A. There is no notion of distances (from a few cm to several kilometers) and I can't find the right solution with QGIS.


texts:
This question was posted in stackexchange, where there were no views/comments/answers.

I have a raster file and would like to take a stratified sample.

library(raster)
f = system.file("external/test.grd", package="raster")
r = raster(f)
r[r < 500] = 1
r[r >=500 & r <=1200] = 2
r[r > 1200] = 3

plot(r, legend = F)
s = sampleStratified(r, 20, sp = T)
points(s, pch = "+")

This has 2 problems for my intended use. First, the same amount of points are assigned to each area, although the surfaces differ a lot. Second, the points are clustered, and I would like to see a minimum distance between them.

I can solve the first problem by isolating individual layers and assigning a number of samples that is proportional to the area:

r[r != 1] = NA
plot(r, legend = F)
s = sampleStratified(r, 20, sp = T)
points(s, pch = "+")

But that does not solve the clustering problem.

How can I use raster to take a stratified sample that is proportional to the areas and in which some distance between sampling points is maintained?


texts:
I have this paper. In equation 2 they are calculating surface emissivity.
I want to ask what is  v plant proportion?

https://www.researchgate.net/file.PostFileLoader.html?id=5760fe20b0366df3a5237081&assetKey=AS%3A373137431646209%401465974304318

 = **  + [ −]  +
  is the vegetation proportion in a pixel

Was it calculated from NDVI?


texts:
So I am trying to do a where clause on a select and my code is

ST_Within(shape::GEOMETRY, ST_Buffer(ST_GeogFromText('@Value(_location)'), 5, 'endcap=round join=round')

I am trying to buffer the _location attribute by 5 meters so when the ST_Within is used it will help match better.  This code is suppose to look for for input points and see if they match anything in the data set already. via spatial match.  Just FYI both points for ST_within are points.

If anyone has any idea how i used ST_Bufffer wrong that would be great help.


texts: How to Split building outlines polygon into individual features/buildings?

texts:
Like the title suggests I have ONE Building Outlines Polygon where instead of all the buildings being split into their own features they are all in one feature. How can I split these into individual features? I would like all of them to be their own entry in the attribute table.

I can't seem to find anything.  The picture below shows my situation. Any help would be appreciated!


texts:
I am currently working on a project where I try to reconstruct historical borders of the austro-hungarian empire between 1815 and 1918. Therfore I use the comination of digitized old maps an legal texts, as well as already existing shapefiles.

I would like to create a homogenous dataset and afterwords split it to the different time spans where administrative changes took place.

For some regions and timespans we use already existing shapefiles. E.g. from HGIS or Census Mosaic Project (https://censusmosaic.demog.berkeley.edu).

The problem is now that the polygons won't suit together. I fight with great inaccuracies, inconsistencies and data heterogenity. The files are mismatching and contraicting.

At the moment I revise the polygons manually, but I am looking for a way to adapt and harmonice the data more automatized. Creat differences are especially at the external borders, as you can see in the picture below.

Does anybody know a way how I can deal with this problem or a tool i could use to match the borders and the polygons?

I work with QGIS and do not have a lot of programming/Python skills though.


texts: Kernel density estimate with physical boundary

texts:
I'm am using adehabitatHR to estimate the size of the foraging range of seabirds using telemetry data. Is there anyway to take into account a physical boundary, where it is impossible for the animals to travel? In my case, the seabirds do not fly over land, therefore the tortuous coastline and islands are not 'available habitat'. Obviously, a simple kernel density estimate does not consider this, and I end up with a polygon that overlaps land.

this is my current code for the home range estimate. Is there anyway I can take into account land, using a polygon file? (I have a polygon shapefile called coastline)

dfrelocations is a dataframe containing the coordinates of animal relocations

relocations_xy <- as.data.frame(dfrelocations[,c("longitude","latitude")])
relocations_xysp <- SpatialPoints(relocations_xy)
proj4string(relocations_xysp) <- CRS("+proj=longlat +datum=WGS84")
relocations_xysp <- spTransform(relocations_xysp, CRS("+proj=utm +zone=9 ellps=WGS84"))

bird_ud <- kernelUD(relocations_xysp, h=5000)
bird_homerange <- getverticeshr(bird_ud, 95)


texts: What do the following eigenvalues/eigenvectors tell me regarding the bands in the following LANDSAT TM image?


texts: Plotting lines graduated based on a speed attribute

texts:
Is there a way of turning point data (ship vessel tracks.csv) into lines, whereby the lines inherit a attribute (speed)? I have tried points to one plugin which doesn't carry over any attributes from the original point data.

The goal is to show vessels operating at a set speeds e.g 0-4 not trawling, 4-6 trawling, >6 steaming as tracks(lines) with each speed category using a different coloured line.

The end shapefile would be a continuos line with different colors indicating the speed of the vessel.

Im using QGIS and ARCGIS so a solution for either platform would be great.


texts:
I have a few choropleth maps in one Leaflet map using folium. However, in my output, my legends all congregate underneath one another and on one side where it blocks the 'control zoom' icon.

Is there a way to alter the location of my legend(s) or only have one activated at a time?

Image of result:

Code of

NYC_coords = (40.7, -73.8) # specifying our location data
flm_map = flm.Map(location = NYC_coords, zoom_start = 10.3, tiles = 'cartodbpositron', zoom_control = True)# zoom level, tiles
dangerous_int = flm.FeatureGroup(name='Dangerous Intersections', control=True)
neighborhoods = flm.FeatureGroup(name = 'Neighborhoods', control = True)

dscore = flm.Choropleth(
        geo_data = nbhood, # json file, contains geometries
        data = nbhood_df, # dataframe used
        name = 'Danger Score',
        columns = ['Neighborhood','Danger Score Avg'], # columns to refer to for choropleth colors
        key_on = 'feature.properties.Neighborhood', # the key we are looking for to base
        fill_color = 'Reds', # color ramp
        fill_opacity = 0.7,
        line_opacity = 0.5,
        highlight = True,
        popup = flm.GeoJsonPopup(fields = ['Neighborhood', 'Danger Score Avg']),
        legend_name = 'Average Danger Score',
        smooth_factor = 0.2).add_to(flm_map)
motor_inj = flm.Choropleth(
        geo_data = nbhood, # json file, contains geometries
        data = nbhood_df, # dataframe used
        name = 'Motor Injuries',
        columns = ['Neighborhood','Motor Vehicle Injuries'], # columns to refer to for choropleth colors
        key_on = 'feature.properties.Neighborhood', # the key we are looking for to base
        fill_color = 'BuPu', # color ramp
        fill_opacity = 0.7,
        line_opacity = 0.5,
        popup = flm.GeoJsonTooltip(fields = ['Neighborhood', 'Motor Vehicle Injuries']),
        highlight = True,
        legend_name = 'Motor Vehicle Injuries(total)',
        smooth_factpr = 0.2,
        show = False).add_to(flm_map)

flm.LayerControl('topright', collapsed = True).add_to(flm_map)
display(flm_map)


texts: Snapping line to polygon border using PyQGIS?

texts:
I want to snap the line (completely present inside the polygon) to polygon boundary. I tried create 10 m buffer on line start and end point to get polygon and check intersection between line and polygon it returns all the line vertices as intersected. I want only intersection of line to polygon boundary. I use the following code to get an intersection. My aim is if line is 10 m distance to border snap to line border.

 intersection = line.geometry().intersection(polygon.geometry())
 intersectstr= intersection.exportToWkt()
 gem=QgsGeometry.fromWkt(intersectstr)
 intersetpnt=gem.asPoint()


texts: PostgreSQL password setup in QGIS (disallowing multiple authentication to individual databases

texts: I have a QGIS project where I am connected to about 10 databases. The issue I have is that I cannot seem to save the password to disallow repetitive password authentification. After logging in individual databases and working on the project when I pan around a request comes to log in after a few seconds. I have seen somewhere it suggest that I save the password when I do the connection to the database, however in QGIS 3.6, this seems to be different or I cannot get it right, after storing credentials, and not working on the project for a few minutes, when you come back it will ask you for the password again. This link here suggest it is a bug? From QGIS 3.4.5 and 3.6 asks password repeadly to connect to PostgreSQL 10 .


texts:
Has anyone ever got the cutline to work in C++?

I have it working off the command line no problem, but it doesn't seem to do anything when I populate that argument in C++ as below:

*//Create a OGRGeometry variable
OGRGeometry *thisGeom;
//parse an input WKT to char*
char* wktIn = (char*) pIn_wkt_geom;
//Create and polulate thisGeom with the input WKT geometry
OGRErr err = OGRGeometryFactory::createFromWkt(&wktIn, poSRS, &thisGeom);

//Add thisGeom as the cutline argument psWarpOptions->hCutline = thisGeom;*


texts: How to identify Sliver triangles in a TIN surface and count them?

texts:
Is it possible to count slivers in a TIN?

I generated a Delaunay triangulation surface model and i want to know how many slivers the TIN has.


texts: Qgis: QgsLogMessage.logMessage to output message in a specific tab

texts:
I am writing a Qgis plugin.
I would like the plugin to write log message in a specific tab, in the Log Messages window. I tried:

QgsMessageLog.logMessage("File "+ii.source()+": all pixels set to no data", 'Raster Processing',QgsMessageLog.INFO )

but QGis keeps on writing under tab "General". I suppose I'm missing an initialisation part to instantiate a new tab and connect it to the message logger?


texts: Adding constrain/check after adding first point of line and after last point of line

texts:
I want to customize the QGIS editing tool in the following way.

Scenario: The first and last point of the line layer must have a point layer.

Steps:

User selects a line layer.
User starts/toggles the line editing.
When user starts clicking on the map (first point) of line there should be a check that it must be snapped to the one of the point layer point loaded on the QGIS if not then message must be displayed that no point found at the start of the line.

I tried the following:

signal geometryChanged works on existing geometry
signal featureAdded works after the geometry and attribute is added
signal committedFeaturesAdded works after the edited geom and attribute is commmited
class qgis.gui.QgsMapToolEdit does not have signal for geometry add trigger
cannot implement multiple MapTool as one canvas can have only one active maptool
qgis.gui.QgsMapCanvas have mouseReleaseEvent event but I cannot figure out how to use it as its not a signal

I cannot find a signal that triggers when the actual geometry editing starts and in the signal I want to receive the coordinates of the mouse release point.

I am planing to implement the mouse release event using PyQt5.QtWidgets.QGraphicsView its the base class of qgis.gui.QgsMapCanvas any suggestion or direction is appreciated.