texts: Why do I get negative values for the topographic position index (TPI)?

texts:
I run into issues, when trying to calculate the topographic position index (TPI).
I get values ranging from -13 to 17. And the resulting layer looks a bit odd, too:

If I understand it correctly, this index should range from 1-10 (Poster of TPI-creator)
My DEM as a resolution of 30m and I used the default values fo the Tpi-Tool.
I'm running QGIS 3.16.3 in Windows 10 and I tried to use "SAGA > Terrain Analysis - morphometry > Topographic position index (tpi)" from the toolbox. (I also tried "GDAL > Raster analysis > Topographic Position Index (TPI)", resulting in values ranging from -80-70ish).
I made sure I have the right projection assigned to my DEM.
I also made sure there are no holes or anything in the DEM.
It gives me the error message:
The following layers were not correctly generated.
• C:/Users/xxxx/AppData/Local/Temp/processing_flYZSV/1abd3e1d82594e88b04d28cd346343e3/TPI.sdat
You can check the 'Log Messages Panel' in QGIS main window to find more information about the execution of the algorithm.

However, in the log-panel I can't find something helpful:
021-04-16T17:31:35     INFO    SAGA execution commands
   ta_morphometry "Topographic Position Index (TPI)" -DEM "C:/Users/xxxxx/AppData/Local/Temp/processing_flYZSV/10afe001939746859ef23972335147d5/DEM30mSavaSofia.sgrd" -STANDARD true -RADIUS_MIN 0.0 -RADIUS_MAX 100.0 -DW_WEIGHTING 0 -DW_IDW_POWER 1 -DW_IDW_OFFSET true -DW_BANDWIDTH 75.0 -TPI "C:/Users/xxxx/AppData/Local/Temp/processing_flYZSV/642542ca03774ce59275236d376a8fc0/TPI.sdat"
2021-04-16T17:32:00     INFO    SAGA execution console output

   D:\xxxx>set SAGA=C:/OSGEO4~1/apps\saga-ltr
   D:\xxxx>set SAGA_MLB=C:/OSGEO4~1/apps\saga-ltr\modules
       D:\xxxx>PATH=C:\OSGEO4~1\apps\qgis\bin;C:\OSGEO4~1\apps\Python37;C:\OSGEO4~1\apps\Python37\Scripts;C:\OSGEO4~1\apps\qt5\bin;C:\OSGEO4~1\bin;C:\WINDOWS\system32;C:\WINDOWS;C:\WINDOWS\system32\WBem;C:\OSGEO4~1\apps\Python37\lib\site-packages\pywin32_system32;C:\OSGEO4~1\apps\Python37\lib\site-packages\numpy\.libs;;C:/OSGEO4~1/apps\saga-ltr;C:/OSGEO4~1/apps\saga-ltr\modules

   D:\xxxx>call saga_cmd ta_morphometry "Topographic Position Index (TPI)" -DEM "C:/Users/xxxx/AppData/Local/Temp/processing_flYZSV/10afe001939746859ef23972335147d5/DEM30mSavaSofia.sgrd" -STANDARD true -RADIUS_MIN 0.0 -RADIUS_MAX 100.0 -DW_WEIGHTING 0 -DW_IDW_POWER 1 -DW_IDW_OFFSET true -DW_BANDWIDTH 75.0 -TPI "C:/Users/xxxx/AppData/Local/Temp/processing_flYZSV/642542ca03774ce59275236d376a8fc0/TPI.sdat"
   ____________________________

   ##### ## ##### ##
   ### ### ## ###
   ### # ## ## #### # ##
   ### ##### ## # #####
   ##### # ## ##### # ##
   ____________________________

   SAGA Version: 2.3.2 (64 bit)
   ____________________________
   library path: C:\OSGEO4~1\apps\saga-ltr\modules\
   library name: ta_morphometry
   library : Morphometry
   tool : Topographic Position Index (TPI)
   author : O.Conrad (c) 2011
   processors : 8 [8]
   ____________________________

   Load grid: C:/Users/xxxx/AppData/Local/Temp/processing_flYZSV/10afe001939746859ef23972335147d5/DEM30mSavaSofia.sgrd...

   Parameters

   Grid system: 30.284487; 5718x 5048y; 231561.120511x 8322967.608079y
   Elevation: DEM30mSavaSofia
   Topographic Position Index: Topographic Position Index
   Standardize: yes
   Radius: 0.000000; 100.000000
   Weighting Function: no distance weighting

   D:\xxxx>exit

So I wonder, what could be the cause? And what else could I try?
PS: I'm sorry I don't know how to attach sample data, otherwise I would do.


texts:
Is there a way to export my customized settings or copy some xml or config file, in order to install my settings to another computer without me manually changing everything? If this is possible it would make the installation process much easier for people who are not that familier with GIS.

Many people at my workplace are going to be using Qgis in the future and this would make the transition easier :)


texts:
I am trying to do kriging on mortality rates using QGIS (I have version 3.0 but I downloaded the 2.18.20 to be able to follow the web documentation and most of the answers here).

From what I have understood, before running kriging it is necessary to run a variogram model to get the estimated model. I know both variogram and kriging can be run in QGIS using SAGA tools and that is what I am trying to do.

However, most of the "tutorial" I found online provide an empirical variogram graph (y-axis, semivariogram and x-axis, distance) that I am not able to draw with any command on QGIS.

Should I run this part of the analyisi in R as most of the people do or there is a way to do the same staying in QGIS?


texts: Finding vertices of polygons for calculating offsets in R

texts:
I am working on an application where a user can select a land plot (sf_POLYGON) by clicking on it. Once selected (highlighted plot no. 6 in the figure), I then need to calculate some "offsets" and plot the resulting polygon. These offsets are dependent on the adjacent polygons. For example, if the selected plot has a road adjacent to it, then the offset is say 6m inside the plot. The side adjacent to a plot will have an offset of 4m. The same is displayed in the inscribed polygon. The blue lines are for sides touching the road and the orange lines are for the sides touching other plots.
The way I was thinking to solve this was :

Find the vertices of the selected polygon.
Draw lines between the vertices.
Use st_intersect or st_touches to find out the polygons that the line touches.

I tried to convert the polygon into LINESTRING using st_cast(test_poly, 'LINESTRING') to figure out if I can get the lines automatically and then see which polygons they intersect/touch. However st_coordinates on the polygon cast to LINESTRING reveals coordinates which are not all vertices. So my question is, how do I extract the vertices of the selected polygon? Or if there is a better approach to solve this problem?
The second problem is of creating the offset polygon. Not sure how to do that, but i am guessing buffer is something that would help?
library(sf)
library(leaflet)
library(dplyr)[![enter image description here][1]][1]
test_data <- readRDS('test_data.Rds')
leaflet() %>%
  addTiles() %>%
  addPolygons(data = test_data, fillOpacity = 0.8, stroke = TRUE,
              weight = 1,
              color = 'black',
              label = test_data$Plot_Numbe,
              layerId = test_data$OBJECTID,
              fillColor = test_data$color)

test_poly <- test_data %>% filter(OBJECTID == 1977)
test_line <- st_cast(test_poly, 'LINESTRING')
st_coordinates(test_line)

Can't seem to attach the data. Here's a link to it.


texts: Are there free maps available that can show me land planted with oil palm?

texts: Currently I am taking gps coordinates of oil palm plantations in Africa. I'm putting the waypoints in a GIS program but the free basemap provided isn't always giving high resolution images, making it difficult to adjust the waypoints to the real boundaries of a plantation. On Google Earth you can detect plantations easily. So I was wondering which free imagery is available to detect plantation boundaries/palm trees. Or is there a way to detect oil palm plantations from low resolution imagery (landsat).


texts:
Is there a way to prefer straight roads in network analyst rather than going through zigzag routes (i.e creating straight forward routing instead of really going through the shortest path)?

I am creating a wayfinding map and want to keep directions straight.

In this example, I want to go from point 1 to point 2 but I don't want to route through the diagonal road(cross). I am not allowed to manipulate time on the network and I cannot block the diagonal road. I prefer to use the other road (with more vertices and longer) is there a way to do this?


texts:
I already tried this http://hg.orfeo-toolbox.org/OTB-Qgis-plugins/ repository and couldn't get  it...
This guy seems to have it: http://www.elance.com/samples/orfeo-toolbox-qgis/24111664/


texts:
I'm having some trouble parsing the following WFS response into a feature:

<?xml version="1.0" encoding="UTF-8"?>
<wfs:FeatureCollection numberOfFeatures="1355" xsi:schemaLocation="http://www.emsa.europa.eu/ais ./STIRES_WFS.xsd http://www.opengis.net/wfs http://schemas.opengis.net/wfs/1.1.0/wfs.xsd" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:ais="http://www.emsa.europa.eu/ais" xmlns:gml="http://www.opengis.net/gml" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns:ows="http://www.opengis.net/ows/1.1" xmlns:wfs="http://www.opengis.net/wfs">
    <gml:featureMember>
        <ais:feature>
            <ais:MMSI>219009397</ais:MMSI>
            <ais:callsign>unknown</ais:callsign>
            <ais:name>unknown</ais:name>
            <ais:trackProperty>
                <ais:track>
                    <ais:objectStatus>
                        <gml:Point srsName="EPSG:4326">
                            <gml:pos>8.422975 55.475732</gml:pos>
                        </gml:Point>
                        <ais:time>2013-09-18T23:42:34.000</ais:time>
                        <ais:speedOverGround>0</ais:speedOverGround>
                        <ais:courseOverGround>306</ais:courseOverGround>
                        <ais:heading>244</ais:heading>
                        <ais:navigationalStatus>7</ais:navigationalStatus>
                        <ais:orderDetailId>127889</ais:orderDetailId>
                    </ais:objectStatus>
                </ais:track>
            </ais:trackProperty>
        </ais:feature>
    </gml:featureMember>
</wfs:FeatureCollection>

Everything works fine up to the trackProperty, which is undefined in the parsed feature's properties. The rest of the properties (MMSI, callsign, name) all show up fine. So I assume the problem is with nested properties.


texts: Processing Modeller manual / use-case for Multiple Input and Model Variables

texts:
Is there a decent manual/examples available somewhere for the Processing Modeller?

Official help is very barebones and doesn't explain the different inputs/model variables.

Specifically I am interested in a use-case for: Multiple Input and Model Variables.


texts: Why can't I save new points created in QField?

texts:
We are using QField 0.7.5 that seems a really great project to bring QGis on a mobile device.

I prepared a QGis 2016 project : orthophoto + 1 polygon layer (casdatre) + 1 point layer to input data.
I transferred the project and the data to a Galaxy Tab A (Android 5.0.2).
I open the project with QField and everything displays fine, I can move and query the data, really great.

But when I create a point, I can't save the new point (save don’t save the new point).

Do you have any idea what’s wrong?


texts: Converting point layer to raster grid that shows frequency of points per cell uisng QGIS

texts:
I have a point layer that shows the distribution of a species.

Is it possible in QGIS to convert this point layer to a raster grid in which the value of each cell corresponds with the amount of points within that cell?

Until now, I only managed to attach attribute values to the grid cells. I have already added a column in my attribute table that has the value '1' for each point, hoping that there is a way to use a sum of attributes for each grid cell.


texts: Calculate location on line with known azimuth to point

texts:
I have the following situation (consider distances up to 300m - without great circle etc. affects):

How can I calculate the unknown location when I have all other locations in the triangle and azimuth to the unknown location from some source point (green)?


texts: How to execute ol.extent button from a new button in OpenLayers 3?

texts:
I'd like to move the ol.extent button that it shows in div.map to new button in another div.

