texts: Using Copy Traversed Features to count turns for a route result in Network Analyst

texts:
I have the output of the route layer. And I want to determine the number of right and left turns for each route. For that, I used the COPY TRAVERSED ROUTE FEATURES tool from Network Analyst. However, this gave an empty output table for the turns. Are there any prerequisites for the network dataset to use the COPY TRAVERSED SOURCE FEATURE tool?
Or is this possible by using Global Turns in the analysis by any way?

I am looking for the result similar to below (from a study of Dalumpines and Scott).


texts: I am having trouble finding resources for night lights layers for use in QGis. It seems these are available for ArcGis (Earth at Night Package), but unfortunately I don't have access to this and can't find a way of transferring ESRI packages into a form usable by QGis. Are there any packages or tutorials about for loading night light images into QGis?


texts: ST_LineSubString Query Not Giving Expected Results

texts:
I have a multilinestring that is approximately 9700 meters long. My goal is to merge the multilinestring into a single 9700m linestring, then split the merged linestring into 97 linestrings of 100m lengths. When using the query below (taken from the docs for ST_LineSubstring), I get the expected results, but I cannot view the results in QGIS. I suspect I am not understanding something conceptual about the coordinate system?

The SRID for public.subdivide_lineTesting is 4326 and I'm using PostGIS 2.2.1.

SELECT ST_LineSubstring(the_geom, 100.00*n/length,
         CASE
           WHEN 100.00*(n+1) < length
           THEN 100.00*(n)/ length
           ELSE 1
         END) AS the_geom
FROM (
  SELECT ST_LineMerge(geom) AS the_geom,
         ST_Length(geom::geography) AS length
  FROM public."subdivide_lineTesting"
) AS t
CROSS JOIN generate_series(0,10000) AS n
WHERE n*100.00/length < 1

Also, a side observation:
When I count the number of points in the geometry produced using the below query, every tuple is null except for the last one. That record shows an expected 2 points.

st_numpoints(ST_LineSubstring(the_geom, 100.00*n/length,
         CASE
           WHEN 100.00*(n+1) < length
           THEN 100.00*(n)/ length
           ELSE 1
         END)) as points


texts: Import Shapefile to PostGIS in QGIS' DB Manager without using its ID column

texts:
I have a shapefile with a "id" field in which some non-unique ID numbers exist (don't ask).

When I try to import this to a PostGIS database using QGIS' DB Manager, it tries to re-use those IDs are unique primary key and fails.

The "Primary key" option is not ticked. "id" is shown in grey next to it.

How can I make DB Manager ignore an existing "id" column and just create a new SERIAL field or something like that?

QGIS 2.18.


texts: What is an efficient Algorithm to draw a straight line from points?

texts:
I have some cross section survey data, which is like this:

X Y Z pt_dt
2970969.635 359725.0088 83.4242 1-x
2970968.278 359722.2182 83.2591 1-x
2970941.771 359670.127  83.0655 1-x
2970961.369 359708.6424 83.4785 1-x
2970950.888 359688.048  85.2994 1-x
2970955.084 359696.2922 83.6479 1-x
2970953.47  359693.1165 84.9628 1-x
2970948.813 359683.977  83.1451 1-x
2970931.904 359650.7413 83.1042 1-x

When I use the X & Y to make points, I can see the points on the map like this:

What I need to do now, is to get the points in linear order (i.e. from left to right, or right to left), so that I can create a 3D line, which can serve as the Cross sectional profile.

My current algorithm, is to calculate the distance between every pair of points(For one given cross section, which is indicated by pt_dt). The pair which has the maximum distance, are the outermost points. I then arrange the points based on distance from one of the outermost points. I have implemented this algorithm in python, and it works.

This algorithm is O(n2); My entire data consists of about 5000 such surveys, so any complexity that I can reduce will be beneficial.

So Is there a more efficient way, or algorithm to do this?


texts: Connect signal to function with arguments (Python)

texts:
I have this function in Python that take 2 arguments:

     def niveau(controName,idNiveau):
              ........

I want to connect it to this signal:
QObject.connect(dialog.findChild(QDialogButtonBox, 'buttonBox'), SIGNAL('accepted()'),niveau(control,1))
But I got this error: Qt signal is not callable.

Does someone could help me please?


texts: How does ST_SnapToGrid work using nearest neighbor interpolation?

texts:
I have rasters in a PostGIS database that are not necessarily aligned (some rasters offset from the grid). I want to use ST_SnapToGrid to align all rasters to the same grid but have some questions about how it works using nearest neighbors.

Which of these ways is how nearest neighbors works?

If the pixels are aligned then it is just the value of the source pixel, but otherwise the target pixel takes on the value of the closest source pixel. This would be akin to "moving" the pixel into the correct place.
The target pixel takes the average of the nearest neighbors, such that it might overlap 99% with source pixel A and 1% with source pixel B but take the averages of these pixels to be the value of the target pixel.
Some other way?


texts: Using PostGIS to find which county in California has most settlements in forest?

texts:
Using PostGIS how do I find which county in California has the most settlements in the forest?
I have point settlements settlements_CA, polygon forests forests_CA and polygon counties counties_CA in PostGIS Database.
I tried this:

SELECT COUNT(st.*) as count, ct.name
FROM settlements_CA AS st, forests_CA as fo, counties_CA AS ct
WHERE
ST_Contains(fo.geom, st.the_geom)
GROUP BY ct.name
ORDER BY ct.name
;

this query was very slow and the result is wrong because it reports that all counties have the same numbers of settlements.

I also tried JOIN and LEFT JOIN but failed.


texts: Retrieving raster/coverage width and height in pixels from GeoServer

texts:
Is it possible to retrieve the original width and height in pixels for a coverage in GeoServer using either the Python API or REST?  I could not find anything in the documentation regarding how to retrieve height or width of a coverage in GeoServer.

Here is my current solution interfacing with GeoServer.  It needs GDAL, and it will probably fail if the coverage is not a simple file readable by GDAL (i.e. a PostGIS raster)

import gdal
import os
from geoserver.catalog import Catalog

#...code

#store is a GeoServer store
def get_raster_dimensions(store):
    """ Get the size of the raster from the catalog. """
    gdal.UseExceptions()
    #this, of course will fail if the data isn't a regular file
    base_url = store.url.replace('file:', '', 1)
    #what if the GEOSERVER_DATA_DIR is unset?  How to find it then?
    coverage_path = os.path.join(os.environ['GEOSERVER_DATA_DIR'], base_url)
    #if the coverage store is a raster readable by GDAL,
    #we can get the dimensions with GDAL
    try:
        dataset = gdal.Open(coverage_path)
        return (dataset.RasterXSize, dataset.RasterYSize)
    except RuntimeError:
        return (None, None)
    finally:
        #no close method in GDAL docs, this was referenced as the way
        #to clean up the file after finished
        dataset = None


texts: How can I calculate the center point inside a polygon in ArcGIS 9.3?

texts:
I'm using the JavaScript API 1.6 on ArcGIS 9.3. I have a polygon geometry in an ArcGIS map, and I'm trying to find a way to get the point that lies in the exact center of the geometry. What is the best way to do this?

Thanks!


texts: How to determine correlation between two databases?

texts:
I need to perform some tests and verifications on databases and/or datasets.
I need for example to analyze two databases and then determine if there is enough correlation between both.
My main issue is how to deal with correlation...
Do I need to dice/grid the map and analyze precisely each square?
Is it possible to compute statistical techniques (i.e. mean, median...) to achieve this goal or not at all?
My datasets are essentially made of GeoTIFF files for elevation terrain information and Shape files for my objects (roads, trees...).
I had a look on some APIs to work with those kind of files:

Shapefile C Library
PySAL
GDAL - Geospatial Data Abstraction Library

Finally I read in some research papers about softwares determining databases correlation, but couldn't find any of them:

LightBox SAIC
ZCAP

Thanks for help.


texts: Counting area percentage of overlapping polygons

texts:
I have a set of points (meteorological stations) with 5km buffers around them. Additionally, I have a set of polygons (municipalities) in which those stations are located.

I would like to have data on how much of the buffer is filled by each polygon but I was not able to find a function or a way to do it. So for example for point A I would have that in the buffer of 5km there is 80% of polygon B and 20% of polygon C as the output.

Edit:
The anwsers provided do anwser this question in a setting in which I know what polygons and buffers I need to have in the code or when I can do it manualy. Problem is that in my setting with 1 shapefile having 500 polygons and one vector layer having 13000 buffers I am unable to do this manually and I am asking if there is an easier way to do it. I tried implementing the approaches in the questions but I was unable to use them on my example in big scale. Moreover I had problems with coding that I also noted in comments down below.

Main difference to the other similar question provided is that I have 2 vector layers (Buffer with 500 buffers and Municipalities a shape file with 500 different polygons). Previous questions were analyzing this question in situations in which I know precisely the polygons and buffers analyzed, at least that is how I understood them, so after trying to use those approaches on my datasets I was not able to replicate them.

Another way of approaching this problem for me is to use join attributes by location -one-to-many and then getting averages of the duplicate results, only problem is that then I take average of polygons instead of their weights corresponding to % of intersection with the buffer. In my case I do not know what polygons of the shapefile are going to intersect what buffers, or even if the buffer is going to be intersected by 1,2 or more polygons


texts: Display DegreesMinutesSeconds on status bar

texts:
I have used status bar to display x and y coordinates in esriDirectionUnits.esriDUDegreesMinutesSeconds form.

I have created the code like this

statusBarXY.Text = string.Format("{0}, {1}  {2}", e.mapX, e.mapY, esriDirectionUnits.esriDUDegreesMinutesSeconds);

but it will not display any format like this X=70°0'20.797"E  Y=  22°28'.072" N
                                                      X= 700'20.797"E                                     Y= 2228'.072"N

How could I make this?  If any one has a solution please send me link or a sample.


texts: How to merge overlapping and duplicate lines from vector layer to one single line?

texts:
I used the "polygon to line" tool from QGIS, to get a line-layer from this polygon-layer:

Unfortunately this causes overlapping and duplicate lines, after I broke the qgis-export via GRASS into primitive lines (2 vertices). I need to have just one line, because two identical or overlapping lines are causing problems in CityEngine.
I found no way to identify those features that are within other longer lines. I know, in QGIS there is a function called "within", but i have no idea if it works with line layers or how to use it.

many thanks in advance


texts: subset FeatureCollection by position - Google Earth Engine

texts:
I have a FeatureCollection such as:

// Make a list of Features.
var features = [
  ee.Feature(ee.Geometry.Rectangle(30.01, 59.80, 30.59, 60.15), {name: 'Voronoi'}),
  ee.Feature(ee.Geometry.Point(-73.96, 40.781), {name: 'Thiessen'}),
  ee.Feature(ee.Geometry.Point(6.4806, 50.8012), {name: 'Dirichlet'})
];

// Create a FeatureCollection from the list
var fromList = ee.FeatureCollection(features);

How can I extract Features from the FeatureCollection using their index position within the FeatureCollection?

For example, fromList[0] would return the Feature named "Voronoi".


texts: I want to consume a erdas apolla WMS service using ArcGIS Javascript API on my GIS application

texts:
Where do I type in the username and password.

<!DOCTYPE html>
<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8">

<meta name="viewport" content="initial-scale=1, maximum-scale=1,user-scalable=no">
<title>Map with WMS</title>
<link rel="stylesheet" href="https://js.arcgis.com/3.20/dijit/themes/claro/claro.css">
<link rel="stylesheet" href="https://js.arcgis.com/3.20/esri/css/esri.css">
<style>
    html,
    body,
    #map {
    height: 100%;
    width: 100%;
    margin: 0;
    padding: 0;
    }

    body {
    background-color: #FFF;
    overflow: hidden;
    font-family: "Trebuchet MS";
    }