I have the following code but without success:

    Set Position

Any idea?


texts:
I am moving an object from one point to another and would like to update the heading of the dynamic display properties. I am using this code:

m_dynamicSymbolProperties.set_Heading(esriDynamicSymbolType.esriDSymbolMarker, heading);

But I am having trouble getting my heading right.  Does anyone know how to calculate a heading based on two IPoints which have standard esriSRGeoCS_WGS1984 spatial references?


texts: SQL Server - ESPG:4326 data stored in Geometry column

texts:
It seems that the Microsoft SQL Server Geography data type is better suited for lat/lon because Geography "stores ellipsoidal (round-earth) data, such as GPS latitude and longitude coordinates."

While the Geometry data type "represents data in a Euclidean (flat) coordinate system."

I have data stored in SQL Server Geometry data type column, and it is lat/lon coordinates - the SRID is EPSG:4326.

The reason being, is the data is used by Geoserver, and Geoserver only recognizes Geometry data types for SQL Server. But the data is 4326.

Is there anything wrong with doing this?  Will there be accuracy or performance issue?  Should/could the data be converted?


texts: How to split lines with defined split number by user in pyqgis?

texts:
Actually my question is on header. I want to split (or explode whatever you say) every line in my lineshape layer. For every line, I must define the split number and it should split them with equal length according to these numbers. Is it possible in pyqgis? I am working with QGIS 3.0.

Here is my layer:

After defining numbers it should be split like that (I wrote some of them split numbers and expected explosion:


texts: Leaflet + GeoServer + British National Grid

texts:
I have followed along some examples and was able to display BNG 27700 data nicely on the map, however now I want to perform getFeatureInfo function on click. Again the code I used is as below:

var BBOX=map.getBounds().toBBoxString();
var WIDTH=map.getSize().x;
var HEIGHT=map.getSize().y;
var X=map.layerPointToContainerPoint(e.layerPoint).x;
var Y=map.layerPointToContainerPoint(e.layerPoint).y;
var SRS=map.options.crs.code;
var URL="http://localhost:8080/geoserver/OSM/wms?INFO_FORMAT=application/json&REQUEST=GetFeatureInfo&EXCEPTIONS=application/vnd.ogc.se_xml&SERVICE=WMS&VERSION=1.1.1&WIDTH="+ WIDTH+"&HEIGHT="+HEIGHT+"&X="+X+"&Y="+Y+"&BBOX="+BBOX+"&SRS"+SRS+"&LAYERS=OSM:BuildingBoundaries&QUERY_LAYERS=OSM:BuildingBoundaries";`

However no feature is returned. If I make use of sample data from Geoserver in leaflet, the click event works perfectly. Am I missing something?


texts: Field calculator returning ERROR 000989?

texts:
I'm trying to split a field of street addresses (field name "PROP_STREE") so that it returns only the house number (so "123 Main St" would become "123").

Following the solution provided Field Calculator: Split address into house number, street name, I am using the following script:

def addressParser(inString):
  splitString = str(inString).split(' ')
  houseNumber = splitString[0]
  return houseNumber

PROP_NUM = addressParser(!PROP_STREE!)

But it keeps returning "Failed to execute. Parameters are not valid. ERROR 000989: Python syntax error: Parsing error SyntaxError: invalid syntax (line 6)."

What mistake am I making in my code?

Here's what the Field Calculator window looks like, for reference.

I'm using ArcGIS 10.6.1.


texts:
I'm working with an inventory of a species in Southern Sweden. I have recorded my findings (together with coordinates, name of the species and year of finding) and created a CSV-file with i have imported to QGIS.

Example from my data

Decade;species;ort;lat,long
1880;falcata;Borgholm;56.882657,16.655931
1880;falcata;Algutsrum;56.674716,16.53471
1880;falcata;Borgholm;56.882657,16.655931

In order to make a map that's easy to grasp I would like to collect all findings from e.g 1890 into an area, all findings from 1900 into a second area an so on. Someone that have a clue on what I can to to achieve this?


texts: SQL Spatial view using a Query Layer in ArcMap 10 is SLOW

texts:
We have a number of SQL Spatial views each made up of approximately 10 tables (so relatively complex views) and each including a table with a SQL Geometry column. A view could return up to 300-400 thousand records. No SDE involved. We are using a Query Layer (QL) for each spatial view (Select * from the view) to display in ArcMap. Each spatial table has a primary key, spatial grids and uses an SRID of 0. Just to load a QL might take 10-15 seconds. To display with a 100k - 250k data view (maybe 100s of features) takes 10-20 seconds. We previous had these as SDE views (same underlying SQL views w/o the SQL spatial column) in Arc 10 and they loaded and displayed in 1-3 seconds.

Anyone have ideas about why these are SO MUCH SLOWER without using SDE?


texts: Openlayers: authenticate users to send any wfs requests to geoserver

texts:
I can see that this question is asked before but not answered properly.
I'm trying to authenticate with geoserver before trying to send any wfs requests.
I need to use the geoserver users.

I already tried this and the login was succesfull but i can't find a way to check if the login is succesfull or not cause the cookie that geoserver returns is http-only.

Any ideas?


texts: Within returns results despite there is no points inside the selected polygon

texts:
I have a method that uses within function, in order to check if a point is within a polygon:

public List<PointsEntity> findPointsIntoAnExistingPolygon(PolygonEntity polygon) {
                StringBuffer sb = new StringBuffer();       sb.append("select p from PointsEntity p, PolygonEntity poly "       + "where within(p.coordinates, poly.coordinates) = true");
                Query q = em.createQuery(sb.toString());
                        List<PointsEntity> pointFound = new ArrayList<PointsEntity>();
                try {           pointFound =  (List<PointsEntity>) q.getResultList();           }catch(NoResultException nre){          }
                return pointFound;  }

But it returns all the points in the DB that there are into a polygon.

For example, if I pass as a parameter a polygon A that doesn´t contain a point, the method returns the point that is contained by polygon B instead of return nothing.


texts: Stack ImageColleciton in a multi-band Image and name fields after date in Earth Engine

texts:
In Earth Engine how can the bands generated through a stacking of an ImageCollection to a a multi-band image be named after the corresponding month-year?

I use the following code to stack the ImageCollection, but the bands are simply named sequentially with integer numbers.

def stackCollection(collection):
first = ee.Image(collection.first()).select([])
def appendBands(image, previous):
    return ee.Image(previous).addBands(image)
return ee.Image(collection.iterate(appendBands, first))

stacked = stackCollection(imageCollection)


texts: PostgreSQL : finding all lines intersecting with a polygon

texts:
I am totally new in the world of GIS but as a student surveyor I am working on a project using postgreSQL and PostGIS.

I have one database with all the roads of my country (lines), and one database with a lot of area's (polygons). And for every polygon, I need to find all the roads (lines) that intersect with it (contain in partially or fully).

I do know the ST_Intersects function, but is this the best function to use ? right now I am trying to make a loop function where I loop all the roads to see if they intersect, but there must be a quicker way..

If anybody could get me started, that would be great!

Thanks


texts: Assign additional style to GeoServer GetMap request to view those layers in CKAN

texts:
I am using GeoServer 2.12.x & CKAN 2.8 on Cent-OS.
I have added GeoServer layers in Ckan. I am using plan Getmap request to add layers in Ckan
GetMap request : https://domain.name/geoserver/workspace/layer/ows?service=WMS&version=1.1.0&request=GetMap
I have layer called 'US State' for which I have assigned 3 SLD's in Geo-server. I am passing "style-2" using additional parameters such as "styles", "SLD_BODY" in URL its giving me error "Error parsing style" only white space  content is allowed.
Passing sld name as follows : https://domain.name/geoserver/workspace/layer/ows?service=WMS&version=1.1.0&request=GetMap&styles=workspace%3Alayername

In GeoServer layer preview URL I can change style easily & its working but the thing is that GeoServer layer preview URL is not working in CKAN
I researched on my side a lot, but I am not getting any specific direction/hint. I have used 'Geo-view' extension in CKAN for display GeoServer layers.


texts: Converting from 16-bit TIFF image to 8-bit with GDAL does not work properly

texts:
I am a beginner in the area of GIS and remote sensing image analysis. I am dealing with the satellite images from the **VNP46A1 Daily At-sensor TOA Nighttime Radiance Product. I want to get the TIFF images from the DNB_At_Sensor_Radiance_500m scientific dataset layer for further analysis.
I used the HDF to GeoTIFF convert script from NASA which I slightly changed. I need to reproject the images from the WGS84 to the S-JTSK coordinate system. The original arrays from which the images are produced contain UInt16 values. When I want to display the image, I got the black image. I need to convert them to the 8-bit images. In order to do this, I used the parameter scaleParams from gdal.TranslateOptions. However, I am still seeing the black image. Below is the part of the code that I used:
gdal.TranslateOptions(gdal.ParseCommandLine(translateoptions),
 format='GTiff', height=2400, width=2400, scaleParams=[[0, 65535, 0, 255]], resampleAlg=gdal.GRA_Bilinear, outputType=gdal.GDT_Byte)

How can I solve this problem?


texts: Seeking free application or web service for reverse geocoding

texts:
I'm looking for a fast and free web service or local application (this it would be better, a local service with database of street, numbers, city, etc.) for reverse geocoding.

I could use google but is limited to 5000 query per day...

What I m trying to do is convert positions acquired from a gps (1 per second) in addresses. These positions could be more than 7000 per track, and if I have 10 track for example, the mole of data to "translate" becames "huge"..


texts: QGIS Georeferencer; CRS selector doesn't open, destination coords 0

texts:
I've tried this in QGIS 2.14 LTR, 2.16 and 2.18, and a fresh install of 2.14 LTR on another PC.

My questions are:

How do I get the CRS selector dialog to open? When opening a raster image, I get no CRS selector dialog, as these two guides clearly show it should. Instead I get an internal error message "CRS was undefined: defaulting to project CRS". I can select the right CRS from the project properties, but the second issue is unchanged.
How do I make the georeferencer output destination coordinates? When adding control points in the Georeferencer, their coordinates just show 0.00. Predictably, the image will fail transformation.


texts:
Whenever I go to the Esri website to download ArcGIS Desktop, it always directs me to download ArcGIS Pro.

Can someone provide me with the Esri link to the ArcGIS Desktop download?


texts: Why is SQL input missing features from original Origin-Destination Matrix using QGIS?

texts:
I am using QGIS 3.16.
I computed an Origin Destination Matrix for my studies using a total of 2478 (x13) routes. I then computed an SQL query in order to match each start point with its nearest end point. But this only computed 2458 routes (when it should have computed 2478).
Does anyone know why this is the case or what may have happened?
My SQL input is below. The algorhitm's ran smoothly and no error messages appeared, so I'm not quite sure why it cut layers out.


texts: Converting 32639 spatial reference to 3857 in PostGIS


texts: Which function/condition should I use to get the opposite of "PropertyIsLike" in GeoServer?


texts: SQL for spatial view of line start/end points (Oracle ST_GEOMETRY)

texts:
How do I create a spatial view that generates dynamic point geometry of the start point and end point of a line feature class?

I have successfully made a similar view that is a buffer of a line feature class with this SQL:

SELECT OBJECTID, SDE.ST_BUFFER(SHAPE, 1) AS LINE_BUFFER FROM ENG.ROAD

I thought I'd be able to use similar SQL to generate start points of lines using the PointN  function:

SELECT OBJECTID, SDE.ST_POINTN(SHAPE, 1) AS START_POINT FROM ENG.ROAD

The view creates successfully, however it only creates a tabular view, not a spatial view. What have I done wrong?

I'm also unsure how to get the end point index position, but I'm not even thinking about that yet.


texts: If do not want to use pgRouting then what?

texts:
I am interested in creating an online routing application. Supposing I have all the data I want into a postGIS database and loaded onto GeoServer and I am willing to use OpenLayers for front end, and of course I know routing algorithms like Djikstra or shortest path algorithms, what are my choices when it comes to do the actual routing.

I know many use pgRouting, and as many members of this community have told me I tried to learn it and use it. But what I want in my web app is more freedom. I want to change the algorithms and be able to whatever I imagine(academic reasearch-change the algorithms at some points).

So what are my choices ?

SQL user defined functions?

How to manipulate the data and find shortest paths?

In few words how can I create a library/extension just like pgRouting?

What languages do I need(if any) and other material?

Any good tutorials?


texts: Converting .txt containing two columns ID and GeoJSON Geometry?

texts:
I am working with an unusual format file - it came as a .txt and is separated by tabs.  In the first column is an id number and in the second column is the geometry stores as geoJSON.

When I attempt to import into QGIS as a text layer, breaking it up by tabs and setting the header works.  When I specify the geometry as the geoJSON field I get the error that it is an invalid file.  I am not sure how to approach converting this file into a Shapefile (ideal solution) or something that QGIS can display.

Id      geoJSON
49518   {"type":"Polygon","coordinates":[[[29.9596,-2.3272],[29.9189,-2.7028],[29.7235,-2.8189],[29.3706,-2.8401],[29.1497,-2.5917],[29.0405,-2.7447],[28.8619,-2.5309],[29.1194,-2.2488],[29.1361,-1.8602],[29.3616,-1.5092],[29.8237,-1.3089],[29.9149,-1.4825],[30.0524,-1.4314],[30.4704,-1.0525],[30.8394,-1.6514],[30.8567,-2.3151],[30.1356,-2.4382],[29.9596,-2.3272]]]}


texts:
I'm trying to write some code to skip over empty raster sets when converting rasters to point (the empty sets stop the conversion, leaving me with a half complete model run).

I know that I can identify empty rasters using the ALLNODATA result from "Get Raster Properties" Create conditional to check if raster has value.

I know I can delete rasters using the Delete.management tool.

I know I should be able to create a branch using calculate value.

But I'm inexperienced with python and modelbuilder, and haven't been able to put this all together. Could someone help me figure out what I am doing wrong? My code and an image of the modelbuilder setup are below.

def fn(Property):

#Import modules
     import arcpy

#ifelse statement
if Property == "1":
           return "true"
elif Property == "0":
           return "false"


texts: How to retrieve the corners of the bounding box from ST_Extent?

texts:
Is there any way to fetch the lat/long of the corners of the bounding box we got from ST_Extent?

select ST_Extent(geom) from tableName;


texts: How to join hydrological data from a big area with a small area with better pixel resolution?

texts:
I use TAUDEM to make hydrologic analysis.
I have a hidro 1K dem file from South America (big area) pixel resolution=1000x1000 meters.
I have a specific dem file located in Brazil (small area) pixel resolution=30x30 meters.
There is no problem to calculate contributing area, slope, streams, etc for each dem file.
My problem is that I need to calculate the small area including data from the big area, for example the contributing area and the stream upslope length.
What is the best way to do this?


texts: Seconds of the day to Seconds of the GPS week

texts: I have a UTC time as 3hr 21min 50sec and this refer 12110 seconds of the day. I want to convert this time to seconds of the GPS week but I'm not sure how to do this. I want to calculations for 2nd of May for 2021 years. Should I apply leap seconds correction for GPS time? Can you share the formula for how to implement the transformation?


texts: Adding elevation from geoserver in Cesium.JS

texts:
I'm trying to view elevation data loaded from Geoserver on cesium globe. But the added elevation is viewed as 2D layer on the globe.

My Code:

import Viewer from "cesium/Source/Widgets/Viewer/Viewer";
import WebMapServiceImageryProvider from "cesium/Source/Scene/WebMapServiceImageryProvider";
this.viewer = new Viewer(this.cesiumContainer);

var provider = new WebMapServiceImageryProvider({
  url : 'http://localhost:8080/geoserver/nyc_roads/wms?',
  layers : 'nyc_roads:DEM_Maharashtra',
  parameters : {
    transparent : 'true',
    format : 'image/png'
}
});
this.viewer.imageryLayers.addImageryProvider(provider);

Do I need to change the function to view it on 3D. Please help.


texts: How to measure correct distances in a project with OpenLayers background?

texts:
I have 3 lat/long points as follows:

A = {35.257246,139.721256}
B = {35.256541,139.721491}
C = {35.256326,139.7203271}

1) Measurement result by  Movable Type Scripts:

By using the distance calculation method at Movable Type Scripts, the distance results are:

distance(A,B)= 80.88745771063131(m)
distance(B,C)= 101.58773310547808(m)

I believe these results are correct.

2) Measurement result by QGIS measure line tool:

The QGIS desktop 2.0.1 and the Projection CRS is set to [+proj=longlat +datum=WGS84 +no_defs]. I use the Plugins/OpenLayers plugin/Add Google Satelite and import 3 lat/long points into QGIS.

The results are:

distance(A,B) = 100.374(m)
distance(B,C) = 135.114(m)

I need to generate a map from a large set of lat/long points as above and Google Satellite background on QGIS.

Could you show me how to solve this problem?


texts: Mapnik: mapnik.has_pycairo() returns False. How to fix this on Windows for Python27?

texts:
I'm trying to modify rundemo.py (comes canned with Mapnik Windows binaries) to render output files at a higher resolution, but it seems I still have issues with my Mapnik-to-Cairo setup.

I've got mapnik and its python bindings working, and I've also got cairo and these python bindings working (..I think). If I start Python, I can do the following without issue:

E:\x_FOSS\mapnik\v2.0.0\mapnik-2.0.1rc0\demo\python>python.exe
Python 2.7.1 (r271:86832, Nov 27 2010, 18:30:46) [MSC v.1500 32 bit (Intel)] on win32
Type "help", "copyright", "credits" or "license" for more information.
>>> import mapnik
>>> import cairo
>>>

Meanwhile, back in rundemo.py line 351, it fails on this condition:

if HAS_PYCAIRO_MODULE and mapnik.has_pycairo():

And this..

print mapnik.has_pycairo() # returns False.

But this conflicts with what I'm reading in the Mapnik blog which says these Windows binaries are built with Cairo backend support.

Does anyone know what conditions have to be in place for mapnik.has_pycairo() to return True?

[Update, are these build, or runtime instructions?]

Per Zero_Qualms' advice, I downloaded the mapnik source and found reference to HAVE_PYCAIRO in %\mapnik-v2.0.1\bindings\python\mapnik_python.cpp. I'm not familiar with C++, so I'm not sure if these are build instructions or runtime instructions. First, line 86 has this block following several #include statements:

#if defined(HAVE_CAIRO) && defined(HAVE_PYCAIRO)
#include <pycairo.h>
static Pycairo_CAPI_t *Pycairo_CAPI;
#endif

Then there are a couple insignificant calls to defined(HAVE_PYCAIRO), followed by this has_pycairo() property at line 347:

// indicator for pycairo support in the python bindings
bool has_pycairo()
{
#if defined(HAVE_CAIRO) && defined(HAVE_PYCAIRO)
    Pycairo_CAPI = (Pycairo_CAPI_t*) PyCObject_Import(const_cast<char *>("cairo"), const_cast<char *>("CAPI"));
    if (Pycairo_CAPI == NULL){
        /*
          Case where pycairo support has been compiled into
          mapnik but at runtime the cairo python module
          is unable to be imported and therefore Pycairo surfaces
          and contexts cannot be passed to mapnik.render()
        */
        return false;
    }
    return true;
#else
    return false;
#endif
}

That block comment seems to imply my mapnik-pycairo linkage just isn't correct, and that it can be fixed. Next I searched for PyCObject_Import and found the following block in this file:

%\mapnik-v2.0.1\bindings\python\python_cairo.cpp:

void register_cairo()
{
    Pycairo_CAPI = (Pycairo_CAPI_t*) PyCObject_Import(const_cast<char *>("cairo"), const_cast<char *>("CAPI"));
    if (Pycairo_CAPI == NULL) return;

    boost::python::converter::registry::insert(&extract_surface, boost::python::type_id<PycairoSurface>());
    boost::python::converter::registry::insert(&extract_context, boost::python::type_id<PycairoContext>());
}

Blast! It seems reasonable that I could create the registry entries manually, but ..what entries?


texts: How to create a .lyr file pointing to a non-existent shapefile?

texts:
I've been asked to write code that creates a .lyr file pointing to a non-existent shapefile.  I expected to be able to do this utilizing IDataLayer2, however I get COMException FDO_E_TABLE_NOT_FOUND when I call dLayer.DataSourceName = (IName)fcn in the c# below.

Update: As suggested by Vince, here's the use case.  We have an end-user with limited experience in ArcGIS who's offered to test an app we're developing.  We'd like the user to verify our app works with their ArcSDE instance.  At this point the app uses pre-configured lyr files that we haven't set up yet.  So instead of trying to explain to the user how to set up the lyr files, I've been asked to provide the user lyr files, that point  their ArcSDE instance.  I'm working off-site, with no access to their network.  I'm assuming that if I can't do this with a shapefile, then I can't do it with an ArcSDE featureclass either, so I've posted code that keeps it simple.

public static void Test()
{
    string folder = @"C:\projects\";

    string shpfile1 = Path.Combine(folder, "mypolygons.shp"); // exists
    string lyrfile1 = Path.Combine(folder, "testsave1.lyr");

    string shpfile2 = Path.Combine(folder, "junk.shp"); // doesn't exist
    string lyrfile2 = Path.Combine(folder, "testsave2.lyr");
    TestSaveLyr(shpfile1, lyrfile1); // works
    TestSaveLyr(shpfile2, lyrfile2); // fails
}

public static void TestSaveLyr(string shpFilepath, string lyrFilepath)
{
    // make a simple featurelayer for a shp file and save it
    FileInfo fi = new FileInfo(shpFilepath);
    var wsn = new WorkspaceNameClass() as IWorkspaceName;
    wsn.PathName = fi.DirectoryName;
    wsn.WorkspaceFactoryProgID = "esriDataSourcesFile.ShapefileWorkspaceFactory";

    var fcn = new FeatureClassNameClass() as IDatasetName2;
    fcn.WorkspaceName = wsn;
    fcn.Name = fi.Name.ToLower().Replace(".shp", "");

    if (File.Exists(shpFilepath))
        Debug.Print("exists: " + fi.Name);

    var dLayer = new FeatureLayerClass() as IDataLayer2;
    ((IFeatureLayer)dLayer).DataSourceType = "Shapefile Feature Class";
    dLayer.DataSourceName = (IName)fcn; // COMException here FDO_E_TABLE_NOT_FOUND -2147220655
    ((ILayer)dLayer).Name = fcn.Name;
    Util.Save((ILayer)dLayer, lyrFilepath, true);
}


texts: Connecting own OSRM server on leaflet routing machine?

texts:
I have my OSRM running on UBUNTU.

How do I change the demo server of routing machine to my own OSRM server on Leaflet?


texts: Want to get an SHP file for my coutry Algeria

texts:
H everyone , i'm running a clustering test on my country with R, but the problem the clustering doen't look good because of the map i have, i hant to get a bigger map of my country 'ALGERIA' in an SHP format!

Can you suggest links or Ideas of how can i get the Map of my coutry in a big format ?

Thanks


texts: Merging two vector layers and removing overlaps in one step using PyQGIS

texts:
I have two overlapping polygon vector layers.

I want to combine these two layers into one vector layer, but I don't want to keep overlaying parts - the blue layer is of higher importance (weight), so I want to keep only blue features and if there is no blue feature, than keep the red polygon (or said in other words: fill the empty space in blue layer into the boundary of red polygon - the red layer doesn't hold important information, I use it only to know the boundary box to fill).
I know how to do it in multiple steps: use Difference tool (to clip out the blue parts from red polygon) and then Union to combine both layers, but is there please any way how to do it in only one step?
I want to implement it into Python (PyQGIS) script and I don't want to create unnecessary files that take up space and increase computing time.
processing.run("native:red", {'INPUT':blue,'OVERLAY':grid,'OUTPUT':difference})
processing.run("native:union", {'INPUT':difference,'OVERLAY':blue,'OVERLAY_FIELDS_PREFIX':'','OUTPUT':combined_layer})


texts: Output of an algorithm in QGIS Graphic Modeler

texts:
I use two Processing tools in a graphic Model in QGIS:

1-Select by Location

2-Select by Attribute

The problem is that I can't get the output of the selection by location so that I can't use it as input for other algorithms.


texts: R plot raster: changing breaks also changes colour scale

texts:
I am trying to change the breaks of a raster in R that is show with a continuous scale. However, changing just the breaks also changes the colour scale. I am not sure how I can fix this?

library(raster)
r <- raster(xmn=1, xmx=5, ymn=1, ymx=5, nrows=5, ncols=5)
r[] <- 1:length(r)

Normal plot:

plot(r)

Now just change the breaks, say to 21 instead of 20:

plot(r, breaks=c(0, 5, 10, 15, 21, 25))

Changing breaks also changed the colour!? Why? My guess is that now plot.raster() uses the 5 (or 25?) first colour of the colour ramp, which are all pretty much the same? The issues seems that it handles this now as a categorical colour palette? Indeed, giving a custom ramp equally spaced with 25 colours, it seems it only uses the first 6?

cols <- rev(terrain.colors(255))
plot(r, breaks=c(0, 5, 10, 15, 21, 25),
 col=cols[round(seq(1,255,length.out=25))])

cols <- rev(terrain.colors(255))


texts: Dynamically set mapview based on the overlay feature of L.geoJson coordinates in Leaflet

texts:
I want to dynamically setView of the map based on the overlay feature which I extract from L.geoJson in Leaflet. I tried using fitBounds, which doesn't seem to work. Here is the code snippet:

 var counties = $.ajax({
            url:"http://localhost:8080/geoserver/........../&maxFeatures=50&outputFormat=application/json",
          dataType: "json",
          success: console.log("data successfully loaded."),
          error: function (xhr) {
            alert(xhr.statusText)
          }
        })

       $.when(counties).done(function(){
        var map = L.map('map')
           .setView([52.30, 8.23], 15); //I want to set the coordinates dynamically here like setActiveArea
         // .setActiveArea('activeArea', true, true);

        // Map layer

        var Esri_WorldImagery = L.tileLayer('https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}', {
       attribution: 'Tiles &copy; Esri &mdash; Source: Esri, i-cubed, USDA, USGS, AEX, GeoEye, Getmapping, Aerogrid, IGN, IGP, UPR-EGP, and the GIS User Community'
       }).addTo(map);

      //Style layer for counties
      var geojsonMarkerOptions = {
            radius: 8,
            fillColor: "#ff7800",
            color: "#000",
            weight: 1,
            opacity: 1,
            fillOpacity: 0.8
        };

    // Add requested external GeoJSON to map
     var ky = L.geoJSON(counties.responseJSON, {
                style: geojsonMarkerOptions
      }).bindPopup(function (layer) {
                return ("Name of the field:"+(layer.feature.properties.name)+'<br/>'+"Name" +(layer.feature.properties.sid)+"sid");
      }).addTo(map);

  //Just a try, doesnt seem to work either
     ky.on('load', function() {
     map.fitBounds(ky.getBounds(), {maxZoom: 3});
        });
    });

I tried using Leaflet plugin https://github.com/Mappy/Leaflet-active-area but got a blank output. Is there any other way to dynamically use the overlay coordinates and set it as a map view?


texts:
I need to do some calculation between two Raster, and I first find the arcpy.plus_3d etc., but I find that if I can use Numpy to do this calculation, it would be more simple.

I get stuck in the first step, the only command I type in is:

**

import arcpy,numpy
arr1=arcpy.RastertoNumpyArray("C:/test.tif")

**

**
error message:

_RasterToNumPyArray(*args, **kwargs)
RuntimeError:<exception str() failed>

**

and the error is shown again and again:

since I had built the pyramids before, I think it's not its falut.

I find that all the others can use this freely, but I can not solve this problem, I use Arcgis 10.4 and the version from arcgis is 2.7.13, maybe I miss some package to install?


texts: Why is export of joined table half empty (all zeroes) even though joined table was fine?

texts:
I have joined two tables (excel 10 but saved as 03 to import it into ArcGIS) in ArcGIS 10 based on a column. This column has some empty records in both tables. I need the agreeing records to be added to the right and the others to be on the bottom of the table. After fiddling with it I got it to work; we filled in the empty cells of that column with dummy numbers. Anyway, the problem now is that when I export the joined table, the second (right) half of it is empty (all zeroes) even though the joined table was just fine.
I tried exporting it as dbf and txt and neither worked.

Does anybody have an idea how to fix this?
Thanks a bunch!


texts: How to use lastrees on a catalog with a function requiring a CHM raster?

texts:
I'm still fairly new to using LiDAR data and R, so please let me know if this isn't clear.

I am using the lidR package to segment individual trees over a large landscape. I have been using the catalog function to process the data in reasonable chunks, and have created and saved the CHMs for each chunk.

What I can't figure out is how to use the lastrees function on the catalog and call the correct CHM to use with the algorithm. Is there a way to do this within the package, maybe using catalog apply (still haven't figured that out), or am I better off attacking it some other way?

If there is a way to process the CHMs to make this easier, I'm willing to try that as well!


texts: Splitting line into segments using GRASS not working

texts:
Condensed version: I wish to divide all the streets of a city into 200-meter-long segments. I tried using v.split.length in the GRASS toolbox in QGIS 2.18.13, but it didn't work for me: instead of creating 200 m segments, it divided all streets at intersections.

Long version: I obtained a shapefile of centerlines of all streets in Oakland CA from the city's database. There were multiple entries of the same street name occurring in the database, one entry for each segment of street between two intersections. For e.g., if Adeline Street goes through 3rd to 35th Streets, there would be 32 instances of "Adeline" in the database. I first processed this shapefile using the "Dissolve" feature in QGIS to merge all same-streets into one (i.e., to create one long "Adeline" street entry). However, when I attempt to run v.split.length to divide these dissolved streets into 200 m segment (to go with previous example) instead of splitting Adeline into 200 m segments, it breaks it into 32 segments, occurring between Adeline+3rd and Adeline+35th.


texts: ArcGIS Versioning - Conflicts option greyed out / disabled

texts:
I'm relatively new to the world of Versioning and using SDE databases, and have come across a problem that I can't solve.

Summed up - the Conflicts option/button is greyed out/disabled, and I know there are conflicts in the versions.

I'm sure I am probably missing out a step, so hoping someone can educate me on what I'm doing wrong please.

Background - Using ArcGIS v10.3 and data is held within ArcSDE Database. There are two primary versions - Default and FieldEditing.  We have just started testing Collector, so each time a sync occurs, another version is created.  I am trying to reconcile, review and send changes manually, feature by feature, and post to the FieldEditing version.

My steps are:

Add data from Default database (where I have full permissions) to ArcMap
Change the version to one of the Collector versions
View version changes (at this point I can see that there is a difference between the Collector versions and the FieldEditing version)
Reconcile (Target version: DATA.FieldEditing, Define conflicts: by object, conflicts to be resolved: in favor of the Target Version)
Last week I am sure that I was then able to click the 'Conflicts' button and go through the conflicts one by one, but this is now greyed out/disabled.

I've also got an issue whereby last week I had one version which had over 600 conflicts and out of the 26 versions, today there is only one which shows that there are changes.  I have just had an error to say "License Connection Lost - Lost connection with license server for ArcGIS for Desktop Standard".


texts: QGIS crashes after script is run for the second time

texts:
QGIS Version - 3.0.1
I am trying to create a shapefile with geometry set to 'Point' and SHPT set to 'Point'. But before that I am facing following problems.

When I run the below script through QGIS, it gives me following error the first time.

    exec(open('C:/Arundhati/Automation Test Script/Radio Signal/QGIS_Test.py'.encode('utf-8')).read())
    Traceback (most recent call last):
      File "C:\PROGRA~1\QGIS3~1.0\apps\Python36\lib\code.py", line 91, in runcode
        exec(code, self.locals)
      File "<input>", line 1, in <module>
      File "<string>", line 12, in <module>
    AttributeError: 'QgsVectorLayer' object has no attribute 'layer'

When I run the same script second time, QGIS crashes.

When I run the same script through OSGEO4W shell, it runs without any errors but I never prints the messages in the print statement. So I am not sure if it is running the script properly.

C:\Program Files\QGIS 3.0\bin>python-qgis.bat "C:\Arundhati\Automation Test Script\Radio Signal\QGIS_Test.py"
C:\Program Files\QGIS 3.0\bin>

What am I doing wrong?
Following is the entire script which is causing trouble.
from qgis.core import *
import qgis.utils
import sys

QgsApplication.setPrefixPath("C:/Program Files/QGIS 3.0", True)
qgs = QgsApplication([], False) ###
qgs.initQgis() ###

_layer1 = QgsVectorLayer(
             "C:/Arundhati/Automation Test Script/Radio Signal/SIGNAL.shp",
             "SIGNAL",
             "ogr")

print("I am here" + _layer1.layer.wkbType())

_writer1 = QgsVectorFileWriter(
             "C:/Arundhati/Automation Test Script/Radio Signal/SIGNAL.shp",
             'utf-8',
             _layer1.fields(),
             Qgis.WKBPoint,
             QgsCoordinateReferenceSystem(),
             'ESRI Shapefile')

#add a feature
fet = QgsFeature()
fet.setGeometry(QgsGeometry.fromPoint(QgsPoint(10,10)))
fet.setAttributes([1,"text"])
_writer1.addFeature(fet)

del _writer1

qgs.exitQgis()


texts:
I am using Python 3.6 and still having this problem,, although I have tried to uninstall both numpy and GDAL and installed them again in the order numpy then GDAL
I am using Windows 10.
This was my code:
import numpy
import PIL
import osgeo
from osgeo import gdal
import raster2xyz
from raster2xyz import Raster2xyz
from PIL import Image

and this is the error?
Traceback (most recent call last): File "C:/Users/user/PycharmProjects/Tutorials/19august.py", line 3, in import osgeo File "C:\Users\user\AppData\Local\Programs\Python\Python36\lib\site-packages\osgeo_init_.py", line 32, in gdal = swig_import_helper() File "C:\Users\user\AppData\Local\Programs\Python\Python36\lib\site-packages\osgeo_init.py", line 31, in swig_import_helper return importlib.import_module('gdal') File "C:\Users\user\AppData\Local\Programs\Python\Python36\lib\importlib_init.py", line 126, in import_module return _bootstrap._gcd_import(name[level:], package, level) ModuleNotFoundError: No module named '_gdal'

Process finished with exit code 1


texts: Hatch pattern instead of solid fill to highlight selected polygons?

texts: Can I highlight selected polygons with a hatch pattern instead of solid fill (I know that fill transparency is an option)? Either globally or on a project basis?


texts: How to merge 48 separated US States layers into 1 layer, using MapInfo?

texts:
They are Census Block layers. The person on the website separated them into states so they can be download easier.  I uploaded each states separate layer into MapInfo.

Now I need to merge the files into one US Layer.

They all have the same table structure and fields.  I just downloaded MapBasic but I am unsure how to really use it.
I tried to laso the layers to select them all and then just save the selected. Only works on the top layer.

Please help! Thank you all for your help in advance!


texts: ArcGIS ModelBuilder and Inline Variable Substitution Error

texts:
I have a feature class consisting of about 100 or so lines. I'm trying to use ModelBuilder to generate a profile table for each line so that I can generate elevation profiles for each one. My model is below.

I want each profile table to have the same name as the Name variable so the path ends with %Name%. Whenever I try to run the model like this, it returns Error 999999: Something unexpected caused the tool to fail. Contact Esri Technical Support (http://esriurl.com/support) to Report a Bug, and refer to the error help for potential solutions or workarounds. I would rather not contact ESRI support.

I've come to the conclusion that the inline variable substitution is what's causing the problem as when I run the model giving Profile Table a fixed name at the end of the path, it runs to completion and one profile table is present in the destination folder. It's the profile table for the last feature that gets iterated. I'm not very experienced with ModelBuilder so it may be something obvious that I'm missing.

Does anybody have any idea what the root of this issue may be?


texts: Error in installing QGIS 2.6 using Osgeo in Windows 8.1

texts:
I was trying to load QGIS 2.6 on Windows 8.1 using OsGeo installer. Though the files are loading fine while launching the software it is showing a series of errors. Some on these are -

The procedure entry point xmlShemaNewDocParserCtxt could not be located in the dynamic link library C:\OSGeo4W\bin\Spatialite4.dll
Package: Unknown package
qgis-common.bat exit code -1073741515
qgis-dev.bat exit code -1073741515
Package: Unknown package
qgis-common.bat exit code -1073741511

I don't know how to get across these issues.


texts:
I have a situation where I'm getting small polygons in the data I'm receiving and I would like to merge these using Shapely so that they would be combined with the larger ones. Here is an example of what I'm talking about.
The image is generated from Shapely's GeometryCollection using just polygons.svg(). Using cascaded_union seems to merge the small polygons, but the problem is that after that all the walls are just merged and I'm left with this single polygon which is not what I'm looking for.
How could this be achieved?


texts: How to adjust QGIS attribute table view column width using Python?

texts:
I'm programmatically opening the Attribute table via Python through a plugin. All of the columns have a default width, but some of the data that I'm using do not fit.

It's not enough to be able to click the column separator to make the column re-size to fit the data as there are several columns and several layers.

Can I, through Python, change this in any way, or is my only option to create my own table in Qt?


texts: Linking points within a specified tolerance

texts: I have a set of points, and I want to build lines between pairs of points whenever their separating distance is within a specified tolerance. For example, 3 or 4 meters. Is there a way to do it in FME?


texts: How to only label features in layer 1 that intersect with buffer polygon in layer 2?

texts:
I'm using QGIS 2.82 and 2.14 (different production workflows but same problem).

I have two layers:

A layer of polygons (layer1)
A 'buffer' layer (layer2), i.e. a single circular feature centred on a point of interest

What I want is to:

Render both layers
Only label layer1 polygon features that intersect the buffer polygon

I'm trying to use the Label -> Expression string builder with:

CASE
WHEN intersects(layer1,layer2) THEN layer_1_label_column
END

I'm getting a parse error:
Column 'layer1' not found

Anyone know if this is because:

I haven't got the syntax right to point to the geometry in the layer (I've tried dbname.layer1.the_geom and layer1.the_geom with and without quotes);
It's not possible to use a geometry function across two layers in a label expression.

I don't want to do anything that results in a permanent output (i.e. additional columns added to the layer).


texts: Microsoft Excel 2016 3D Maps - better vector data?

texts:
Microsft Excel 2016 comes with a feature called 3D Maps, which charts your spreadsheet data on a vector map base. See Get Started with 3D Maps.

The data seems to be derived from Bing Maps, which in turn gets its vector data from somewhere else. The coastline vectors are crude, too crude for local mapping. Would it be possible to somehow swap in better coastline vectors?


texts: How to print selected page from Data Driven Pages using ArcPy?

texts:
I am looking to print DDP (Data Driven Pages) to a printer through a Python (arcpy) script which I will create in Arc Catalog. In other words, I would like for the script to have the following parameters:

Select the MXD,
Select the feature class,
SQL Expression (to select an attribute),
Then a Target Field (such as the name of the column in a field).

I have a small script which I wrote, but it prints the entire contents DDP in the library layer instead of letting a user select only a specific page in DDP.

How can I or a user select one page and print the selected page?
Additionally, what parameters should be set in the script?

import arcpy
from arcpy import env
#
# Set the workspace
env.workspace = arcpy.GetParameterAsText(0)
dir = env.workspace
#
#
# Local Variables
input = arcpy.GetParameterAsText(1)
expression = arcpy.GetParameterAsText(2)
#
#
# Make a layer from the input feature class
arcpy.MakeFeatureLayer_management(r'Database Connections\gis_sql3.sde\gis_parcels.MP.Districts\gis_parcels.MP.Library', "lib_lyr")
#
#
# Within selected features, further select based on a SQL query within the script tool
arcpy.SelectLayerByAttribute_management("lib_lyr", "SUBSET_SELECTION", expression)
#
#
#
mxd = arcpy.mapping.MapDocument (r"M:\CCAO_GIS_Projects\Library_Districts\MXD\Lib_Test_new.mxd")
mxd.dataDrivenPages.getPageIDFromName("LIBRARY")
# As of now this prints only the first page
mxd.dataDrivenPages.pageNameField.name = targetPageName
mxd.dataDrivenPages.printPages(r"HP Color LaserJet 2600n (Copy 1)", page_range_type = "CURRENT")
mxd.dataDrivenPages.refresh()
del mxd


texts:
What is the correct way to fetch vector features from the server using GeoJson?

{
  "type" : "FeatureCollection",
  "features" : [ {
    "type" : "Point",
    "coordinates" : {
      "x" : 14.8802,
      "y" : 22.1043,
      "z" : "NaN",
      "m" : "NaN"
    },
    "properties" : {
      "date" : 1472936467273,
      "device" : null,
      "longitude" : 14.8802,
      "latitude" : 22.1043,
      "mi_PRINX" : 1,
      "id" : 1
    }
  }]
}

The map shows no points.

  var dynamicFeatures = new ol.layer.Vector({
          source: new ol.source.Vector({
          format: new ol.format.GeoJSON(),
          url: "http://localhost:9090/map/layers"
        })
      });

      // Map
      this.mapObject = new ol.Map({
        target: 'map',
        controls: [],
        layers: [
          new ol.layer.Tile({
            source: new ol.source.OSM()
          }),
          dynamicFeatures
        ],
        view: new ol.View({
          center: [949282, 6002552],
          zoom: 4
        })
      });

URL is correct.

Could it be because the z and m values? these are generated automatically from the server...

Regards,
Ido


texts: Show map extent as a polygon onto another map - OpenLayers 3

texts:
How can one add a polygon that represents at all time the extent of map1 onto map2 as in the following page?

http://www.alexisathlani.com/paris.html


texts: Calculating shortest distance between multiple Polygon borders to a common Polygon border?

texts:
I have 2 sets of Polygons and I need to calculate the shortest distance between the polygons on one layer to the single polygon (one feature) on the other layer! How can I do that? Any help on how to accomplish this in QGIS?

The real world case is to calculate the shortest distance from the edge of a land parcel to the water/ocean. I will need to do this once or twice, but the multiple polygon layer has 10k's of features).

Being not of many words I have included an image of what I need to do!


texts:
When you edit OpenStreetMap you have access to the Maxar Premium Satelite Imagery. That is more recent than other imagery in certain areas.

How can I get Maxar Premium Satelite Imagery as a layer into QGIS?


texts: Configuring WorldImage data store starting from folder containing tif and tfw files

texts:
I have a folder of tif and tfw files and I would like to show the raster in a map using geoserver. Even reading the geoserver tutorial for this kind of datastore (http://docs.geoserver.org/stable/en/user/data/worldimage.html) i coldn't find the way to show the raster. In the same folder i have also other files (shape, aux, properties, ..) but i don't know if they are useful for publish the map or not.

I could not figure out how to configure the pyramid of tiff.

Is there a way to configure geoserver or directly OpenLayers in order to show the folder of tif?


texts: Create polygons enclosing same number of points on PostGIS

texts:
I  have a geodatabase on Postgis with points representing customers, what i need is give a number of customers, lets say 100 and it need to generate polygons enclosing 100 clients each one. I only have generated a grid but its uniform and don't take in count the number of clients.

Yes @John Barça and @underdark your solution of using cluster points and then convex hull can work bot i need to create clusters with same number of points, do u have an idea how this can be done on PostGis?? sample i have 100 points and i want bunches of 10 based on proximity.


texts: Add drawing functionality in type ui.Linker Map

texts:
is there drawing functionality in the ui.Linker Map object of Earth Engine?

var linker = ui.Map.Linker(maps)

function display3 ()   {
    var map = new ui.Map();
    map.addLayer(toFeature);
    map.style().set('cursor', 'crosshair');
    map.centerObject(prorestaura_mg, 9)
    map.addLayer(select_2016,visPar432, 'Com nuvem - 2016');
    map.addLayer(select_2016.mask(select_2016.select('BQA').eq(2720)),visMosaic, 'Sem nuvem -2016',false);
    map.addLayer(prorestaura_mg.draw({color: 'F8E9CA', pointRadius: 5}), {},'PRO-RESTAURA/MG',false);
    map.add(ui.Label('DATA:26-08-2016', {position:'bottom-center', fontWeight: 'bold'}))
    map.ui.enableDrawing(false)
    return map}    } });


texts:
I'm quite new to python and hoping for some assistance with a script.

I have a shapefile covering over 100 map sheets and would like to clip the data from this shapefile to the area of each map sheet. I have used a separate script to create a personal geodatabase named for each map sheet and each map sheet clip will be saved as a feature class in the corresponding personal geodatabase.
I know I need to query for each map sheet to perform the clip, but am unsure how to proceed from here.


texts: How to upload GeoJSON to CartoDB in a CRS that is not lat/lon?

texts:
I'm pretty sure I need to convert the coordinates in my GeoJson file to EPSG:3857 for them to play nicely with CartoDB, but I don't know what coordinate system they are currently in. The original data is from Ordnance Survey. Here's how I got that data and how I transformed it:

I downloaded the UK OS shapefile for "lower layer super output area" boundaries.
I used the mapshaper command line tool to heavily simplify the polygons and export the result as GeoJson.
(I then used a Ruby script to discard all the features that were not in Greater Manchester, but I hope that's not relevant.)
Finally, I imported the result as a CartoDB dataset but saw no polygons on the map view.

Here is an example feature from my GeoJson file:

{
"type": "Feature",
"properties": {
    "LSOA11CD": "E01004766",
    "LSOA11NM": "Bolton 005A",
    "LSOA11NMW": "Bolton 005A"
},
"geometry": {
    "type": "Polygon",
    "coordinates": [
        [
            [
                370664.009,
                411365.069
            ],
            [
                370452.531,
                411428.156
            ],
            [
                370344.061,
                411733.869
            ],
            [
                370510.978,
                411885.617
            ],
            [
                370774.826,
                411942.479
            ],
            [
                370795.626,
                411526.363
            ],
            [
                371172.854,
                411677.763
            ],
            [
                371456.953,
                411961.328
            ],
            [
                371525.9,
                411866.915
            ],
            [
                371452.615,
                411272.858
            ],
            [
                371336.156,
                411156.394
            ],
            [
                371161.151,
                411170.972
            ],
            [
                371052,
                411185.906
            ],
            [
                370664.009,
                411365.069
            ]
        ]
    ]
  }
}

If I can find out what the current coordinate system is then I am confident that I can convert using QGIS.


texts: Vector Data for Disputed Island Areas (Paracel and Spratley Islands)?

texts:
I am looking for a free vector data set for the Paracel and Spratley Island archipelagos.

Natural Earth, suggested on StackExchange, does not include these features, due to their small size.

Has anyone seen such a layer, ideally including island names?

My interest is not in their sovereignty, but rather in showing their location and potentially spatial relationship with coral reefs. As such, data could be in any language (although English is preferred). Other island groups in the region are easy to find, but these two seem to be "missing" from many data sets because they are not assigned to any one country.


texts:
I wish to add a new value in my domain and I get this message :

An invalid SQL statement was used [UPDATE GDB_Domains SET DomainName = 'Statut_Pygargue',Description = Statut de l'activite du nid', MergePolicy = 3, SplitPolicy = 3 WHEREID = 20]

I am the owner of the file geodatabase.


texts:
I would like to create centerlines for many glaciers I am doing work with in ArcMap. I feel like I have a conceptual model that could work but can't seem to figure out how to make it happen with the tools available. I have polygons for each glacier outline, and DEMs for elevation data.

My idea is as follows:

Using zonal stats find max/min elevation for each glacier using the outline polygons
Turn the max and min pixel into a point file
Then find either a least cost path using the DEM between the max/min for each glacier(seems complicated) or find the path between the two points that lies in the center of either side of the polygon.

It's the last step I need help with.


texts:
When I export to webmap with qgis2web, it creates a new folder with a default name (for example "qgis2web_2018_11_14-17_24_47_123987") that relays on the precise time when it is created, but I need to know the exact name of the folder where the webmap it's going to be before being exported.

Is there a way to change these default folder names?


texts: Improving Maplex placement when ArcMap callouts cluttered

texts:
I use ArcMap 10.2 to produce map document containing a lot of label boxes with line callout. Sometimes, I have so many labels that leaders come very matted. Then I use Illustrator to unravel.

Is it possible to tune Maplex to prevent leaders from superposing labels?

Unfortunately I can't reduce the size of boxes (graphic charter) neither its content. And I need the data exhaustivity.

Here are my Maplex parameters:

View unplaced labels
Placement quality to best
The red, gray and black lines label weight ranking to 1000
Leader tolerance = 0
Regular placement
Label offset = 15 mm with secondary offset between 15 mm and 200 mm
Label buffer to 0%
Background label and 'Never remove' options unchecked


texts: Calculating population within a circle excluding intersects with shapely

texts:
I have a data of a centroid and a count:
    centroid               population id
0   POINT(38.9081 77.0349)  9004.0     0
1   POINT(38.8521,77.0377)  517.0      1
2   POINT(38.9081 77.0377)  10369.0    2
...

The count is a count of the population within a 15km radius from the centroid.
I need to:

take each centroid, buffer it with a 15000 m radius and find out which other circles (radiuses around centroids) in the dataset it intersects.
If it does intersect with one, or more, then calculate the total count of the population within the intersects(one or more) and then subtract those population counts from the count of the circle. In other words, for each buffered centroid, I need to get a shape of the circles without the intersects and count the population within that shape.

Here is my code:
import geopandas as gpd
import geopy
import geopy.distance
from shapely import wkt
from shapely.geometry import Point, LineString, Polygon
import pyproj

df['centroid'] = df['centroid'].apply(wkt.loads)

for i, r in df.iterrows():
a = gpd.GeoSeries(r.centroid, crs = 4326).to_crs(crs = 3857).buffer(15000)
df_dropped = df.drop(i, axis = 0)
l1 = [] # an empty list to collect intersections
for j, s in df_dropped.iterrows():
    b = gpd.GeoSeries(s.centroid, crs = 4326).to_crs(crs = 3857).buffer(15000)
    tf = a.intersects(b)
    l1.append(tf)
    print(l1)

My idea is to get the list populated with intersections with the other circles (True) and then use their indexes to subtract those buffered centroids from the initial buffered centroid to get the semi-moon for the original buffered centroid. If there are more than 1 intersect, then I would subtract it from the previous result, etc., until there are no more intersects. I would then calculate the population for that shape.
What would be the code to achieve that?


texts: CRS Error when clipping NetCDF using a shapefile

texts:
When trying to mask a NetCDF (a daily time series of 31 days) with the code below, I get the following error related to the part of writing the CRS to the rio NetCDF data:
CRSError: The WKT could not be parsed. OGR Error code 6

Here is the code I am using (taken from snowman2's answer here: How to mask NetCDF time series data from a shapefile in Python?):
import rioxarray
import xarray
from shapely.geometry import mapping

ERA5_daily = xarray.open_dataarray('ERA5_VIMD_2018_08_Daily.nc')
ERA5_daily.rio.set_spatial_dims(x_dim="longitude", y_dim="latitude", inplace=True)
ERA5_daily.rio.write_crs("epsg:4326", inplace=True)
shapefile = geopandas.read_file('Cluster_34586.shp', crs="epsg:4326")

clipped = ERA5_daily.rio.clip(shapefile.geometry.apply(mapping), shapefile.crs, drop=False)

Does anyone know how to fix it?


texts: Adding Here WeGo Satellite imagery to QGIS?

texts:
Here WeGo has more updated imagery for Cambodia than Google.  I would like to add the WeGo image tile to QGIS.  I have followed the instructions here and can add the WeGo street layers to my map of Cambodia.

I have followed this link to add satellite imagery.  I took the original XML file posted in the link and altered the URL but used the same key as the street layer version of the XML.
Here is a screen grab of how/where I got the app ID and app key.

This resulted in a layer that repeats a tile of southern Spain.  I am unaware of how to alter my XML code to get the WeGo image tile to draw for southern Cambodia.

How can I alter the code to get WeGo imagery to draw in QGIS?

<GDAL_WMS>
<Service name="TMS">
<ServerUrl>https://2.aerial.maps.api.here.com/maptile/2.1/maptile/newest/satellite.day/5/15/12/256/png8?app_id=tEu0WBZmhFgA0aSvAJfM&amp;app_code=V32dyBBowZHxPHGMUTQH2g</ServerUrl>
</Service>
<DataWindow>
    <UpperLeftX>-20037508.34</UpperLeftX>
    <UpperLeftY>20037508.34</UpperLeftY>
    <LowerRightX>20037508.34</LowerRightX>
    <LowerRightY>-20037508.34</LowerRightY>

<TileLevel>20</TileLevel>
<TileCountX>1</TileCountX>
<TileCountY>1</TileCountY>
<YOrigin>top</YOrigin>
</DataWindow>
<Projection>EPSG:3857</Projection>
<BlockSizeX>256</BlockSizeX>
<BlockSizeY>256</BlockSizeY>
<BandsCount>3</BandsCount>
<Cache />


texts: Using memory management and error handling to improve performance of ArcPy script?

texts:
I have an ArcPy script that adds new fields and then deletes unnecessary fields from a group of feature classes. I am a novice at ArcGIS and even more so at ArcPy, but have been able to piece together the code below from other working examples. I am currently using ArcGIS 10.0, and expect to get version 10.3 before long.

    import arcpy, os

try:
    # Set the workspace to edit
    arcpy.env.workspace = r'E:\Python Testing\NavajoAUM_SiteScreens.gdb\Central'
    # For loop to go through feature classes in the workspace
    for fc in arcpy.ListFeatureClasses():
    # add new fields
        arcpy.AddField_management(fc, "GammaMeas", "DOUBLE", field_precision = 9, field_scale = 0)
        arcpy.AddField_management(fc, "BkgLevel", "DOUBLE", field_precision = 9, field_scale = 0)
        arcpy.AddField_management(fc, "MineID", "TEXT", field_length = 10)
        arcpy.AddField_management(fc, "OrigFileName", "TEXT", field_length = 50)

    # Set the cursor
        cursor = arcpy.UpdateCursor(fc)
        row = cursor.next()
    # Save values for the MineID and OrigFileName
        desc = arcpy.Describe(fc)
        mineValue = desc.name[desc.name.rfind("_") + 1:]
        fileValue = desc.name
    # While loop updates values in each row
        while row:
            row.setValue("GammaMeas", row.getValue("Text"))
            row.setValue("BkgLevel", "12500")
            row.setValue("MineID", mineValue)
            row.setValue("OrigFileName", fileValue)
            cursor.updateRow(row)
            row = cursor.next()
except:
    pass

# Set the workspace to edit
arcpy.env.workspace = r'E:\Python Testing\NavajoAUM_SiteScreens.gdb\Central'
# Local Variables to keep
keep = ['OBJECTID','FID','Shape', 'Text', 'Channel', 'Northing', 'Easting', 'GammaMeas', 'BkgLevel', 'MineID', 'OrigFileName']  # be SURE to include FID/OBJECTID field here!
for fc in arcpy.ListFeatureClasses():
    try:
        discard = []
        for field in [f.name for f in arcpy.ListFields(fc)if f.type <> 'Geometry']:
            if field not in keep:
                discard.append(field)
        arcpy.DeleteField_management(fc, discard)

    except:
        pass

It seems like using memory management would greatly improve the performance of the script, and that error handling could be improved, but I'm at a loss on how to implement either of those correctly here. There might be other obvious improvements to make that I don't know about, or haven't attempted yet (e.g., maybe adding fields and deleting fields in one loop would be better).

As far as errors go, it appears that sometimes the script runs fine, sometimes it will throw an error and stop, and other times no error is raised but the script seems to stop early (e.g., not all fields added or filled, no fields deleted). I'd like to know where / which feature classes are causing the errors so I can QC the feature class / attributes table and fix.


texts: Merging multiple LAS files using Python?

texts: Is it possible to merge multiple LAS files into one using Python? Currently I'm parsing LAS data from a website. So the way I'm currently getting this to work is by, getting a bounding box around a certain address which you can input. I give this coordinates to the website where the las data is located. I then check which las files are located at that location and download the different tiles available. I do this by checking the website via beautiful soup and parse the download link to a request. It's not ascii data and i'm not reading the las files in Python for visualization. But my question is : When I get the different las files that are representing the tile from the specific location, is there a possible way to merge it ? I would love to merge these into one big LAS file that I can open in CloudCompare. Is there a library I need to know about. Or should I look at LasTools? If so how do  implement this into my Python script.


texts:
I'm unable to install FUSION in QGIS 3.2 on Ubuntu 16.04.

Andre Silva on Creating DEM from LAS file without using LAStools? says to use the same procedure as LAStools (up to QGIS 2.18 version). I've previously installed LAStools which work great. But I don't get how to do it for Fusion.

I've installed the Fusion .exe with Wine correctly. Then, as said in how to install LAStools toolbox in QGIS, unzip the FUSION plugin folder and paste it in QGIS plugin folder even if it is said to be only used for QGIS 2.2, but there is none for 3.2 (it is a 2013 post).

Nothing appears in QGIS toolbox bar or advanced configuration. Maybe it is because I've already installed LAStools. Any clues?


texts:
I've been using QGIS 2.6 (Brighton) happily for awhile but suddenly it became extremely slow during certain operations, after being reasonably fast in past use.  I'm merging NAIP JP2 files to virtual raster, reprojecting the raster to another vrt, then clipping that raster and using gdal_retile to cut it into tiles.  All files are on an external USB hard drive.

It takes more than 24 hours for the Mac (2.9ghz Core i5, 8gb ram, Yosemite 10.10.1) to cut a single tile.  Meanwhile, duplicating the exact same steps, a quad core i5 Windows machine outputs a tile every 29 minutes.  I certainly didn't expect the Mac to be faster - the Windows machine is speedier, but it is slower by a factor of 10 or more.

What is very strange is that I used to be able to cut a single tile using a vector shape as a mask layer in about an hour.  Now, that also takes over a day.  I've tried deleting and re-installing QGIS and the KyngChaos frameworks but the problem persists.

Any suggestions about what may be causing this?


texts: Calculate Field tool in ModelBuilder only returns first value of string

texts:
I am using ArcGIS Pro 2.6.1.
In ModelBuilder, I am using a Calculate Field tool to determine a compliance status based on a color value derived from a Bentley DGN file. Here is the expression used:
Compliant(!Color!)

and here is the Code Block:
def Compliant(Color):
    if Color == 2:
        x = "Compliant"
    else:
        x = "Non Compliant"
    return x

The code block correctly determines the status, however instead of writing "Compliant" or "Non Compliant" into the table, only a "C" or an "N" is captured.
This also works perfectly outside of the ModelBuilder, but behaves in the way described within ModelBuilder.
Edit - Add Field Parameters
As requested by @Hornbydd, here is an image displaying the Add Field tool parameters within ModelBuilder:

Edit - Screenshot of Attribute Table
As requested by @Hornbydd, here is an image displaying Attribute Table after executing the Model:


texts:
Ultimate goal is to pull data from a Feature Class, group out the different types of Lease Status (field), run a dissolve and sum calculation on each new layer/FC then merge all new layers/FC's together into one FC.

This is my manual process for what I am trying to accomplish: I went about creating a table that includes one Tract_ID, the lease_status and a sum calculated on the Percent_interest field for all the owners. (probably very rudimentary LOL)

(create initial  tables to work from)

Export DBLeaseOR8312m into a new layer/FC for each lease_status  for example DBLeaseOR8312m_open, DBLeaseOR8312m_ThirdParty, etc.
From the DBLeaseOR8312m_open FC select all other lease_statuses and delete
Run dissolve on the new layer

Select Tract_ID as the dissolve field
Select percent_interest and “SUM” in the statistics area

(modify new tables to include new field)

This creates a new table that only has one Tract_ID, one type of Lease_status and the percent_interest calculated for all owners and included in percent interest field.
Creat new field in this table for Lease_status and add the specific type to that field
Do this for each Lease_status type

(merge data together into one final table)

Then merge all tables into one table to include all lease_status types.

I'm thinking I can use Iterate Feature Selection in model builder but not quite sure how to bring it all together and set my outputs...


texts: Tools for discovering spatial data defects

texts:
Spatial data often have defects that may include incomplete polygons, sleathers (non-uniform administrative boundaries), overlaps, having points that belong to more than one polygon, etc. I have a collection of spatial data in .shp format, which I wish to clean up. They include administrative boundaries and population/point data.

To clean the data, I'll have to first discover the defects. This process is hectic. I'm hoping that there are tools that I can use to automate this process.

I'm using GeoServer and PostGIS. I'm looking for tools that can I can use for the defect discovery.


texts: Composite bands in ArcGIS Pro and Raster Stack in R produce different multiband raster

texts:
I have 8 rasters (GeoTIFF; LANDFIRE layers) that I need to stack into a multiband raster for use with downstream analyses (i.e., within FlamMap). I have a workflow using Python and ArcPy and would like to continue to do so for this task. I need to repeat that task for many locations, hence the desire to find a more automated solution.
All of the TIFFs have the 'correct' values before I stack them with arcpy.CompositeBands_management, but after compositing some of the nodata values change for the bands. But this only matters for one of the bands. This band is the only one that has categorical data, which it reports as text and a numerical code. This band reports nodata as 255, but in order for the downstream analyses to work this needs to be remain nodata or be reported as 99.
When I use the stack function from the raster package, the resultant stacked raster works and nodata remains nodata.
I'd prefer to keep everything within Python and ArcPy if possible. I've tried removing the other attribute information from this categorical TIFF, leaving just count and value. This did not work. Any ideas on why compositing changes nodata values or how to fix it?


texts:
I am having troubles with the encoding of the osm data in R. Here is a reproducible example using the osmar package:

osmData <- osmar::get_osm(osmar::center_bbox(23.334360, 42.693180, 100, 100))
osmData$nodes$tags[80:100, ] #the output is not UTF-8

I have also downloaded a planet file from https://download.geofabrik.de/europe/. The filename is "bulgaria-latest.osm.bz2".

After unzipping it, installing osmosis, and reading the file with osmar::get_osm(source = osmsource_osmosis(fileName)) I still have the same issue. The cyrilics letters are not readable.

I am using a 64bit Windows 10 with R version 3.3.3. My default encoding is Bulgarian, but I have also tried English encoding in R with no help.

Any ideas how can I counter this?


texts: Custom info window code to display youtube videos no longer working

texts:
We have some info window code that displays youtube videos that just stopped working. We haven't made any changes and the videos in the info windows had been working great.

Anyway, we used the following basic instructions to get the videos in the info window:

You can embed YouTube or any other web content that can be published with an IFRAME directly in your infowindows. Just like embedding your images, you can embed your videos. For example, if you have a column containing YouTube video URLs called youtube_embed_link, the HTML would look like this,

Our youtube_embed_link looks like this:

http://www.youtube.com/embed/UJDlDN9LGKU

Any reason this would stop working?


texts: How do I attach a click event to a geoJSON that then executes an Ajax function when clicked. I looked into onEachFeature but that executes when the geoJSON is loaded, not when clicked, thus executing a ton of ajax calls!


texts: Can I style a feature point based on its id?

texts:
I wonder if its possible to style a feature point which belongs to an OpenLayers vector layer with points by using its id.

For example, I have in this layer three points with ids: a,b and c.
Right now they have all of them fillColor:red.

Is it possible to grab the feature point with the id:b and change its fillcolor to green while the other points hold the previous style?

Is there an example for something like this?


texts: Convert a Toolbox (.TBX) to a Python script (.PY)?

texts:
I have a .tbx tool box, which runs  a script that is incredibly slow.
So what I want to do is to open the python script that is running below the .tbx and make some changes or perhaps write something similar and faster. Is there someway to convert a .tbx ( and not just the form ) into a .py so I can inspect the code that it runs?
to clarify,

it is a third party tool created as a .tbx, and specifically here:
http://www.arcgis.com/home/item.html?id=c1c649ad439442189547c12685fcc9c3

it is slow in running in arcmap, but that is not the central idea, here. I want to open the python underneath the .tbx. As suggested by others i have tried
to right click and export (grayed out) and
also i have tried using a .tbx to .pyt converter (only returns the input form) and one more thing
i has also right clicked and opened in wordpad, but it was nearly illegible.


texts: MicroDEM changes the aspect ratio of DEM

texts: I downloaded a DEM from Open Topography (SRTM) and opened it in MicroDEM which changed the aspect ratio. The original DEM was 2329 px * 1679 px but that in MicroDEM when it was saved on disk was 2329 px * 2134 px. No clipping happened and it just changed the aspect ratio. Why is that?


texts: How to intersect lines and polygons in R?

texts:
Following on from an answer about intersecting polygons with lines to chop up a polygon into smaller polygon units (in QGIS), I wanted to try the same thing in R. However, I cannot seem to find a method that works!

over() doesn't have a method for polygons intersecting with lines; I found gIntersection() from rgeos but it fails:

require(sp)
require(rgeos)

poly <- readShapePoly("polygon.shp")
lines <- readShapeLines("lines.shp")

chopped <- gIntersection(poly, lines)

Giving:

Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, "rgeos_intersection") :
  UnsupportedOperationException: GeometryGraph::add(Geometry *): unknown geometry type: N4geos4geom18GeometryCollectionE

Update: Here's a link to the files in question.

Update 2: PaulG notes that it works and after updating rgeos and R I got rid of the error above. Thanks PaulG ...

However, gIntersection results in a SpatialLines object no matter whether I put in (poly, lines) or (lines, poly) - whereas the operation I did in QGIS (or Arc, back in the bad ol' days) will divide a polygon with the lines and result in a polygon object, not line.

So, how do I chop up my polygon with the lines and get polygons out?


texts: Exporting map to jpg (multiple times with python)?

texts:
I set symbology and exported each raster in my mxd (one by one) to jpg file. My script is ok and it works (it is shown at the bottom).

The question is: how can I customly set my legend?

Now I have those standard labels ("High: #max.number" and "Low: #min.number"), see the image link.

Is it possible to use a python script that sets labels different, for example "Max. value: #maximum number" and "Min. value: #minimum number"?

I have 200 images and renaming each legend is time consuming...

image link (example):

MY SCRIPT:
mxd=arcpy.mapping.MapDocument("CURRENT")
for lyr in arcpy.mapping.ListLayers(mxd):
    arcpy.ApplySymbologyFromLayer_management(lyr, "my_first_layer")
for lyr in arcpy.mapping.ListLayers(mxd):
    #turn off all layers
    lyr.visible = False
for lyr in arcpy.mapping.ListLayers(mxd):
    #turn one layer on at a time and export ...
    lyr.visible = True<BR>
    arcpy.mapping.ExportToJPEG(mxd, "D:\\Users\\JPG_FILES\\"+ lyr.name
                                     +".jpg",resolution=300)
    lyr.visible = False


texts: Adding featureSet into a featureLayer (dynamic renderer)

texts:
I'm trying to use a renderer and a legend on data from a featureSet, it displayes symbols and popups but the renderer fails.

The data comes from  a featureLayer but I display it with view.graphics so I think that may be the problem. Still I can't find an add graphic (maybe addAttachment but I could not make it work) .

This is my code:

var layer = new FeatureLayer({      url:'https://services.arcgis.com/V6ZHFr6zdgNZuVG0/arcgis/rest/services/la_county_labor_centroid/FeatureServer/0'});

  function succ(featureSet) {
  layer.applyEdits({addFeatures:featureSet})
    map.add(layer)
  };

  function fail(error) {
    console.error('An error ocurred in the query: ', error);
  };

  on(dom.byId('population'), 'change', function (e) {
    var population = e.target.value;
    if (population.length > 0) {

      var query = new Query();
      query.returnGeometry = true;
      query.outFields = ["*"];
      query.where = "TOTAL_POP > " + population;

      layer.queryFeatures(query).then(succ,fail)
    }
  });

  });

heres my code pen:

https://codepen.io/segev-salman/pen/MdmrLB

the first time i use the select it refreshes the layer and adds the feature set but it doesnt query it(returns the full layer instead of asked data), but thats it, the promise returns the FeatureSet but it doesn't apply the edits
Do I need a new featureLayer or can I use the one with the URL to the service to be displayed with the renderer with queryed data and if not do I really have to insert the data manually after the query does it just bring me back graphics or changes the layer?

I know about  definitionExpression this is an example to help me learn how to edit and add data to layers.


texts: Where to find Advanced Python Field Calculator qgis syntax documentation

texts:
I am trying to follow some tutorials using Advanced Python field calculator, but when I enter the formula value = $geom.area() / < area> I get syntax error as in picture, while if I substitute the field name in pointy brackets with a number (eg. 1000) it works fine.

I have already checked the link http://hub.qgis.org/issues/9867 but still cannot find any documentation on the tool syntax


texts: Why are representations functions accessible in Table of Contents not in Catalog window of ArcMap?

texts:
It's really frustrating that adding a new representation form the Catalog window of ArcMap  does not work while it works for TOC when adding the required lyr file as a symbology then right click and "convert symbology to representation", why it's not working from the Catalog window?

However, from the Catalog window:


texts: QGIS problem with PostgreSQL/PostGIS column permissions

texts:
If a create this table:

CREATE TABLE streetmap.ways (
    pk bigserial not null,
    name character varying(150),
    geom geometry(LineString,32628),
    rowdate timestamp with time zone NOT NULL DEFAULT now(),
    username character varying(50) NOT NULL DEFAULT "current_user"(),
    CONSTRAINT ways_pkey PRIMARY KEY (pk)) WITH (
    OIDS = FALSE ) TABLESPACE pg_default;

And assign all permissions (GRANT ALL ON TABLE streetmap.ways TO qgis_editor) to qgis_editor user, everything works fine with QGIS 3.

But if instead I set:

GRANT SELECT(pk, rowdate, username) ON streetmap.ways TO qgis_editor;
GRANT ALL(name, geom) ON streetmap.ways TO qgis_editor;

Or:

GRANT SELECT ON TABLE streetmap.ways TO qgis_editor;
GRANT ALL(name, geom) ON streetmap.ways TO qgis_editor;

QGIS 3 can´t edit the table. Any suggestions?


texts: Highlighting features in qgis2leaf export

texts:
I have exported a QGIS map to Leaflet using the qgis2leaf plugin. I am looking for some assistance as to how to highlight features (polylines) when the mouse hovers over it.  I just want the polyline to stand out (become bold); i only want popups to appear when the feature is clicked (which I have already customized).
I have seen plenty of examples such as the Leaflet Cloropleth tutorial but I have two geojson files in my map that users can click between (so i'm guessing that any 'highlighting' code would need to be entered twice.

The two geojson files are exp_AllRoutesAM and exp_AllRoutesPM.  Here is the leaflet code that shows how the exp_AllRoutesPM file is handled.

 function pop_AllRoutesPM(feature, layer) {
            var popupContent = 'Survey: ' + feature.properties.Survey + '' + '
' + 'Section: ' + feature.properties.Section + '
' + 'Speed In 2013: ' + feature.properties.OUT_2013 + ' km/h' +'
' + 'Speed In 2012: ' + feature.properties.OUT_2012 + ' km/h' +'
' + 'Speed In 2011: ' + feature.properties.OUT_2011 + ' km/h' + '
' + 'Speed In 2010: ' + feature.properties.OUT_2010 + ' km/h' + '
' + 'Speed In 2009: ' + feature.properties.OUT_2009 + ' km/h' ;
            layer.bindPopup(popupContent);
        }

    var exp_AllRoutesPMJSON = new L.geoJson(exp_AllRoutesPM,{
            onEachFeature: pop_AllRoutesPM,
            style: function (feature) {
                return {weight: feature.properties.radius_qgis2leaf,
                        color: feature.properties.color_qgis2leaf,
                        opacity: feature.properties.transp_qgis2leaf,
                        fillOpacity: feature.properties.transp_fill_qgis2leaf};
                }
            });
                feature_group.addLayer(exp_AllRoutesPMJSON);
                //add comment sign to hide this layer on the map in the initial view.
                exp_AllRoutesPMJSON.addTo(map);


texts:
I am using mapbox.js 1.6.1 to serve .mbtiles created in tilemill (via mbtiles-server.php) on a private host.

my knowledge of html and js is very (very) basic, I can't figure out how to use L.mapbox.legendControl properly to add a legend to a L.mapbox.map map.

I am using mbtiles-server.php to serve .mbtiles created in tilemill. Here is what works:

//create map, centered on papua
  var map = L.mapbox.map('map', 'server/mal.tilejson').setView([-6, 140], 6);

this creates a map with the tile- and gridlayer (hover tooltips) as well as the legend.

I have not been successful with creating an empty map and adding layers with their legends showing:

//create map, centered on papua
  var map = L.mapbox.map('map').setView([-6, 140], 6);

//baselayer
  var srtm = L.mapbox.tileLayer('server/neuguinea_zl3-14.tilejson', {
                        maxZoom: 14,
                }).addTo(map);

//overlay
 var malgroup = new L.LayerGroup();
 var malariatile = L.mapbox.tileLayer('server/mal.tilejson').addTo(malgroup);
 var malariagrid = L.mapbox.gridLayer('server/mal.tilejson').addTo(malgroup);
 map.addControl(L.mapbox.gridControl(malariagrid));

//create baselayers and overlays
var baseLayers = {
  "SRTM v3": srtm
  };
var overlay = {
  "Malaria": malgroup
  };

//add layer control
 L.control.layers(baseLayers, overlay).addTo(map);

// create legend
  map.addControl(L.mapbox.legendControl());

I was hoping that map.addControl(L.mapbox.legendControl()); would activate the overlay's legend automatically, but it doesn't.

Thanks for looking into this!


texts: Selecting POIs around specific buildings using PostGIS

texts:
The problem: I need to select,  for each residential building in my table that has say at least 2 pharmacies and 2 education centers within a radius of 1km, all POIs (pharmacies, comercial centres, medical centers, education centers, police stations, fire stations) which are within 1km of the respective building.
table structure->

building (id serial, name varchar  )

poi_category(id serial, cname varchar) --cname being the category name of course

poi(id serial, name varchar, c_id integer)-- c_id is the FK referencing poi_category(id)

all coordinate columns are of type geometry not geography (let's call them geom)

here's the way i thought it should be done but i'm not sure it's even correct let alone  the optimal solution to this problem

    SELECT r.id_b, r.id_p
    FROM (
         SELECT b.id AS id_b, p.id AS id_p, pc.id AS id_pc,pc.cname
         FROM building AS b, poi AS p, poi_category AS pc
         WHERE ST_DWithin(b.geom,p.geom, 1000) AND p.c_id=pc.id
         ) AS r,
         (
         SELECT * FROM r GROUP BY id_b
         ) AS r1

     HAVING  count (
                       SELECT *
                       FROM r, r1
                       WHERE r1.id_b=r.id_b AND r.id_pc='pharmacy'

                    )>1
                 AND
                 count (
                       SELECT *
                       FROM r, r1
                       WHERE r1.id_b=r.id_b AND r.id_pc='ed. centre'

                    )>1

Is this the way to go for what i need ? What solution would be better from a performance point of view?  What about the most elegant solution?


texts: Examples where customized vector Google Maps rendered via WebGL are used?

texts:
Some of the new features are exciting in Google Maps. Creating custom map layers in the cloud and specifically rendering vector maps in WebGL.
https://cloud.google.com/blog/products/maps-platform/introducing-new-maps-customization-features
Do you have real-world examples where I can see this in action?


texts: Create new shapes by merging an existing shape

texts:
I want to create a new shapefile from an existing one.

I would like to create a new shapefile for next level of geo hierarchy (say regions) by grouping some of features (states) from the existing shapefile.

I have a CSV file which contains the geo hierarchy relationship (i.e., which all states to be combined to form a region).

Imported the CSV file as a text delimited layer
Joined it with the existing shapefile
Using "Select Features using Expression" I am able to select all the sales which forms the particular region
Then I went to Edit-->Merge Selected features to merge the selected featues (states) to form the region

Doing so I was able to create one region but since I want to create 100 regions by joining several states, is there not an automated way like using some SQL, etc., wherein I can create a group by query to group various states (features)? And also one region could overlap with another i.e. one state can belong to two-three regions, is it possible to have a overlapping shares layered over each other?


texts:
I have an outlet (point inside a polygon, see the image below) and using global mapper I have computed the basin leading to this outlet.

I want also to get the main water channel of this basin in order to calculate the time of concentration, but I could not figure out how to do that. Below is an image of the data:

I want to use QGIS.


texts: DBSCAN and CONVEXHULL: how to obtain coordinates of the convexhull polygons

texts:
Iam using DBSCAN to cluster coordinates together and then using convexhull to draw 'polygons' around each cluster. I then want to extract the coordinates of the convex hull objects however the code im using (hull.simplices) does not match the coordinate system im using.

here is the full code:

import pandas as pd, numpy as np, matplotlib.pyplot as plt
from sklearn.cluster import DBSCAN
from scipy.spatial import ConvexHull

Lat=[10,10,20,23,27,28,29,34,11,34,66,22]
Lon=[39,40,23,21,11,29,66,33,55,22,11,55]

D=list(zip(Lat, Lon))
df = pd.DataFrame(D,columns=['LAT','LON'])

X=np.array(df[['LAT', 'LON']])

kms_per_radian = 6371.0088
epsilon = 1500 / kms_per_radian
db = DBSCAN(eps=epsilon, min_samples=3)

model=db.fit(np.radians(X))
cluster_labels = db.labels_

num_clusters = len(set(cluster_labels))

cluster_labels = cluster_labels.astype(float)
cluster_labels[cluster_labels == -1] = np.nan

labels = pd.DataFrame(db.labels_,columns=['CLUSTER_LABEL'])

dfnew=pd.concat([df,labels],axis=1,sort=False)

z=[] #HULL simplices coordinates will be appended here

for i in range (0,num_clusters-1):
    dfq=dfnew[dfnew['CLUSTER_LABEL']==i]
    Y = np.array(dfq[['LAT', 'LON']])
    hull = ConvexHull(Y)
    plt.plot(Y[:, 1],Y[:, 0],  'o')
    z.append(hull.simplices)
    for simplex in hull.simplices:
        ploted=plt.plot( Y[simplex, 1], Y[simplex, 0],'k-',c='m')

plt.show()

print(z)

Please help to advise the code required to get list of coordinates for each convexhull object(to be same as original coordinates ass see in the plot)

This is output iam getting:

[array([[2, 0],
       [5, 4],
       [5, 2],
       [3, 0],
       [3, 4]], dtype=int32), array([[0, 1],
       [2, 1],
       [2, 0]], dtype=int32)]

As can be seen, it doesnt match the coordinates i used.

Appreciate your help on this. thanks


texts: Installing Qgis 2.16 on Windows 10 - MSVCP120.dll?

texts:
For some problems in my PC it was need to restore Windows 10 and when I tried to launch Qgis (2.16 and 2.14) a message from Windows appeared (see image in attached; sorry, the message is in italian):

So I decided to Unistall Osgeo directory and install it again using OSGeo4w install. Everything worked fine untill it was installing qgis-common. At that moment the previous Windows message appeared again. After OK, installation continued and at the end appeared:

  OSGEO4W message - Postinstall script error "Package: qgis-common,
  qgis-common.bat code-1073741521".

So I decided to install Qgis 2.14 but I had always the same problems.
Somebody can help me? .


texts: Convert points inserted to line and stored in a column in sql server in a trigeer

texts:
I'm basically making a script convert points inserted to line and stored in a column in sql server in a trigeer but I don't find it.
I tried:

DECLARE @BuildString NVARCHAR(MAX)
SELECT @BuildString = COALESCE(@BuildString + ',', '') + CAST([Lng] AS NVARCHAR(50)) + ' ' + CAST([Lat] AS NVARCHAR(50))
FROM dbo.ligneTemp

SET @BuildString = 'LINESTRING(' + @BuildString + ')';
DECLARE @LineFromPoints geometry =geometry::STLineFromText(@BuildString, 4326);
SELECT @LineFromPoints
END