</style>
<script src="https://js.arcgis.com/3.20/"></script>

<script>
    var map;
    require(["esri/map", "esri/layers/WMSLayer", "esri/config",
        "dojo/domReady!"
    ],
    function(Map, WMSLayer, esriConfig) {
        esriConfig.defaults.io.proxyUrl = "/proxy/";

        map = new Map("map", {
        basemap: "streets",
        center: [-98, 37],
        zoom: 5
        });
        var wmsLayer = new WMSLayer(
        "http://103.203.139.72/erdas-apollo/coverage/District_Jai?REQUEST=GetCapabilities&SERVICE=WMS&VERSION=1.3.0", {
            format: "png",
            visibleLayers: [2]
        });

        map.addLayer(wmsLayer);
    });
</script>
</head>

<body>
<div id="map">
</div>
</body>
</html>


texts: Where can I find all Styled Layer Descriptor vendor options for geoserver?

texts:
I am working with GeoServer and use SLDs to style my data. I noticed some tags you can use in your SLD which seem to be GeoServer specific.
Is there a site which lists the different vendor options which are available for GeoServer? I have seen the cookbooks on the GeoServer site, but do they list all the different options?


texts: ERROR Validate Network Topology tool in ArcGIS pro

texts:
I get the error in the image that I added dem image in the Validate Network Topology tool in ArcGIS Pro program:

ERROR 001793: Network topology must be enabled. Failed to execute
(ValidateNetworkTopology)

How can I fix this error?


texts: QGIS Dissolve tool erroneously fills rings in polygons

texts:
I have a shapefile layer with >27,000 polygons. I want to dissolve these together to create a footprint layer. The Dissolve tool works fairly well, but it has filled in what should be rings in the new layer. (Same result in both QGIS 2.14 and 3.4). They are gaps in the original layer, surrounded by multiple polygons, (i.e. these are true spaces, not rings in a single polygon, and they should not be filled). Can anyone suggest a way to avoid this happening, please? Or a workaround to fix the layer now that the rings have been filled?

Starting layer;

Post Dissolve tool - with rings erroneously filled

Desired output overlaid on the first dissolved layer showing the filled rings.

Edit

I have now managed to get it to work, so perhaps it was a bug or an issue with my input layers? Would still like to know why it was happening to avoid a repeat in the future.


texts: Postgis Spatial query using pyqis between two spatial table

texts:
I try to query two spatial table using st_within function, then I want to display the query result as a QGIS Layer. But the layer did not add to the QGIS Canvas.
Here my code :

conn = psycopg2.connect("dbname='db_pbb' host=localhost port=5433 user='postgres' password='septin'")
    cur = conn.cursor()
    cur.execute("""SELECT wa FROM gis.kecamatan""")
    rows = cur.fetchall()
    row_list = []
    for row in rows:
        #row_list.append(row.)
        self.dlg.cboKecamatan.addItems([row[0]])
    self.dlg.show()
    # Run the dialog event loop
    result = self.dlg.exec_()
    text = str(self.dlg.cboKecamatan.currentText())
    txtnop = self.dlg.lineNOP.text()

    # See if OK was pressed
    if result == 1:
        #Do something useful here - delete the line containing pass and
        # substitute with your code.
        #Versi satu - sudah berhasil, yeay, you know the key right ?

        sql = "(SELECT gis.tm_bidang.numb, st_within(gis.tm_bidang.geom,gis.kec_palu_utm.geom) as within_geom, gis.tm_bidang.d_nop,gis.tm_bidang.d_luas, gis.tm_bidang.aktor from gis.tm_bidang,gis.kec_palu_utm where gis.kec_palu_utm.wa = '"+text+"' and st_within(gis.tm_bidang.geom,gis.kec_palu_utm.geom))"
        uri = QgsDataSourceURI()
        uri.setConnection("localhost", "5433", "db_pbb", "postgres", "septin")
        uri.setDataSource("", sql, "within_geom", "", "numb")
        tmadminkec = QgsVectorLayer(uri.uri(), "tm_bidang", "postgres")
        QgsMapLayerRegistry.instance().addMapLayer(tmadminkec)

Any ideas ?,


texts: Setting annotation feature character spacing?

texts: I'm having an issue with an annotation feature in ArcGIS desktop 10.5 where the character spacing is not being applied to the last word in multi-word features (see "ocean" in Pacific Ocean below) when exported to pdf (displays fine in layout view). Other multi-word features don't have this issue. Only features where character spacing <> 0. Any ideas on troubleshooting this issue?


texts:
I am making a map but my histogramm doesnt seems to be good.
 sgf_var <- "V"
tab <- "REC_T234"

label <- labels_rec %>%
 filter(VAR == sgf_var & TABLEAU == tab)

 map <- tm_shape(sgf_polyg_t234) +
     tm_polygons(col = "V8", style = "jenks", title="sgf_var",
     legend.hist = TRUE,
     legend.format=list(fun=function(x) formatC(x, digits=0, format="f"), text.separator= "-")) +
     tm_borders("#9E9E9E") +
 tm_layout(main.title = label$LIBELLE, main.title.size = 0.5,legend.outside = TRUE,  legend.title.size = 0.7, legend.hist.width = 1)

Here the result :
   map


texts: GeoPandas ValueError: A LinearRing must have at least 3 coordinate tuples

texts:
I am trying to read a multi-polygon GeoJSON file in GeoPandas, but I am getting the following error:

'ValueError: A LinearRing must have at least 3 coordinate tuples'

command used is gpd.read_file(layer_file)

What can I do to filter bad geometry? I converted the dataset from pbf to GeoJSON before reading with GeoPandas. The same GeoJSON file works with QGIS.


texts: Splitting single line based on 3D length?

texts:
I have a transect line to which I have already added the z value using the Interpolate Shape tool. I need to split this line into 100m in 3D length, which due to elevation must be different than 2D. However the split tool in 3D editor only gives me the 2D length.

Is there anyway to resolve this?

It will be great if I can do it in batch process.


texts: Error selecting geometry from Oracle DB with ArcSDE

texts:
I'm trying to select a geometry column as well-known text from a feature class in an Oracle 11g/ArcSDE 10.1 geodatabase. (This is in a Python script using pypyodbc as a driver). This is my SQL:

SELECT SDE.ST_AsText(SHAPE) FROM OWNER.FEAT_CLASS

When I run this I get a few Oracle errors:

ORA-29900: operator binding does not exist\nORA-06553: PLS-306: wrong number or types of arguments in call to 'ST_ASTEXT

I found this support article which says to fully qualify all references to ST_GEOMETRY operators, but I have the SDE prefix so I'm not sure what's missing. Does anyone have any ideas regarding this?

UPDATE: If I just select the SHAPE field without any SDE functions I get Decimal('214'). Interestingly, that's also the object ID.


texts: Dump all spatialite tables to shp and dbf

texts:
Given a spatialite database i want to export/dump all tables to shapefile (or dbf for those without geometry).

I want to do it programatically so acceptable solutions will be using a (bash) command line tool like spatialite, ogr2ogr, ... or Python.

Is there a command that export all database in one shot or should i iterate over each layer and dumping it individually?


texts: Selecting south-facing slopes with slope > 10° using ArcGIS Desktop?

texts:
I have to make a Model in the Model Builder which will select the areas with slope bigger than 10° and which are facing to the south. It should be used for environmentalist who will just insert a DEM and the model will automatically display only these wanted areas.

How do I do that?

I already used tools Slope to and Aspect, but I don't know how to select only that values which I want to display and make a Model for that.


texts: Correcting incorrect altitudes of points recorded by my Android or how to calculate Geoid Height in QGIS

texts:
I tried to edit my question with more details and I also think that lot of people can benefit from this question!

I was recording a track with my Android and later I found out that I end up having incorrect altitudes. To be precise I was recording a terrain of an island situated in Croatia.
First I started asking myself why is that, and I got answer right here: Inaccurate altitude GPS data with my Android? as well on some other pages, including one also mentioned in the linked above: stackoverflow.com/questions/11168306/.
So it turns out that I recorded ellipsoidal height(h) and in order to get Orthometric height(H) (real height that I need), I need to calculate first Geoid Height(N), which can be seen here:

[source: http://www.geomobileinnovations.com/]

So I have manually calculated Geoid Height from this page: http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm96/intpt.html and as alternative for not transforming decimal to degrees from third party there is also this page: https://geographiclib.sourceforge.io/cgi-bin/GeoidEval?input=44.97775%2C+14.73826&option=Submit.

Now I am wondering how to calculate geoid height automatically in QGIS?
People were mentioning that there is some sort of documentation for calculating geoid height on earth-info.nga but it is given in Fortran/C/?. In addition to that Google Earth have option that pops up on opening a file: Adjust altitudes to the ground level, which I assume it might be a possible solution for Google. So if Google Earth has it I guess there is some sort of plugin from QGIS too.

Mobile phone that I have: Sony Xperia m4 aqua and I had High accuracy mode turned on. Specification for my phone can be found here: imei.info/phonedatabase/16958-sony-xperia-m4-aqua-e2303 (GPS: A-GPS/Glonass, GLONASS: No).

EDIT

My question was marked as duplicate, if it refers to my old question, I dont see it to be the same. My last question was about "why", this question is about "how".


texts: Enable WMS GetFeatureInfo request in Leaflet using R?

texts:
I've been using R with Leaflet and Shiny to produce a web map that is accessible to the public online. It is a real-time air temperature map of Tasmania, Australia based on current Bureau or meteorology weather station observations. You can view a prototype here: http://144.6.226.103/shiny/

Each air temperature map is referenced from GeoServer via WMS. I have produced the R code below for reference.

I'm trying to enable the GETFeatureInfo function so users are able to query each map spatially by clicking on any location thereby producing a popup and showing a temperature value.

I'm wondering if this is possible to set up in R?

library(leaflet)
library(raster)
library(htmlwidgets)
library(leaflet.extras)
    m<-leaflet() %>% addTiles() %>% setView(147, -42, zoom = 8) %>%
            addWMSTiles(
              "http://144.6.226.103:8080/geoserver/TempPred/wms?service=WMS",
              layers = paste0("TempPrediction_2019090511AEST"),
              group = paste0("TempPrediction_2019090511AEST"),
              options = WMSTileOptions(format = "image/png", transparent = TRUE, version ="1.3.0",opacity = 0.5))%>%
            addWMSTiles(
              "http://144.6.226.103:8080/geoserver/TempPred/wms?service=WMS",
              layers = paste0("TempPrediction_2019090510AEST"),
              group = paste0("TempPrediction_2019090510AEST"),
              options = WMSTileOptions(format = "image/png", transparent = TRUE, version ="1.3.0",opacity = 0.5))%>%
            addWMSTiles(
              "http://144.6.226.103:8080/geoserver/TempPred/wms?service=WMS",
              layers = paste0("TempPrediction_2019090509AEST"),
              group = paste0("TempPrediction_2019090509AEST"),
              options = WMSTileOptions(format = "image/png", transparent = TRUE, version ="1.3.0",opacity = 0.5))%>%
            addWMSTiles(
              "http://144.6.226.103:8080/geoserver/TempPred/wms?service=WMS",
              layers = paste0("TempPrediction_2019090508AEST"),
              group = paste0("TempPrediction_2019090508AEST"),
              options = WMSTileOptions(format = "image/png", transparent = TRUE, version ="1.3.0",opacity = 0.5))%>%
            addWMSTiles(
              "http://144.6.226.103:8080/geoserver/TempPred/wms?service=WMS",
              layers = paste0("TempPrediction_2019090507AEST"),
              group = paste0("TempPrediction_2019090507AEST"),
              options = WMSTileOptions(format = "image/png", transparent = TRUE, version ="1.3.0",opacity = 0.5))%>%
            addWMSTiles(
              "http://144.6.226.103:8080/geoserver/TempPred/wms?service=WMS",
              layers = paste0("TempPrediction_2019090506AEST"),
              group = paste0("TempPrediction_2019090506AEST"),
              options = WMSTileOptions(format = "image/png", transparent = TRUE, version ="1.3.0",opacity = 0.5))%>%
            addWMSLegend(uri=paste0('http://144.6.226.103:8080/geoserver/TempPred/ows?service=WMS&request=GetLegendGraphic&format=image%2Fpng&width=20&height=20&layer=TempPrediction_2019090511AEST'))%>%
            addLayersControl(
              baseGroups = c(paste0("TempPrediction_2019090511AEST (°C)"),paste0("TempPrediction_2019090510AEST (°C)"),paste0("TempPrediction_2019090509AEST (°C)"),paste0("TempPrediction_2019090508AEST (°C)"),paste0("TempPrediction_2019090507AEST (°C)"),paste0("TempPrediction_2019090506AEST (°C)")),
              position = 'bottomleft',
              options = layersControlOptions(collapsed = FALSE))


texts: Find LineStrings with common segments in PostGIS 2.3

texts:
I would like to find LineStrings with common segments ('collinear', i.e. LineStrings that have 2 or more consecutive points in common) with PostGIS.

I'm not shure if I'm missing something, but ST_Touches as well as ST_Overlaps seem to return True for 2 geometries that allready share only one common point.

Is there any PostGIS function I'm missing or any other solution to this?


texts: Using multiple conditional factor in QGIS?

texts:
I have a column filled with multiple informations in each row :

> e.g. row 1 :( 'Sécurisation des intersections;Reprise du revêtement;Ajout de signalisation (M12, feu vélo, etc.)' )

> e.g. row 1 :(   'Plateau traversant;Sécurisation des intersections;Aménagement de traversées;Ajout de signalisation (M12, feu vélo, etc.)'  )

I would like to fill another column with the calculation of related value.

I will try to explain. The first field is filled with different type of road work, and it's possible that some sections needed multiple work. Therefore there is multiple type of road work in the same row separate with point comma.

Now I want to affect some price to those road work. I know the value for each type of road work and would like to calculate the sum of thoses different road work for each row.

I know how to write this :

  IF  "MESUACCOMP" LIKE '%Sécurisation des intersections%' THEN 1000
  IF  "MESUACCOMP" LIKE '%Reprise du revêtement%' THEN 5000
  IF  "MESUACCOMP" LIKE '%Ajout de signalisation (M12, feu vélo, etc.)%' THEN 12000

But how could I get a sum?

When those 3 case are encountered in the same line I want to have 1000+5000+12000 = 18000

This is a bit difficult to explain in English.


texts: Option to save as GeoPDF greyed out 3.14.16 MacOS Catalina

texts:
I've been looking through post asking the same question but can't get it work.
This post Why can't I make GeospatialPDFs using QGIS 3.14? mentioned a "PDF driver there for creating Geospatial PDF's" under Preferences > Options > GDAL but I can't located it.
I have the exact same about screen as in the post.
Has this issue been solved yet?


texts:
Trying to compute chart for NDWI

var dataset = ee.ImageCollection('LANDSAT/LE07/C01/T1_8DAY_NDWI')
                      .filterDate('1999-01-01', '2002-12-31')
                      .filterBounds(mumbai);
    var colorized = dataset.select('NDWI');
    var colorizedVis = {
      min: 0.0,
      max: 1.0,
      palette: ['0000ff', '00ffff', 'ffff00', 'ff0000', 'ffffff'],
    };
    Map.setCenter(72.88, 19.154, 10);
    Map.addLayer(colorized, colorizedVis, 'Colorized');
    print(ui.Chart.image.series(colorized, mumbai));

Error:

Error generating chart: Error in map(ID=19990525):
Image.reduceRegion: The default WGS84 projection is invalid for aggregations. Specify a scale or crs & crs_transform.


texts: Looking for the Coords of US Map Boundaries

texts:
I am writing a program that needs a list of coordinates for the boundaries of the US. I can only find states individually but not the whole country. This is an example of the boundaries for Alabama.

{"Alabama": [coordinates: [35.0041, -88.1955],
                             coordinates: [34.9918, -85.6068],
                             coordinates: [32.8404, -85.1756],
                             coordinates: [32.2593, -84.8927],
                             coordinates: [32.1535, -85.0342],
                             coordinates: [31.7947, -85.1358],
                             coordinates: [31.5200, -85.0438],
                             coordinates: [31.3384, -85.0836],
                             coordinates: [31.2093, -85.1070],
                             coordinates: [31.0023, -84.9944],
                             coordinates: [30.9953, -87.6009],
                             coordinates: [30.9423, -87.5926],
                             coordinates: [30.8539, -87.6256],
                             coordinates: [30.6745, -87.4072],
                             coordinates: [30.4404, -87.3688],
                             coordinates: [30.1463, -87.5240],
                             coordinates: [30.1546, -88.3864],
                             coordinates: [31.8939, -88.4743],
                             coordinates: [34.8938, -88.1021],
                             coordinates: [34.9479, -88.1721],
                             coordinates: [34.9107, -88.1461]
                            ]
                }

Does anyone have a list for the whole country of the US?


texts: Converting large data with lat and long into X and Y?

texts:
I have 9888562 records in dataframe and I would like to convert my lat, long to UTM x,y. according to my code, I have used pyproj package but because my data are too much it takes a long time and finally, it doesn't work.
I wonder whether you know another way or package that I can use for my data?

def rule(row):
    p = Proj(proj='utm',zone=10,ellps='WGS84', preserve_units=False)
    x,y = p(row["LON"], row["LAT"])
    return pd.Series({"X": x , "Y": y})
My_data = My_data.merge(My_data.apply(rule, axis=1), left_index= True, right_index= True)


texts: GRASS r.reclass QGIS 2.14. error: "the following layers were not correctly generated"

texts:
I am getting a strange error message when running the GRASS r.reclass algorithm.
The error states that:

"The following layers were not correctly generated
       *reclassified
 You can check the log message to find more information about the execution of the algorithm. "

After having read the log-file a couple of times, I am none the wiser.

All that I found were a couple of error messages...
"illegal latitude for north"
"Invalid input"
"Raster map not found"

If any of you have any idea what the reason is why the algorithm is not functioning correctly, please...let me know.
See below for screenshots of log-file and re-class file.

New screen shots. I ran the r.reclass algorithm using different input parameters, but to no avail. (xmin,xmax,ymin,ymax = -180,179,-64,90)


texts:
I am wondering if there is a way to do field value maps that are filtered based on another attribute field?  I know how to associate value maps, but not that are connected to another field.

I have two fields in a table Type and Subtype.  Type will have three values: Trees, Bushes, and Flowers.

The Subtypes will have different values.  ie:

 Subtype for trees would be conifer and deciduous
 Subtype for flowers would be roses, daisies, and sunflowers
 Subtype for bushes would be blueberry, raspberry, and blackberry

What I want is when Trees is selected from Type, only conifer and deciduous are available from subtype; etc...

Thanks...


texts: Are FOSS4G conference presentations available online?

texts:
I'm wondering if anyone has knows of links to a website(s) where previous FOSS4G conference presentations are available online.  I would like to:

watch recorded video presentations,
listen to audio recordings, or
read through detailed PPT presentations (more than just abstracts)

In particular, I am interested in videos from the FOSS4G-North America 2013 conference, but this would be a good archive of locations to find any past FOSS4G conference presentations.


texts: Can I use ST_Within or ST_Contains to update columns based on a geom being within another?

texts:
Based on the following tables:
Street

Municipality

I'd like to populate Street.NOM_ENT with Municipality.NOM_ENT and Street.NOM_MUN with Municipality.NOMBRE any time Street.geom is within Municipality.geom.
In other words, I have a table with all the streets in a country and I need to add to which state and county that street belongs.
I know how to check if Street.geom is within Municipality.geom, but I'm kind of new working with postgis/postgres and I'm not sure how to use this as a condition for an UPDATE query.


texts: How to Find the Nearest Node from a Given Coordinate

texts:
I have already developed a table manually. It contains the coordinates and the distance of all the nodes, there is no geom field or any field related to postgis. But I want to find out the nearest nodes to a given coordinates. Is this possible, or is there any other way to find out?

Kindly help me so that I can proceed further.

Thanks.


texts: Confusion with Postgis Topology and Geometrytype

texts:
in the moment I try to understand the principles of postgis topology. While testing I discovered a behavior, which I don't understand.
I have a table public.test with a standard geometrycolumn type LINESTRING.
Then first I added a topology with:

select topology.CreateTopology(’testtopo’, 25832);

Next I added a column for the topologygeometry to the table test:

select topology.AddTopoGeometryColumn('testtopo','public','test','topogeom','LINESTRING');

At least I made an update to fill the column topogeom:

update test
set
topogeom = topology.toTopoGeom(geom, 'testtopo', 1, 2.0);

Till this point everything works fine. What I now doesn't understand is, that if I test the geometrytype with

select ST_GeometryType(topogeom) from test;

the result is a MULTILINESTRING.
When I query the relatet topology table with

select ST_GeometryType(geom) from testtopo.edge_data;

it is a LINESTRING.
Maybe someone has an explanation for that or is it possible to force the topology to be a LINESTRING?


texts: Reducing monthly averaged MODIS data over several ROIs

texts:
I want to create monthly averages for the MODIS Land Surface Temperature (LST) product for several ROIs. I am able to create an image collection containing the monthly averaged MODIS data, but I have not figured out how to reduce this data for the ROIs.

Here is the first part, with help from here and here:

var modis = ee.ImageCollection('MODIS/006/MOD11A2');
var months = ee.List.sequence(1, 12);

Map.centerObject(areas, 7);
Map.addLayer(areas, {}, 'ROI');

// First I divide the data into monthly averages
var byMonth = ee.ImageCollection.fromImages(
      months.map(function (m) {
        return modis.filter(ee.Filter.calendarRange(m, m, 'month'))
                    .select(0).mean() // selecting band 1: LST_Day_1km
                    .set('month', m);
}));
print('byMonth:', byMonth);

var check = ee.Image(byMonth.first());
Map.addLayer(check, {}, 'check');

Then comes the next part, where I try to use .reduceRegions() to reduce over multiple regions (my FeatureCollection "areas"). The idea is to loop/map through the 12 months, reducing the monthly data from "byMonth" to a value for each ROI. Unfortunately it doesn't work.

var reduced = months.map(function(mm) {
  return byMonth.reduceRegions({
  collection: areas, // a featureCollection containing multiple polygons
  reducer: ee.Reducer.mean(),
  scale: 30,
  maxPixels: 1e9
}).filter(ee.Filter.eq('month', mm));
});
print(reduced);


texts: Image source of world political map for Google Maps overlay in Android?

texts:
I have searched a lot. But I am not able to find some high resolution image of world map showing countries in different colors without label, which I can slice and use to create tile overlay in my Google Maps.

Any idea from where I can get a world image showing countries in different colors without labels?


texts: Error 999999 when converting a point feature class to a raster in python

texts:
After having converted my xyz-file to a point feature class, I am trying to convert to a grid using a Python snippet. The data represents river bathymetry. However, when executing the script, I get: 'ERROR 999999: Error executing function'.

There must be something wrong with my input data, or with the arguments in the python snippet. I already changed double quotes into singles (worked for me with a previous problem), and before converting my .pts-input file to a point feature class, I checked for any odd lines that might contain string - in this file there are none.

This is my code (the converting-to-point-feature-class section does work):

The variable fn is simply the input file name defined earlier on in the script. RDnew is a coordinate system that was also defined earlier on in the script.

    print 'Converting to point feature class.'
    pts_shp = fn[0:8] + '.shp'  #output file name (using fn[0:8] to get rid of the .pts or .txt extension)
    arcpy.ASCII3DToFeatureClass_3d(fn,'XYZ','.\\shp_points\\' + pts_shp,'POINT','1',RDnew,'#','#','DECIMAL_POINT')

    print 'Converting to grid.'
    output_grid = 'grid_' + fn[3:8]
    arcpy.PointToRaster_conversion('.\\shp_points\\' + pts_shp,'Shape.Z','.\\grids\\' + output_grid,'MOST_FREQUENT','NONE','1')

Error 999999 does not really get me anywhere. Do you have any suggestions on what might be wrong with my input / code?


texts: Selecting features by attribute in Leaflet

texts:
Is there a way to select an individual feature based on an attribute from a FeatureGroup in Leaflet?

For example, when selecting a feature based on a click event, I will do something like,

function myFun(e){
    layer = e.target;
    }

Is there a way to do something along the lines of,

function myFun(uniqueID){
    layer = featureGroup WHERE featureGroup.feature.properties.ID == uniqueID;
    }


texts: Counting points (accidents) per segment of roads layer?

texts: The addresses of the road accidents were known and I did geocoding in order to find their position. But their position are not snapped on the road. I managed to snap them to road segments by following Snapping points to lines in ArcGIS Desktop and automate using VBA? . Now, I want the segments to obtain a new field which counts the number of accidents which occured on them. I want sth like the tool "points in polygon". I have heard that linear referencing is helpful but I haven't found something for me.


texts: ValueError: array larger than output file, or offset off edge Python GDAL

texts:
I have been working on this iterator that would count instances where data conformed to a set of conditions over a number of rasters. I had people help me with other aspects of this script (Counting events from a GRIB raster), but now I have hit an error where the return is "ValueError: array larger than output file, or offset off edge", Which I can't seem to be able to work arround. The code is provided bellow. It requires the user to input the path to folder and provide a name for the output GeoTIFF raster. Files read are GRIB.

import gdal
import numpy
import osr
import os

del_mapa = str(input('Vpiši pot do datotek v obliki "C:\\Mapa1\\Mapa2\\...\\Zadnja mapa": '))
os.chdir(del_mapa)

pix_vel = 1000
#Defines pixel size
ime_rez = str(input("Poimenuj datoteko z rezultati: ") + ".tif")
#User can provide a name for the output picture
format = "GTiff"
driver = gdal.GetDriverByName(format)

vir_mere = gdal.Open('inca_20140101-0100.grb')
trans = vir_mere.GetGeoTransform()
geotrans0 = trans[0]
geotrans1 = trans[1]
geotrans3 = trans[3]
geotrans5 = trans[5]
#Extracts transformation data from one of the GRIBs
dst_ds = driver.Create(ime_rez, 401, 301, 1, gdal.GDT_CInt32)
dst_ds.SetGeoTransform([geotrans0, geotrans1, 0, geotrans3, 0, geotrans5])
#and pastes it directly into the newly created GeoTIFF
projInfo = vir_mere.GetProjection()
srs = osr.SpatialReference()
srs.ImportFromWkt(vir_mere.GetProjectionRef())

dst_ds.SetProjection(srs.ExportToWkt())

raster = numpy.zeros((401, 301), dtype=numpy.uint8)
dst_ds.GetRasterBand(1).WriteArray(raster)

def krog(baza,novi):
    x = 0
    while x < 25:
        baza = gdal.Open(ime_rez)
        if x < 10:
            novi = gdal.Open('inca_20140101-%s00.grb') % ('0' + str(x))
            #This should open the next file in line
            if novi > 15 and novi < 25:
                baza = baza + 1
            else:
                baza = baza + 0

        elif x > 10:
            # novi = gdal.Open('inca_20140101-%s00.grb') % (str(x))
            if novi > 15 and novi < 25:
                baza = baza + 1
            else:
                baza = baza + 0
        x += 1
dst_ds.GetRasterBand(1).WriteArray(baza)
dst_ds = None


texts: How to programmatically start drawing in OpenLayers 4._

texts:
In openlayers there are 2 paradigms for drawing features one is to create a "draw interaction" that lets the user draw a feature and the other is to create the geometry and feature yourself--then add it to the layer.

I would like to programmatically start a polygon "interaction" drawing for the user by placing the first point. After which the user will place the remaining points. It therefore must be an "interaction".

This means I must simulate a user click or draw event that the openlayers api accepts as valid.

I've tried..

map.dispatchEvent(<clickevent>)
map.dispatchEvent(<drawstartevent>)

/*draw being ol.interaction.Draw*/
draw.dispatchEvent(<clickevent>)
draw.dispatchEvent(<drawstartevent>)


texts:
I can't seem to find out how to add labels in my popups by using qgis2web.
Where others can add labels to the popup fields (picture below)

I can only enable or disable popups (picture below). This results in popups with only values, which are difficult to understand without the labels.

Using QGIS Desktop 3.4.8.


texts:
I am working with some data that needs to be visualized in different transparencies levels.
I know how to do it in ArcMap: Symobology --> Advanced --> Transparency.
Now I am working a 3D project and I want to do the same only in ArcScene, so you can look through the surface. I tried to do it just like in ArcMap but it does not work.

Anyone has an idea how to do it?


texts: Why is there a need to train the clusterer in Google Earth Engine

texts:
I'm wondering why before clustering the global Landsat-7 image in this tutorial, the authors first trained the clusterer on a subset region around Egypt?

Since the Kmean algorithm is supposed to be unsupervised, I though there was no need to train it beforehand.

The line of code I'm puzzled about is this one:

var clusterer = ee.Clusterer.wekaKMeans(15).train(training);

Is the train() method above the equivalent of the method fit() in Scikit-learn?


texts: Importing colour map from .text file which will then correspond to geologic shapefile on QGIS 3.10.6

texts:
I have been searching for this solution for over a week now and have tried many of the suggested solutions online but so far northing has worked. I have a Geology vector layer and a corresponding RGB text file which has the appropriate colour values.
I have tried using the Setting > Style Manager > Import. This requires an xml file which I made from the text file but I am met with the following error when I attempt to import it.
Incorrect root tag in style: w:wordDocument

I have also attempted using the Style tab in the Symbology section of the Layer Properties but this method requires a qml or sld file.
I have attempted to import the txt file to the layer tree and then save the style as a qml and import the saved style. I am met with the following error.
Cannot apply style with symbology to layer with different geometry type

The geology shapefile geometry is Polygon (MultiPloygon)
I have spent a few days attempting the various Python solutions but I have zero experience with Python and those attempts have gotten me no where.
Here is an example of the text file I am uploading:
Value,Red,Green,Blue
63,   227, 215, 164
61,   90,  134, 103
60,   109, 167, 133
59,   122, 181, 139

Can anyone shed some light on this issue?


texts: Can Bing base layers be used with the Mapfish / Geoserver print plugin?

texts:
After being able to solve most of my problems just be lurking, I find I now need to come into the light for this one.  The headline is pretty explicit.  I have been unable to use the mapfish print utility to generate my PDF's including the Bing base layers. For a live example, please visit http://gis.lhtac.org/safety_print .  Steps to reproduce follow:

 1.  Select a region.  Click on one of the colored areas.  I suggest district 2.
 2.  Select a Jurisdiction.  A panel will pop with a list of road jurisdictions.  I suggest 'Grangeville'.
 3.  Please observe the layers, including the Bing base layers, and our in-house WMS layers.  Also, observe the 'Print Map' button is now live in the Crash Statistics panel.
 4.  If desired, drag a bounding box to select some of the indicated crashes. This will select some datapoints for a WFS layer derived from the Crash Location WMS layer data.
 5.  Click on the 'Print Map' button.
 6.  Observe the GeoExt print panel, and that the Bing Layers are shown, but possibly only around the edges.
 7.  Zoom in once or twice and observe the Bing layer is now fully represented.
 8.  Click the 'Create PDF' anchor in the tbox of the print panel.  Wait for a bit while the server thinks about it.  Presuming nothing has broken, a generated pdf will be offered for download.
 9.  Open the pdf, and please observe the following:

 A: the Bing base layer is missing.
 B: our in-house WMS layers are shown, although at poor resolution.
 C: If WFS selections had been made, they are shown, with the known filled polygon rendering bug.
 D: Other artifacts are shown, pretty much as might be expected.

My best guess as to the cause of no base layer is that there is some unspecified problem using Bing with the Mapfish utility.  I have been unable to find any reference to this, however, either as a bug, or as a terms of service issue.  The second most likely problem is my config.yaml is incorrectly referencing the Bing servers. That section of the config.yaml follows:  (The SE editor has mangled it slightly)

hosts:
  - !ipMatch
    ip: 127.0.0.1
    port: 8080
  - !ipMatch
    ip: gis.lhtac.org
  - !ipMatch
    ip: www.camptocamp.org
  - !dnsMatch
    host: maps.live.com
    port: 80
  - !dnsMatch
    host: virtualearth.com
    port: 80

Thanks for looking.


texts: Determining direction of traffic in road network?

texts:
I'm working with road segments as LineStrings in GeoPandas, for each segment I would like to determine the direction of travel (they drive on the right).

This information is not encoded, how can I determine it?

I was thinking to query OpenStreetMap, since when I overlay in QGIS I can see this has travel direction information.

Can you think of any other approaches, or guide me in querying OpenStreetMap?


texts: Using Landsat 8 for quantitative analysis - what DNs does GRASS require?

texts:
Hoping someone can short circuit a lot of exploration

Landsat 8 values are stored as unint16 since the sensors work on a range 0 to 4096 whereas the L7+TM was 8 bit and had a digital number range 0 to 255.
USGS link

Looking at the unprocessed bands in QGIS, the value tool reports values of ~8000 for band 3, whereas GRASS reports values of ~0.19 . For the same bands in L7+TM band 3 values are ~30

Understandably, the output from my old workflow in GRASS is giving nonsense. i.landsat.toar does read L8 but with the funny numbers I don't expect it is on the right track.

The end use is NDVI, I don't expect a pansharpening workflow to be the same.

Any tips on a new workflow?


texts: Getting errors using routing function in QGIS OSM Tools plugin

texts:
I thought I found just the right tool for creating routes from all points in one layer to all points in another layer: The "Routing" function in OSM Tools. Unfortunately I am receiving errors trying to execute it.

I am getting the following messages:
INFO: "Limit exceeded : Request limit of 40 per minute exceeded. Wait for 0.003889799118041992 seconds" like ten times.

Then there is the Python warning:

traceback:
warning:sys:1: ResourceWarning:
unclosed <ssl.SSLSocket fd=4060, family=AddressFamily.AF_INET, type=SocketKind.SOCK_STREAM, proto=0, laddr=('10.20.91.24', 8643), raddr=('129.206.7.231', 443)>

And the Python error:

Traceback (most recent call last):
              File "C:/Users/me/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\OSMtools\osm_tools.py", line 145, in run
                route.directions_calc()
              File "C:/Users/me/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\OSMtools\directions.py", line 137, in directions_calc
                layer_out = self._addLine(responses, values_list)
              File "C:/Users/me/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\OSMtools\directions.py", line 183, in _addLine
                values_list[i][0],
            IndexError: list index out of range

Any ideas what could be the problem(s)?

Operating System is Windows 10, x64
QGIS 3.0.0.


texts: How to create a raster from xyz file format

texts:
I have xyz points. I have imported these XYZ points into QGIS and need to create a raster from these points. I have used the raster interpolation plugin and I used IDW as the interpolation method. I got this raster (The following snapshot). Why it is like this? How I can fix it?

The average distance between points are 15m. In the above snapshot, the distance coefficient p is 2.

In the below snapshot, the distance coefficient p is 15. I have to mentioned that I used the IDW method.

What distance coefficient p. I need to use.

Any Idea?


texts: Can custom map tiles be consumed via ArcObjects within my ArcGIS 10 Desktop AddIn?

texts:
Can I programmatically consume and render custom map tiles in ArcMap via our ArcGIS 10 Desktop AddIn?

The AddIn is built using the ArcObjects SDK for C#/.NET and ArcGIS 10. Pulling and rendering map tiles is NOT the AddIn's primary function. But we would definitely like to implement this feature if it can be done.

More specifics:

I received a tip that the ITiledMapServer2 interface is involved but it is unclear to me, and there is a dearth of info in ArcObjects docs regarding this task. The equivalent class within the ArcGIS for Flex API is TiledMapServiceLayer, but again this is an ArcMap AddIn in the desktop space, not a web app.
We have in-house software to generate our own tiles in the same zoom, tile Y, tile X tiling scheme akin to the Google maps JS/Flash APIs.
We are using Apache to serve the map tiles and thus are currently not using ArcGIS Server in any way (i.e. authoring the tiles, publishing as an ArcGIS Server map service, etc).
I've posted the question to SO and received a recommendation to post here.

Thank you in advance.


texts: Compute fractions of feature in grid cell

texts:
I have a shapefile containing features (line) and I have a UTM grid in the same area. The feature is a coastline with different types.
Now I would like to produce the length of coastline and the length of each type within each UTM grid cell.
I prefer using Python or QGIS/SAGA GIS.

I found the post on how to create my UTM grid and I have the shapefile but how do I compute the fractions?


texts: Run mapbasic 32 bit mbx from larger mapbasic 64bit project

texts:
I need to run a 32-bit mbx within a 64-bit mbx -is this possible?

Basically within a larger 64bit mbx routine, I need to call a mapbasic 32 bit mbx so that I can apply some regex code to a table.

Unfortunately Mapbasic regex only works in 32bit currently -https://nyalldawson.net/2013/05/regular-expressions-in-mapbasic/

I know I can add the following code

Run program "performRegex.mbx"

but as the system default is 64bit, how do I force it to use 32 bit?


texts: Memory Error on QGIS rastercalculator from console

texts:
I am experiencing a weird problem using QGIS on Mac.

I have been using RasterCalculator from the python console until yesterday without any problem, but yesterday my QGIS crashed while running a script (well, actually the whole computer crashed), and since then rastercalculator from the consol does not work any longer.

Whatever operation i do, even the most basic ones, return 5 = MemoryError.

Note that RasterCalculator still work without problems if I run it from the raster menu.

I tried to reinstall QGIS twice, but it didn't help.

it looks like something got broken with the crash, or some cache got filled up..

Any idea how to solve this issue?


texts: Extract one WDPA polygon after wdpa_fetch for entire country

texts:
I want to extract one WDPA polygon at a time after downloading protected area data for the entire Democratic Republic of the Congo (IS0 3 = COD) using the wdpar package. Preferably, I would like to extract the polygon by specifying a value in the WDPAID (ex. 124389 is the ID for Okapis Wildlife Reserve). How do I do this? I tried "plot((cod_pa_data[cod_pa_data$WDPAID == 124389]))" but it did not work properly.
library(wdpar)
library(dplyr)
library(ggmap)
cod_raw_pa_data <- wdpa_fetch("COD")
cod_pa_data <- wdpa_clean(cod_raw_pa_data)

cod_pa_data

Simple feature collection with 47 features and 28 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: 1179127 ymin: -1419811 xmax: 2914336 ymax: 686464
CRS:            +proj=cea +lon_0=0 +lat_ts=30 +x_0=0 +y_0=0 +datum=WGS84 +ellps=WGS84 +units=m +no_defs
precision:      1500
# A tibble: 47 x 29
   WDPAID WDPA_PID PA_DEF NAME  ORIG_NAME DESIG DESIG_ENG DESIG_TYPE IUCN_CAT INT_CRIT MARINE
    <dbl> <chr>    <chr>  <chr> <chr>     <chr> <chr>     <chr>      <chr>    <chr>    <chr>
 1 5.56e8 5555120… 1      Yang… Yangambi  Rése… Biospher… National   Ia       Not App… terre…
 2 5.56e8 5555120… 1      Bomu  Bomu      Rése… Wildlife… National   Ib       Not App… terre…
 3 1.67e5 166889   1      Viru… Virunga   Parc… National… National   II       Not App… terre…
 4 1.08e3 1083     1      Gara… Garamba   Parc… National… National   II       Not App… terre…
 5 1.08e3 1079     1      Upem… Upemba    Parc… National… National   II       Not App… terre…
 6 1.08e3 1080     1      Maiko Maiko     Parc… National… National   II       Not App… terre…
 7 1.08e3 1082     1      Kahu… Kahuzi-B… Parc… National… National   II       Not App… terre…
 8 1.08e3 1084     1      Kund… Kundelun… Parc… National… National   II       Not App… terre…
 9 5.56e8 5556978… 1      Salo… Parc Nat… Parc… National… National   II       Not App… terre…
10 3.70e4 37043    1      Okap… Okapis    Rése… Wildlife… National   II       Not App… terre…
# … with 37 more rows, and 18 more variables: REP_M_AREA <dbl>, REP_AREA <dbl>, NO_TAKE <chr>,
#   NO_TK_AREA <dbl>, STATUS <chr>, STATUS_YR <dbl>, GOV_TYPE <chr>, OWN_TYPE <chr>, MANG_AUTH <chr>,
#   MANG_PLAN <chr>, VERIF <chr>, METADATAID <int>, SUB_LOC <chr>, PARENT_ISO <chr>, ISO3 <chr>,
#   GEOMETRY_TYPE <chr>, AREA_KM2 <dbl>, geometry <MULTIPOLYGON [m]>


texts: Spliting a polyline object by line ID based on @data values

texts:
I am trying split a SpatialLinesDataFrame by attributes I merged into @data

I extracted a list of line ID's associated with a points data frame and aggregated them by category and line ID in dplyr

# A tibble: 1,430 x 3
# Groups:   Crime type [14]
   `Crime type`          nearest_line_id                          n
   <chr>                 <fct>                                <int>
 1 Anti-social behaviour 00206205-7878-4F17-BFF3-6086683F041B     1
 2 Anti-social behaviour 00937739-39C7-47E6-888A-4890FA2D47D8     1
 3 Anti-social behaviour 0108B4DD-049A-4A1D-A7CD-8F216AA0DB06     1
 4 Anti-social behaviour 018FB4F0-2DB1-48C4-A627-52A358B99E34     1
 5 Anti-social behaviour 025D1F38-C552-4D84-AC0D-6A1178C10B0F     1
 6 Anti-social behaviour 04542A56-3E24-4F25-AA64-CA874E7E3B5E     1
 7 Anti-social behaviour 04ACF5BC-78A0-475A-9294-F878E53A8625     1
 8 Anti-social behaviour 04AE68AE-E0A9-42B2-908E-29F4EC92051D     1
 9 Anti-social behaviour 05068729-E1D9-4AC6-86BA-06B126395516     1
10 Anti-social behaviour 06E16FFC-F03D-46D3-89A2-E20DDB19DA1C     1

I then merged this dataframe into the @dataattribute of my SpatialLinesDataFrame

Clipped_Street_2@data<-sp::merge(Clipped_Street_2@data,Points_agg,by.x="identifier",by.y="nearest_line_id", all=TRUE )

Now I want to split SpatialLinesDataFrame by the Crime type column so that I can use it in leaflet map control function using the following script.

map<-leaflet(Clipped_Street_2) %>%addTiles()

names(Clipped_Street_2) %>%
  purrr::walk( function(df) {
        map <<- map %>%
      addCircles(data=Clipped_Street_2[[df]],
                 lng=~Longitude ,
                 lat=~Latitude,
                 weight=~n,
                 label=~df,
                 color = ~pal(`Crime type`),
                 group = ~df)
  })

%>%
  addLayersControl(
    overlayGroups = names(Clipped_Street_2),
    options = layersControlOptions(collapsed = FALSE))

My problem is the splitting doesn't appear to work. It does for simple point data sets and I was wondering if there is a function, method or tool I am missing to split this dataset. I am currently working with a simple split function but I get an error (see below).

Clipped_Street_2 <- split(Clipped_Street_2@data, Clipped_Street_2@data$`Crime type`)

Error in split.default(x = seq_len(nrow(x)), f = f, drop = drop, ...) :
  group length is 0 but data length > 0

Not really sure what to do next, but I will post any progress I make myself. This solution seemed sensible, but I couldn't make it work in my situation.


texts: Transferring coordinates unit from Imperial to metric in Arcmap

texts:
I have a raster with NAD_1983_StatePlane_Indiana_East_FIPS_1301_Feet coordinate.

Is there any way to transfer the coordinate to metric in ArcMap?


texts: Overpass Python module query for a way, returns valid response but empty LineString coordinates field

texts:
I use the following query with Overpass Python module, for getting the features of a road in overpass turbo:

[out:json]
[timeout:25]
;
(
  node["name"="John F Foran Freeway"];
  way["name"="John F Foran Freeway"];
  relation["name"="John F Foran Freeway"];
);
out body;
>;
out skel qt;

The response contains a valid array of points and lines, but when I look inside the LineString objects the coordinates field is empty.

The geometry coordinates for the points are correct.

Any idea about this behaviour?


texts: QGIS raster calculator: difference gives NaN for Landsat 8 scene

texts:
I'm strugling with the QGIS Raster Calculator, trying to get a NDVI image from a Landsat 8 scene (path 222, row 84, date 2015-11-15). The problem is, I'm consistently getting NaN values for all the resulting pixels.

As far as I can tell, everything I'm doing should be right: projections, data type, dimensions, origin, etc, are the same (which is expected from a landsat scene)(see below). The formula used is just

(B5 - B4) / (B5 + B4)

or, expressed in terms of my layers:

( "LC82220842015319LGN00_B5@1" - "LC82220842015319LGN00_B4@1" )  /  ( "LC82220842015319LGN00_B5@1" + "LC82220842015319LGN00_B4@1" )

As I mentioned before, all I get is a raster layer with all pixels NaN.

A couple of observations:

I had no problems with the landsat 8 scene for row 222 path 83, for the same date.
The NaNs somehow occur because of the numerator, but not the denominator... In other words: B5 + B4 works good, but B5 - B4 gives NaN.

So, it's really just the difference operation that works badly.

Metadata from the layers:

B5:

Driver
GDAL provider
GTiff
GeoTIFF
Dataset Description
/home/jmb/Ramsar_2016/SIG/Landsat/222_084/LC82220842015319LGN00/LC82220842015319LGN00_B5.TIF
AREA_OR_POINT=Point
Band 1
STATISTICS_MAXIMUM=31568
STATISTICS_MEAN=7466.8914366775
STATISTICS_MINIMUM=0
STATISTICS_STDDEV=7611.703054711
Dimensions
X: 7911 Y: 7981 Bands: 1
Origin
646785,-3.71428e+06
Pixel Size
30,-30
No Data Value
*NoDataValue not set*
Data Type
UInt16 - Sixteen bit unsigned integer
Pyramid overviews
Layer Spatial Reference System
+proj=utm +zone=21 +datum=WGS84 +units=m +no_defs
Layer Extent (layer original source projection)
646785.0000000000000000,-3953715.0000000000000000 : 884115.0000000000000000,-3714285.0000000000000000
Band
Band 1
Band No
1
No Stats
No stats collected yet

B4:

Driver
GDAL provider
GTiff
GeoTIFF
Dataset Description
/home/jmb/Ramsar_2016/SIG/Landsat/222_084/LC82220842015319LGN00/LC82220842015319LGN00_B4.TIF
AREA_OR_POINT=Point
Band 1
STATISTICS_MAXIMUM=24844
STATISTICS_MEAN=4926.3172446242
STATISTICS_MINIMUM=0
STATISTICS_STDDEV=3631.9117068848
Dimensions
X: 7911 Y: 7981 Bands: 1
Origin
646785,-3.71428e+06
Pixel Size
30,-30
No Data Value
*NoDataValue not set*
Data Type
UInt16 - Sixteen bit unsigned integer
Pyramid overviews
Layer Spatial Reference System
+proj=utm +zone=21 +datum=WGS84 +units=m +no_defs
Layer Extent (layer original source projection)
646785.0000000000000000,-3953715.0000000000000000 : 884115.0000000000000000,-3714285.0000000000000000
Band
Band 1
Band No
1
No Stats
No stats collected yet

Update 1: my version of QGIS when first asked this question:

QGIS version    2.8.6-Wien  QGIS code revision  exported
Compiled against Qt 4.8.7   Running against Qt  4.8.7
Compiled against GDAL/OGR   1.11.3  Running against GDAL/OGR    1.11.3
Compiled against GEOS   3.5.0-CAPI-1.9.0    Running against GEOS    3.5.0-CAPI-1.9.0 r4084
PostgreSQL Client Version   9.5.1   SpatiaLite Version  4.3.0a
QWT Version 6.1.2   PROJ.4 Version  492
QScintilla2 Version 2.9.1

Update 2 (SOLVED, sort of): My current version is 2.14.3-Essen

The solution somehow involves dealing with the projections, but I'm not sure what exactly it's going on. To make it short, I just started a new project in qgis, imported the landsat 8 .tiffs layers and used raster calculator to get the NDVI.

It's still intriguing what I got in my main project though. When I make the calculations, I get the correct result, but only in a small region (gray area)(note that on the sides, the region was cropped beforehand to make my layers smaller. This does not explain the southern horizontal limit though):


texts:
I have some big images without geographic information which I want to load them by leafletjs.

I have two ideas:

1 Georeferring images first

Once the image contains geographic informaiton,  I can serve them by geoserver and generate tiles.

2 Use pixel coordinate

Is it possible to generate tiles directly and load them in leafletjs using the pixel coordinates? If yes, what's the projection in this case?


texts: I've built a tile server with data for the whole country, as well as a single state, but I was wondering if it's possible to combine different regions (load multiple pbf files) into the same database?


texts:
I have successfully built a Web server using ubuntu 14.04, Qgis 2.6, with the plug-in lizmap as a local server.  I have several polygons being displayed and some WMS raster's from NOAA's WMS.  I know that this is possible using NOAA's service because it is stated here that "LayerInfo web service provides access to additional useful layer information" and the layer information query is provided (examples of the correct html format are on this link).

In my web browser I can Identify attribute data of the features that I have created (points, polygons, lines).

However I am unable to identify the raster attribute data.  All I want to be able to do is click on the raster and have the legend appear, it does not need to be a pixel specific legend the image provided by NOAA is sufficient (see below image example).

I was unable to get this to work on my WMS layer so as a test I made a custom pop up of this legend and assigned it to my own data and the legend appeared...Why can I not do this with the WMS raster layer and is it possible?  If not what is the work around?

http://nowcoast.noaa.gov/LayerInfo?layer=EOFS_RAS_SPD_NCST&data=legend http://nowcoast.noaa.gov/LayerInfo?layer=RAS_RIDGE_NEXRAD&data=legend


texts: Problems installing the Astun Technology Loader in Windows 7

texts:
I previously asked a question Using OS Mastermap SLDs in QGIS I got several answers suggesting that I use the Astun Technology Loader. There are details on how to do this https://github.com/AstunTechnology/Loader/wiki/Installation.

However this has proved very difficult in Windows 7

The first problem was the PATH environment in the OSGeo4w shell which is not the same as the standard 'cmd' shell. In the OSGeo4W shell having changed to the directory where distribute_setup.py is located the first problem is trying to run 'python distribute_setup.py'. It is necessary to either add C:\python33 to the PATH variable or directly reference it. This then runs OK.

From the same shell and directory and having located easy_install in C:\python33\scripts I have tried running 'easy_install lxml==2.3' this fails giving the errors
 ERROR: b"'xslt-config' is not recognised as an internal or external command make sure the development packages of libxm12 and libxslt are installed error: Unable to find vcvarsall.bat

There was no file called 'vcvarsall.bat' on my PC. From searches on the net I discovered it is part of Visual Studio 2008. Luckily I had a copy of Visual Studio 2008 which I have now installed. 'vcvarsall.bat' is now installed in the directory C:\Program Files (x86)\Microsoft Visual Studio 9.0\VC. I have added this directory to the PATH variable in the OSGeo4W shell without any success.

I am now completely stuck - Any suggestions?


texts: Precision vs Resolution in elevation data

texts:
Update:

A good example of something I am confused about: on the Wikipedia page for "digital terrain model" one can find the following statement about a satellite in the "See Also" section:

  TanDEM-X: generation of a world-wide, consistent, timely,
  high-precision Digital Elevation Model

Can I take this statement and swap the word precision for the word resolution? Would you agree that it means the same after doing so?

Original Question:

A question about Precision vs Accuracy has already been asked and answered, but I believe my question is different.

In the context of files that store elevation models, a TIFF for example, what is the difference between Precision and Resolution in the elevation detail? Are they the same thing?

Normally, when I think of low resolution elevation models I think of terrain that ends up looking blocky:

and high resolution corresponding to terrain that will look smooth and have more detail.

Can these terms (precision/resolution) be used interchangeably when talking about the "detail" in an elevation model? Is it erroneous to say that the terrain in the image above has lost precision from its original representation?

This post almost looked like what I was looking for, but I believe not - because it speaks of the terms in the context of a different domain - the instrument measurements domain. (Correct me if I am wrong here)


texts:
It seems that by default, vector layer labels only render after everything else is loaded.
Some wms/wmts feeds being very slow, this means you have to wait for a while before you can see your labels, even though they're local data and as light as text gets.
Is there a way to prioritize them so they display at the same time as vector points, lines and polygons or should I issue a feature request on github?


texts:
I'm new to HTML and to CARTO VL, hoping to jump right in to make a map with animated lines and points, which combined are unavailable in CARTO Builder. I've tried to run the basic CARTO VL set-up guide here: https://carto.com/developers/carto-vl/guides/getting-started/

But I can't even get this to work. I've pasted the exact code into an IDE, and then opened the .html file with a Google Chrome browser. The webpage is just blank. Is it a matter of subbing in my own API key? Or using a different IDE? Or opening the file a different way?

I think I'm missing something fundamental and obvious to programming - not sure what it is!


texts:
I have numerous shapefiles with two fields, "ID2" (Integer) and "FILENAME" (String). Using arcpy I'm trying to query the attribute table for rows that have a certain value ("SCN_TCV_001_2002") in FILENAME and for those rows find the max value of the ID2 field. I think I'm just missing something in the query syntax. My code is below which returns an error at the last line "AttributeError: 'NoneType' object has no attribute 'getValue'". I'm using ArcGIS 10.2.2 with an Advanced License.

import arcpy, os, sys
from arcpy import env

env.workspace = "C:\\folder1"
fclist = arcpy.ListFeatureClasses()

for fc in fclist:
    field1 = 'FILENAME'
    field2 = 'ID2'

    sWhereClause = '"' + field1 + '"' + ' = ' + "'" +  "SCN_TCV_001_2002" + "'"

    maxValue = arcpy.da.SearchCursor(fc, sWhereClause, "", "", field2 + " D").next().getValue(field2)


texts: Calculate Euclidean distance like ArcGIS does using PyQGIS

texts: I want to build a plugin like to calculate Euclidean distance like ArcGIS, or are there any existing plugins?


texts: I saw that QGIS now works on Android, any chance there will be an iOS version?


texts: I am trying to compare two DEM products ASTER and SRTM v.3, but the ASTER has the black background which i think will interfere with my analysis. I don't have access to ArcGIS where I can do this easily, I am only working with ENVI 5.2


texts:
I have wms layer and after i do cql filter on this. And now i want to see result. I want map zoomed to this result. My code.

ex=layer.getExtent();
map.zoomToExtent(ex);

But getExtent() returns wrong bounds. In Geoserver i have layer settings

Native SRS=70066
Declared SRS=900913
Reproject native to declared

In OL when create map i write projection: EPSG 900913;

What im doing wrong? I ll try do transformation

ex=layer.getExtent().transformation(map.projection, map.displayProjection);

projection=70066
displayProjection=900913

But it nothig to change.

Maybe im gonna do it whis Capabilities?


texts: Conversion to normalized radar cross section of PALSAR data in Google Earth Engine

texts:
I am mapping deforestation driven by palm oil expansion in Malaysia.

I am using ALOS PALSAR data from JAXA available directly on Google Earth Engine. Those data are raw and I'd like to process them to better detect the plantation.

I speckle filtered the data with a code already available. I still need to convert data to normalized radar cross section. Is there a code already available to do this normalization?


texts: Tooltip overlay with cartodb createlayer

texts:
I enabled a tooltip on my viz using the code below.

cartodb.createVis('map', vizUrl, {
    center_lat: 16.55,
    center_lon: -106.5,
    zoom: 6
})
    .done(function (vis, layers) {
        spillLayer = layers[1];
    layers[1].getSubLayer(1).setInteraction(true);
        layers[1].getSubLayer(1).setInteractivity('v1');

        vis.addOverlay({
            layer: spillLayer,
            type: 'tooltip',
            template: '
{{v1}} gallons
'
        });

})

I am now in a situation where I need to load the viz after a L.Map has already been created. How would I enable the tooltip in this situation? When I run the code below I get an error like

cartodb.createLayer(map, vizUrl)
     .addTo(map)
     .done(function (vis, layer) {

        spillLayer = layer[1];
        layers[1].getSubLayer(1).setInteraction(true);
        layers[1].getSubLayer(1).setInteractivity('v1');

        vis.addOverlay({
            layer: spillLayer,
            type: 'tooltip',
            template: '
{{v1}} gallons
'
        });

)}


texts: Sub-meter Bluetooth GPS units and ESRI Collector App

texts:
We have been using Trimble XT/XH for years now, but we are at a time when we need to upgrade and do not want to drop another $8000 for another Trimble that only causes us headaches and not very user friendly. We have begun using the ESRI collector App and Garmin Glo for some of our projects that don't require sub-meter and everybody seems happy with it. However, since we can't post-processes this data it won't be submissible for some of our projects that require documentation stating the data's accuracy level (e.g. wetland delineations).

From what I've gathered from the research of the newer sub-meter bluetooth units, they are capable of real-time differential correction, but the data can't be post-processed like Trimble Pathfinder/ESRI GPS Analyst. Are there any units that will do both real-time differential correction and post-processed differential correction?


texts: In QGIS using Python, how can I get multiple pixels from a raster based on a vector point file

texts:
I just started working with QGIS, and for the life of me I cannot figure out how to solve this problem:

I have several points of interest, these are stored in a shapefile (just simple points). I have (for now) one raster file with values. What I want to do basically is select those cells in my raster that are near the points in the shapefile (I want several window-sizes, like first only one cell, than 9 (3x3), than 25 (5x5) etc). I want to average these cells and get the median and stdev.

I want to do this in a Python script, since I will probably have to do this again (and again) at some point. I am new to QGIS, but in ArcGIS I would have tried converting my points to a raster, and 'grow' this raster while calculating zonal statistics.

However since in QGIS it should be easier to get data-level access, I was hoping to avoid all that and just convert my points to a position in a raster, and than loop over the values using Python. I already figured out how to get to the data using ReadAsArray:

raster = gdal.Open(inputRaster)
band = raster.GetRasterBand(1)
data = band.ReadAsArray()

However I am not sure how (if it is even possible once the data is converted to an array) to link the positions in my shapefile to the position in the array.

If this is not possible is there another strategy I should consider?


texts:
I am using QGIS 3.12 and QField 1.5.3 on a Samsung Galaxy S4 (Android 9).
My project allows collection of photos. This works when stored on Internal Storage, but not when stored in 'QField directory on external storage'.  The native camera app is displayed, you can take a picture, but no photo is created.
I am aware that this was a known issue some time ago, but thought it had been resolved.  I thought I had had it working on external storage, but not now.
Is it still the case that photo storage will only work on Internal Storage?
In app manager on the tablet, I have tried having QField set to internal and external storage, but it makes no difference. I have also cleared the cache and rebooted several times, to no avail.


texts: How do I make a field of summarized zip codes, 1 for each, put it into my data table and join the table to shapefile without losing data?

texts:
I'm trying to use a field in my shapefile with summarized zip code data. I have a table of hospital admissions data and in this table there is a field of zip codes which are repeated multiple times. What I'm trying to do is join this field of zip codes to my existing US zip codes boundaries shapefile, which mentions the zip codes only once. So the relationship is one:many.

I realize I must summarize the field on my table and then join it to the shapefile that way, but I don't want to lose all of the other data. In lieu of doing just a relationship between the shapefile and my data table, I'd like to do a join so that I can display the data in various ways. Is there a way to do this? Here is a picture of my tables:


texts: Is it safe to delete .qgd files in QGIS 3?

texts:
I am managing some of my QGIS 3 map files in a git repository. With QGIS 2, the map files (.qgs) are text files, which can be easily handled by git; QGIS3 introduces an additional binary file (.qgd) which presumably stores some auxiliary information. A compressed binary format .qgz is also provided with QGIS 3. These binary formats are not text line based, and difficult to manage using text oriented tools such as diff and git. Therefore, I do not want to track them in a repository unless I have to.

How critical is the .qgd file?

Can I delete them or omit them when I use git to commit a map without seriously messing up the map?

(I'm OK with loss of auxiliary information. Just need basic layer information and basic symbology to be committed, i.e. the choices of symbology and logical expressions I entered when defining a layer)


texts: Adding coordinate system to .dwg in order to convert to .shp

texts:
I have a number .DWG files that appear to have no georeferencing or coordinate system attached. They were probably created in line with OS master mapping.

My aim is to try and convert the lines in these drawings to .shp in order to add to an existing project within QGIS. To do this I am assuming that I need to apply a coordinate system as well as reference a point within the dwg or shape to that reference system. (As you can tell, I am guessing)

I am using NanoCad to open the .dwg

Is there a way to do what I want to achieve, and if there is, what is the best way to  apply a coordinate system (In NanoCad, or convert and then in QGIS).


texts: How to programmatically click on Feature or particular region without coordinates?

texts:
How  I can programmatically click on feature using javascript library, lets say in world countries map how can I programmatically click on a particular country ( i don't want to use co-ordinate as my does not deal with them).

Country Name comes form rest services. I will highlight the country through querytask i pass countryname in the country i get the feature object then i add symbol object to it and highlight it.

var queryTask = new esri.tasks.QueryTask("http://10.32.9.106:6080/arcgis/rest/services/KX_Mapping/World_countries_less_Antartica/MapServer/0");
var query = new esri.tasks.Query();
query.where="CNTRY_NAME IN ('Some')";
    queryTask.execute(query);
dojo.connect(queryTask, "onComplete", function(featureSet) {
                 try{
                        console.log(" query complete");
                        var symbol0 = new esri.symbol.SimpleFillSymbol().setColor(new dojo.Color([255, 0, 0, 0.9]));
                        dojo.forEach(featureSet.features,function(feature){
                                //alert('feature'+feature);
                                var graphic = feature;
                                graphic.setSymbol(symbol0);
                                //alert("infoTemplate"+infoTemplate);
                                var infoTemplate = new esri.InfoTemplate("Alert", "ds");
                                graphic.setInfoTemplate(infoTemplate);
                                //alert("graphic"+graphic);
                                //alert("map.graphics"+map.graphics);
                                map.graphics.add(graphic);
                        });
                    }catch(e){
                        alert("e"+e);
                    }
                });

I have associated infotemplate object to graphic object , template will be shown when user click on that country , how can i programitically  make it clickable.


texts: Delete / remove just one feature from a topojson

texts:
I've got a TOPOJSON file that is nearly perfect, but contains one item (a multi polygon) that I'd like to remove. Doing so manually from the code is difficult because it'd require manually counting through hundreds of arcs in the nested arrays.

The feature in question has an ID I can identify it by, so I thought this would be easy using the MapShaper console. However looking in the docs I'm struggling to find any tool that seems to do this. Closest is -erase:

  -erase. Remove features or portions of features that fall inside an area.

...but it seems like all the direct targeting options are for targeting layers, not features, or by drawing a bounding box (which for this geography is tricky). Everything I've tried to select the feature by ID instead tries and fails to find a layer with an error message Missing data layer.

Surely there's a simple way to do something like (fake code):

-delete [id='someId']

I'm open to suggestions using (ideally open-source) tools other than MapShaper, but I'd prefer to not need to convert my TopoJSON to another format then back again, and would prefer a simple solution to such a seemingly simple problem.


texts: How to do a buffer operation for a polygon drawn by the user using python in QGIS

texts:
I have a polygon drawn on the map canvas.Say the polygon would be drawn by using "Measure Area" tool in qgis.Now, I want to do buffer for certain distance around that polygon.I can see the below operation in QGIS in which it uses a layer as input parameter instead i want the polygon as an input in the buffer.How should this be accomplished in QGIS using python.

processing.runalg("qgis:fixeddistancebuffer","layer-name-you-need",0.1,10,False,"Path-to-file-to-store")


texts: Quickly create contour polygons from cost points using OSS for web app

texts:
This is related to a previous question I've asked:
Create filled contour map with GRASS

I've created a point layer using OpenStreetMap data and pgRouting's driving_distance() function, where the points contain a cost attribute which contains the walking distance from a school. I've used QGIS and the contour plugin to create what I want:

However I'd now like to do this programmatically for a web app, and in addition result in polygons rather than a raster. I've got something working using GRASS, however it is quite slow, around 30 secs on my machine versus less than 1 sec using QGIS and the contour plugin. Here the short version (or see Python version):

# read file and set resolution and bounds
v.in.ogr dsn=cost_points.shp output=cost_points
g.region vect=cost_points rows=256 cols=256
# interpolate points and make bands of 0.4km increments
v.surf.rst input=cost_points elev=interp zcolumn=cost
r.mapcalc res='int(interp/0.4)'
# convert to polygons and write out file
r.to.vect -s input=result output=polys feature=area
v.out.ogr -c input=polys dsn=polys.shp type=area

Are there other GRASS functions I should use? r.surf.rst takes the majority of the time, so perhaps it is too robust for my needs, but even then the results aren't smooth enough for my liking. Are there other open source libraries I should be considering for this?

If you'd like to know more about my overall goals, see this page.


texts:
I've just started to use the Time Manage plugin and I was wondering if how compatible it is with PyQGIS. Specifically I was wondering if it was possible to use PyQGis to set up the following features of Time Manage:

Adding layers
Setting start time
Frame length
Setting unit


texts: How to re-calculate area, perim, p2a,width and length in QGIS?

texts:
I cut a part out of one feature and added it to another feature in a polygon shapefile.
My problem now is that both features has the same values in the area-, perim-, p2a-, width- and length-columns?
The features are very much different in size, how do I get QGIS to sort of "re-calculate" the areas etc?


texts: How to obtain statistics within a Buffer/Catchment in MapInfo?

texts:
I need to calculate the population statistics within a 10 minute drivetime catchment from a point and within 1 Km Buffer from a point in MapInfo, how can I do that can someone please elaborate? I have a population thematic map and I have made a buffer on one location and I need to find the total population living within that area? I am also using a Drivetime addon with MapInfo for drive time catchments and Mapinfo to make Buffers.
Many thanks


texts: Determining whether the geometry field in a SQLite table has z coordinates using Python

texts:
In the Shapely package, there is a function that allows you to determine whether coordinates include a Z value. It is object.has_z.

Is there a way that a similar command can be carried on the geometry column of a SQLite table?


texts: comparison operators with Process Modeler's Raster Calculator

texts:
I am trying to process some raster imagery using QGIS 2.18.7 on OS X 10.10.5.
I want to use a comparison operator (<, >, etc...) within the raster calculator in the process modeler (the one located under GDAL[Miscellaneous], not the SAGA one).

I can successfully get the calculator to perform an arithmetic operator by entering "A+1" in the "Calculation in gdalnumeric syntax using +-/* or any numpy array functions (i.e. logical_and())" field (quotes mine). Like I said, this works but produces an error/warning I don't understand.

  Algorithm Recoding1Month starting... Prepare algorithm:
  GDALOGRRASTERCALCULATOR_1 Running Raster calculator 1 [1/1]
  Parameters: INPUT_A
  =/QGISWork/MODFireBurned/Tiffs/2001/y2001NicaBurnedUTM.tif, BAND_A =1, INPUT_B =None, BAND_B =1, INPUT_C =None, BAND_C =1, INPUT_D =None,
  BAND_D =1, INPUT_E =None, BAND_E =1, INPUT_F =None, BAND_F =1, FORMULA
  =(A+1), NO_DATA =, RTYPE =1, EXTRA = GDAL command: gdal_calc.py --calc "(A+1)" --format GTiff --type Int16 -A
  /QGISWork/MODFireBurned/Tiffs/2001/y2001NicaBurnedUTM.tif --A_band 1
  --outfile /var/folders/8y/6924gdvx3jg6k2rz5yrzwz100000gr/T/processing1524180f5e1443d3a74b4fdc0a8bd484/27a58e4ed8ce4200bb71345c1a315e52/OUTPUTALGGDALOGRRASTERCALCULATOR1.tif
  GDAL command output:
  /Library/Frameworks/GDAL.framework/Versions/2.1/Programs/gdal_calc.py:272:
  FutureWarning: comparison to None will result in an elementwise
  object comparison in the future.  myNDVs=1*numpy.logical_or(myNDVs==1,
  myval==myNDV[i])  0 .. 10 .. 20 .. 30 .. 40 .. 50 .. 60 .. 70 .. 80 ..
  90 .. 100 - Done  Converting outputs OK. Execution took 0.352 ms (1
  outputs). Model processed ok. Executed 1 algorithms total Converting
  outputs Loading resulting layers Algorithm Recoding1Month finished

More importantly, I cannot seem to get a comparison operator to produce the expected boolean (1 or 0) output when I enter "(A < 1)".  In fact this crashes and fails to produce any output file. This yields the following log message:

  Algorithm Recoding1Month starting... Prepare algorithm:
  GDALOGRRASTERCALCULATOR_1 Running   [1/1] Parameters: INPUT_A
  =/QGISWork/MODFireBurned/Tiffs/2001/y2001NicaBurnedUTM.tif, BAND_A =1, INPUT_B =None, BAND_B =1, INPUT_C =None, BAND_C =1, INPUT_D =None,
  BAND_D =1, INPUT_E =None, BAND_E =1, INPUT_F =None, BAND_F =1, FORMULA
  =(A=, RTYPE =1, EXTRA = GDAL command: gdal_calc.py --calc "(A GDAL command output:
  /Library/Frameworks/GDAL.framework/Versions/2.1/Programs/gdal_calc.py:272:
  FutureWarning: comparison to None will result in an elementwise
  object comparison in the future.  myNDVs=1*numpy.logical_or(myNDVs==1,
  myval==myNDV[i])  0 .. 10 .. 20 .. 30 .. 40 .. 50 .. 60 .. 70 .. 80 ..
  90 .. 100 - Done  Converting outputs OK. Execution took 0.381 ms (1
  outputs). Model processed ok. Executed 1 algorithms total Converting
  outputs Loading resulting layers Algorithm Recoding1Month finished

Note that, within these output logs, the command seems to be interpreted the way I'd expect for the arithmetic operator "EXTRA = GDAL command: gdal_calc.py --calc "(A+1)""  but is messed up for the comparison operator: "EXTRA = GDAL command: gdal_calc.py --calc "(A " (my external quotes).  So clearly I'm making some sort of error but I cannot find helpful examples of how to format for "gdalnumeric syntax"  (a Google search for this yields only 115 results) or for "numpy array functions" within the raster calculator.  I found this more complex example which seems to suggest the comparison operators are possible: gdal_calc raster calculator syntax for logical operators and other functions. In fact, this (http://www.gdal.org/gdal_calc.html) seems to suggest that (A<1) should be allowed within the numpy syntax.  I have tried adding quotes, but it still doesn't work. So, I'm very confused.

How can I use comparison operators with the process modeler raster calculator?


texts: Output formats for QGIS Heatmap plugin not appearing in drop-drown in Ubuntu

texts:
I've been trying to use the heatmap plugin on QGIS (2.8.2) on Ubuntu(14.04) but the drop down for the output formats doesn't get populated.

Selecting the drop-down allows me to press OK on the plugin popup which gives "GDAL driver error: Cannot open the driver for the specified format".

I've installed GDAL according to this. The drop down gets populated and works as expected in Windows. How can get the same to work on Ubuntu?

Update : I've tested on a different Ubuntu (14.04) machine and it worked as expected, so the problem seems to be occurring on only some machines


texts: Adding, Updating, and Deleting information from arcGIS Server

texts:
I am new to ESRI and GIS. I am starting to get into using client data. I have created maps and put them to the server using arcMap and arcGIS for Server. I need to know how easily I can access these fields and if it is even possible to change them since the mapping data will change from time to time.

Attached is a screen shot of the fields in my map layer. I need to do this somewhat dynamically since each customer is different.

These are the fields for my layer (WVALVE):

I would like to be able to change what is contained in each field through the application and even take the information given in the maps and put it in my own database.


texts: Editing the DBF for a shapefile is yielding totally different results in R and QGIS, but the attribute tables are identical

texts:
I'm having kind of a peculiar issue that I'm hoping someone can shed some light on. I'm actually creating my map in R, but I'm using QGIS as a tool to check my mapping process (R newb here). So let's dig into it.

I'm trying to plot rivers in my area of interest (Tanzania), but the source shapefile (diva-gis.org) includes tons of unnecessary water lines. For the sake of reproducability, I'll be including all edited DBFs and shapefiles. Here is the original DBF file and here is the original shapefile. It looks like this.

As you can see in the DBF, many of the rivers are unknown/unnamed and listed as "UNK" which isn't very useful for me. In QGIS, I am simply able to open the attribute table, sort by name, and remove all UNK values. Alternatively, I can filter out those values through the SQL expression:

"NAM" != "UNK"

Both of these yield the same results. Here is the QGIS edited DBF and here is the QGIS edited shapefile. This is what the map looks like:

Then things get quirky. Again, I'm actually creating this map in R. Using the original source data, I edit the DBF with the "foreign" package, overwrite the DBF, and then map using ggplot. I'm getting a totally different result than the one pictured above. Here's a little snippet of example code (not for the full map, but this is essentially what this portion of the map code will look like):

library(rgdal)
library(foreign)
library(ggplot2)

#set your appropriate working directory
setwd("D:/Mapping-R/Returns-Practice")

#read dbf, sql expression to remove all UNK values in NAM
water.lines <- read.dbf("original-waterlines.dbf")
water.lines <- water.lines[water.lines$NAM != "UNK",]

#overwrite dbf
write.dbf(water.lines, "TZA_water_lines_dcw.dbf")

#read shapefile, fortify to create readable data frame
water.lines.shp <- readOGR(dsn = "D:/Mapping-R/Returns-Practice",
                  layer = "original-waterlines")
water.lines.shp <- fortify(water.lines.shp)

#plot
ggplot() + geom_polygon(data = water.lines.shp,
             aes(long, lat, group = group),
             color = "slategray4")

The result of this code looks like this:

Kind of stretched out and funky, but clearly does not follow the same line pattern as the QGIS-edited image. This has been driving me absolutely nuts. I actually went through and checked the R attributes against the QGIS attributes and they matched up perfectly, but the image result is completely different. After playing around with it some more, I figured that my code was to blame and that's where I messed up. But then I opened the original DBF and manually removed all UNK values and plugged it back into QGIS. Here is the resulting DBF and here is its corresponding shapefile, which matches the R image/plot and looks like this:

What gives?! The attribute tables are identical but the images are different. I feel like there has to be something really obvious that I'm missing but I've been playing with this for literally hours and can't figure it out.


texts: How to fill a field with the ranking of another field for a shapefile attribute table in ArcGIS Pro?

texts:
I want to add a field ("FreqRank") in my attribute table with the rank of the field "COUNT", so it was sorted descendingly (see below).

I tried the code here from another post below (https://community.esri.com/thread/213033-code-field-based-on-rank-of-another-field), but it just give me the OBJECTID.
#Codeblock:

rank = 0
lastValue = None

def getNextRank(value):
  global rank
  global lastValue
  if value == lastValue:
    return rank
  else:
    rank += 1
    lastValue =value
    return rank


texts: Defining next object in street with android?

texts:
I want to put objects in a map and discover which is the closest object from me, supposing that I'm in a moving car with my cell phone. Thus, to define the closest object I need to take the orientation of the movement and the street into account. The following figure illustrate the situation:

So, object 2 is the closest object, but it is not in the same street as I'm (the car). The object 1 is the closest object in the street that I'm , but I'm going to the opposite side. Finnaly, the object 3 is the one I want, since it is the closest object in the same street I'm and also following the orientation of my movement.

Does anybody know how I can obtain the closest object to me (considering that I'm moving) according to the rule I just mentioned?

If you don't know an API, but know what kind of concepts and stuff I need to study would be helpful, because I'm new in android and geographic systems development.


texts: What's the best way to break closed polyline using ArcObjects?

texts:
I have a polyline feature class with several closed lines (the start point = the end point). I need to break these lines in two features, using ArcObjects (doesn't matter Java, VB or C#).

I have already detected these closed lines in this way:

IFeatureCursor featureCursor = lineFeatureClass.search(null, false);
        IFeature feature = null;
        while((feature = featureCursor.nextFeature()) != null) {
            IGeometry geometry = feature.getShape();
            IPolyline polyline = (IPolyline)geometry;
            if(polyline.isClosed()) {
                log.info("closed line detected");
                boolean[] splitHappened = new boolean[1];
                int[] newPartIndex = new int[1];
                int[] newSegmentIndex = new int[1];
                polyline.splitAtDistance(0.5, true, true, splitHappened, newPartIndex, newSegmentIndex);

            }
        }

Now, I get a polyline composed of 2 parts (or path), and I don't know how to get them.

Then, I need to remove the second path from the current feature, and create a new feature with the second path as a geometry and the same attribute values.

Does anyone know what's the best way to do that?


texts: How to correct topology on a shared boundary by removing nodes from a feature?

texts:
Maybe you can help me solve a problem regarding topology cleanup that I recently came a cross.

I have two layers with bunch of polygons (or feature classes as I'm in the ESRI world) that should share a common boundary. (ex. water body polygon and vegetation coverage - the shore line is the common boundary).

Usually Integrate in ArcMap would be the tool for the job, but in this case I'm only allowed to modify one of the layers, the other one is the reference layer and it must not be touched. In other words I must snap polygons in one layer to the polygons of the other at the common boundary.

Here's and illustration of the problem:

The defective border has some extra nodes (marked in blue). If I use Integrate those blue nodes will be introduced to the reference polygon layer regardless of ranks set. So in this case I need a tool that will make the boundary identical by removing nodes from the polygons of the defective layer.

As a side note. The rest of the nodes on the shared boundary are snapped together, just the blue ones mess my border.

Any advice?

I have ArcInfo 10.0 and ET Geowizard available and probably can work with OSGIS soft, but would prefer to avoid converting my layers to shp.


texts: We need to have a login occur before QGIS3 loads

texts: Where is the hook for the invocation of the QGIS3 launch? I have looked at the pyqgis api and did not find a bootstrap class or anything similar.