texts:
I'm trying to figure out an error in GRASS network analysis for several days now. I started to suspect that it has something to do with my OS, namely Windows 7. Hence I installed Linux-Mint 14 (Nadia) on my computer. I would like to run GRASS 6.4 (or latest) on Linux to test the function, yet can't figure out how to do so.

Can someone help me with it? Or direct me to some where on the web where it is explained?


texts: How to calculate how many grid squares a line shapefile covers

texts:
I have a load of transect data (in the form of a line shapefile) covering a very large area (approx. 11,900km2) that has been split into 20m x 20m grid squares on QGIS using:

"Vector -> Research Tools -> Create Grid..."

I would like to know how to calculate (without having to physically count) how many of these grid squares my line shapefile covers.

I can only really do simple things on QGIS and don't really use it extensively, so I'm seeking a step-by-step description.


texts: I have the following map but I don't know how to make a title appear on top of it. Any suggestion on how to?


texts: OSgeo4W64 Windows 7 & ImportError: No Module Named Site

texts:
I realize this question has been asked a lot, but I was not able to resolve my issue through the answers posted elsewhere... I am new to Python, gdal, OSgeo4W, etc. but have experience with ArcGIS.

I'm receiving the dreaded "ImportError: No module named site" error from Windows Cmd after running the following batch file by clicking on it directly from windows explorer (it is save as a .bat file).

  ECHO Off
  C:\OSGeo4W64\bin\gdalbuildvrt.exe depth_max.vrt *.tif
  python C:\OSGeo4W64\bin\gdal2tiles.py -s EPSG:2278 -z 12 -w none "depth_max.vrt" depth_max
  pause

Couple observations/comments:

The file runs fine when running it from the OSgeo4W window
gdalbuildvrt.exe executes fine in either scenario, the error is associated with the python command
Python 2.7 is installed at C:\Python27\ArcGIS10.2\ and C:\OSGeo4W64\apps\Python27
Python 3.6 is installed at C:\OSGeo4W64\apps\Python36
Running Windows 7 (Note: This may be run on Windows 10 as well)
gdalbuildvrt does not have a .py file; only .exe, which I found odd but I guess it is not a python script but a standalone program?
I tried including the SET commands for PYTHONHOME and PYTHONPATH in the batch file and it either did not work or broke python (it would then say DLLs were not found)

My eventual goal is to run this from a SHELL command via another program.


texts: Contours generated in Leaflet losing coordinates?

texts:
I am attempting to use a JavaScript library to generate contours in a leaflet map from a point layer on the map.  I have tried this approach with d3 and turf and am getting the same issue where the output contours (either polygons or polylines) have the appropriate Z value, however the coordinates are not populated:

This is my code for using the d3 library:

m.eachLayer(function(layer){
    if (typeof (layer.options.name) != "undefined") {
        if (layer.options.name == $('#ContourLayer').val()) {
            arr = layer.toGeoJSON().features

            console.log(collection)
            bounds = L.geoJson(layer.toGeoJSON()).getBounds()
        }
    }
});

console.log(bounds)

var width = m.latLngToLayerPoint(bounds.getNorthEast()).x-
m.latLngToLayerPoint(bounds.getSouthWest()).x  //Math.abs(Math.abs(bbox[0])-
Math.abs(bbox[2]))
var height = m.latLngToLayerPoint(bounds.getNorthEast()).y-
m.latLngToLayerPoint(bounds.getSouthWest()).y //Math.abs(Math.abs(bbox[1])-
Math.abs(bbox[3]))
console.log([Math.abs(width), Math.abs(height)])

var contour = d3.contours()
.size([Math.abs(width), Math.abs(height)]);
var contoursGeoData;

contour.thresholds(breaks);

contoursGeoData = contour(arr);

console.log(contoursGeoData)

I did the same process with the isolines function in the TurfJS library and got the same output (just polylines instead of polygons).  What am I doing wrong?

Also, my data (and map) are all in WGS84 (EPSG:4326).  Does the data and map have to be projected for the contours to draw?


texts: Calculating field in ArcGIS Desktop using ArcPy cursor?

texts:
I'm trying to calculate a field based on the values of others in the same layer in a calculation. I can't work out why its not working. I've even cut everything down to just be based on one field and even that won't work, so frustrating!

This is based on one field

rows = arcpy.UpdateCursor(inputFC)
for row in rows:
    val = row.getValue(Infield1)
    row.setValue(outfld , val)
    rows.updateRow(row)
    del rows, row
    arcpy.AddMessage("Completed")

this is the final script. Its meant to calculate calories based on minutes, speed and weight (weight as a number you input rather than a field)
import sys, arcpy, math

inputFC = arcpy.GetParameterAsText(0)      # Input feature class
timeField = arcpy.GetParameterAsText(1)      # Input time in minutes
speedField = arcpy.GetParameterAsText(2)      # Input speed
SlopeField = arcpy.GetParameterAsText(3)      # Input Slope
weight = float(arcpy.GetParameterAsText(4))      #input weight in kg

arcpy.AddField_management(inputFC, 'Calories', 'LONG') #create calorie field
rows = arcpy.UpdateCursor(inputFC)
for row in rows:
    time = row.getValue(timeField)
    speed = row.getValue(speedField)
    slope = row.getValue(SlopeField)
    calories = time * (speed * 3.5 * weight)/200 * slope
    row.setValue('Calories', calories)
    rows.updateRow(row)
    del row, rows

this says completed but the calories field ends up being null


texts: Troubles exporting raster to tiff in ArcGIS

texts:
I am exporting in ArcGIS a raster by doing Data > Export data. But in the Export Raster Data window, I don't have the option for choosing the file format.

Am I missing something? Perhaps must the Coordinate reference system of the raster be a projected one?


texts: Creating 1 hectare rectangle polygon surrounding points and dividing it into equal sized sub-quadrats in ArcGIS Desktop?

texts:
I have a set of points like this:

]

I have only these points' information. Now I want to draw a boundary of this points whose size should be 1 hectare (100m x 100m). After that, I want to subdivide this area into 16 equal quadrats like this: Each side of subplots should be (25m x 25m).

How can I do this in ArcGIS Desktop?


texts: Azimuth (bearing) changes depending on map/project projection

texts:
I'm trying to calculate the initial bearing that is needed for the great circle between 2 points. As an example, let's say from London to New York. There are multiple calculators online and they calculate this as 288.3658° (http://www.movable-type.co.uk/scripts/latlong.html?from=48.609890,-122.426991&to=48.592836,-122.421613).

However, when I try to get this figure I get different results depending on the map projection. For this I'm using the ne_10m_admin_0_countries, and the ne_10m_populated_places shapefiles. Along with a new layer (lets call it LDNNY) with a CRS of EPSG:4326. Finally, I've got 'on the fly' set to TRUE with a project CRS of EPSG:4326.

Step 1. Draw a straight line on LDNNY from London to New York. Using the CAD-tools's Azimuth tool I achieve an initial Azimuth of 261.7193°.

Step 2. Densify this line to achieve a great circle arc (I do this with the Shape Tools plugin). Now I measure the initial part of this line to get the initial bearing (if I measure from London to New York I wouldn't be getting the initial bearing). This time I get 281.5755°.

Step 3. Change the project CRS to EPSG:3857 and I get 288.2342°.
Change the project CRS to EPSG:54030 and I get 284.6549°.
Change the project CRS to EPSG:54002 and I get 292.4153°.

These are just some example - you can see that some are very close to the online calculator but...
Please could someone explain:

Why this is.
How would I know what is 'right'?

N.B. If needed I could add images to help explain better.


texts: Cost Path based on Lowest Maximum Elevation Encountered

texts:
I am trying to locate the lowest elevation point at which a river will overtop its banks, and spill into lower-lying terrain behind them. I'm working in ArcMap version 10.2, and using a high-resolution DEM to represent the surface.

My thought was that I could approach this as a kind of cost-path problem, where the path would only be concerned with finding the smallest maximum elevation along a route between the river and low-lying terrain (like trying to find the lowest pass through a mountain range).

Is it possible to create a least-cost path that only accounts for the maximum value encountered? And if not, are there any other tools I could use to locate the spill point?


texts: Update the geometry of a point through WFS-T

texts:
I am trying to change a point geometry in a postgis table through geoserver.
I use python Requests in order to make the POST request.

My xml looks like this:

xmlstr = """<wfs:Transaction service="WFS" version="1.1.0"
             xmlns:wfs="http://www.opengis.net/wfs"
             xmlns:gml="http://www.opengis.net/gml"
             xmlns:ogc="http://www.opengis.net/ogc">
            <wfs:Update>
             <wfs:Property>
              <wfs:Name>the_geom</wfs:Name>
              <wfs:Value>
               <Point srsName="http://www.opengis.net/gml/srs/epsg.xml#4326">
                <coordinates>-155.7421875 -11.178401873711785</coordinates>
               </Point>
              </wfs:Value>
             </wfs:Property>
            <ogc:Filter>
           <ogc:FeatureId fid="wld_nhr_adameqepic14days_wfp_test_1_1_1.2"/>
          </ogc:Filter>
         </wfs:Update>
        </wfs:Transaction>"""

When I send the request as:

requests.post(url, data=xmlstr, headers=headers, auth=('user', 'pass')).text

I get back the following error:

<ows:ExceptionReport xmlns:xs="http://www.w3.org/2001/XMLSchema" xmlns:ows="http://www.opengis.net/ows" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" version="1.0.0" xsi:schemaLocation="http://www.opengis.net/ows http://localhost:8080/geoserver/schemas/ows/1.0.0/owsExceptionReport.xsd">

java.lang.NullPointerException

I also tried to run the above through the Demos page of geoserver. I get also the same error: java.lang.NullPointerException

Ofcource none of the points are updated.
Any idea what is going on? Do I have a mistake on my XML structure?
I have also checked the WFS OGC documentation and it seems like my XML structure is correct.

EDIT

I have also executed the DesribeFeature for this layer. I see that the_geom attribute is there as:

<xsd:extension base="gml:AbstractFeatureType">
    <xsd:sequence>
      <xsd:element maxOccurs="1" minOccurs="0" name="the_geom" nillable="true" type="gml:PointPropertyType"/>
      <xsd:element maxOccurs="1" minOccurs="0" name="fid" nillable="true" type="xsd:int"/>
      <xsd:element maxOccurs="1" minOccurs="0" name="time" nillable="true" type="xsd:string"/>
      <xsd:element maxOccurs="1" minOccurs="0" name="latitude" nillable="true" type="xsd:double"/>
      <xsd:element maxOccurs="1" minOccurs="0" name="longi


texts: Does PostGIS ConcaveHull exclude points/polygons in some cases, and buffer outward in others?

texts:
In trying to generate concave bounding polygons using PostGIS ConcaveHull, I am seeing what I think is strange behavior.

When running with a 0.99 tolerance, the resulting polygon shrinks tight as would be expected, but it seems to exclude some points/polygons in the input geometry set.  I would not expect there to be polygons fully excluded from the resulting hull.

Subsequently, when dropping the tolerance below 0.99, even to just 0.95 or below, the resulting hull includes all points or polygons, but it seems to kind of buffer out - as in the proximity of the hull boundary to many of the input geometries increases dramatically (for example by as much as 500 meters).  That extra area breaks my particular use case.

In this image, the points are the input geometries, and the lines represent the concave hull outputs.  The lower line is the 0.99 tolerance line, the upper line is a 0.60 tolerance line.  In other instances, the 0.99 is well outside of the 0.60, as it is in theory much closer to a convex hull.  Here however is an area where some concavity is being enforced, albeit not as expected.

If the 0.99 tolerance included all input geoms, I could use it to then erase the additional area 'buffered out' by the lower tolerance runs possibly.


texts:
I have a geopandas dataframe in EPSG:4326 and I'd make a new dataframe consisting of all the rows that fall within a certain bounding box.

First I get the bounding box that I care about (which is actually the bounding box of another dataframe):

print df_sussex.total_bounds
[ -1.57239292  50.57467674   0.14528384  51.27465152]

Then I make a dataframe consisting only of that bounding box:

pts = gpd.GeoDataFrame(df_sussex.total_bounds)

And finally I try to get all the features that intersect with that bounding box:

sac_sussex = gpd.overlay(pts, df_sac, how='intersection')

But this gives me AttributeError: No geometry data set yet (expected in column 'geometry'.

What am I doing wrong?


texts:
First of all I am pretty new with Leaflet, so I´m sorry if I have obvious fails.
My problem is that the layers loaded by leaflet looks stretched from the top and narrower from the sides.
The first pic is how should look and second how it looks.

The layers were created with QGIS in EPSG:4326. I read default CRS of QGIS is the same.
The code how load the layer is:

var singleRoom = this.roomsLayer.features
for (let i = 0; i < singleRoom.length; i++){
            let room = L.geoJSON(singleRoom[i].geometry,{className: "room-"+singleRoom[i].properties})
            this.roomsLayerGroup.addLayer(room)
     }
this.roomsLayerGroup.addTo(this.map)


texts: how to create the replica from a specific version?

texts: I have a default version SDE and I want to make a QC versioning in 2 levels then to create a checkout/checkin Replica from the child version is that possible ?


texts:
I crate a file that represents 16000 edges with information of start node, end node and cost in a table.

In MATLAB i want to solve shortest path with PSO algorithm. How i can solve it?


texts: ERROR 999999 from Project Raster - any advice?

texts:
I have a DEM layer (TIFF raster) for the island of Crete downloaded from the EU Copernicus EU-DEM website and a point data layer for locations of interest around the island.

I cannot seem to get the two layers to match up, whereby the point data lay on top of the DEM. Been searching through previous questions on here and seemingly nothing I have done is able to sort the issue.

Current data frame Coordinate System: *GCS_WGS_1984*
TIF Raster XY Coordinate System: *GCS_WGS_1984*
Point data Coordinate System: *GCS_WGS_1984*
Point data Projected Coordinate System: *WGS_1984_UTM_ZONE_34N*

I have tried to project the raster into WGS_1984_UTM_Zone_34N with the hope that they will match up. However, the Project Raster tool is not working and I keep getting the following error (see image). I have been trying for two days!


texts:
I have working PostgreSQL 9.5.3 with PostGIS.
I have sucessfully created table and column with geom data type.
I am able to execute geometry queries, namley ST_Transform without problem.

# SELECT ST_Transform(location_geom, 26986) FROM table;
                    st_transform
----------------------------------------------------
 01010000206A6900003CB104EFAB3D6241733751CA4EDF5941
(1 row)

However, when trying to use ST_ClusterDBSCAN it says, it does not know function like this.

SELECT ST_ClusterDBSCAN(location_geom, eps := 0,5, minpoints := 2) FROM table;
ERROR:  function st_clusterdbscan(geometry, eps => numeric, minpoints => integer) does not exist

Same result when I append "OVER ()" just after ST_ClusterDBSCAN(...), although I do not know what it does, but I have seen it in documentation.

http://postgis.net/docs/manual-dev/ST_ClusterDBSCAN.html

What needs to be installed, to use ST_ClusterDBSCAN?


texts: Subtracting polygon from another in same layer using QGIS?

texts:
I need to subtract a polygon from another bigger one, in the same layer.

I have a layer, with a polygon. Inside this polygon I need to edit another polygon, smaller, without overlap.

How can I do this?


texts: NA values in my raster after changing the resolution, extent, and origin

texts:
I'm working with two rasters that differ in their origin, extent, and resolution. I have a bathymetry raster, with a very high resolution (x=0.0008333333, y=0.0008333333) and a MUCH great spatial extent. I also have a sea surface temperature raster, which has a much coarser resolution (x=0.04166667, y=0.04166667). Both rasters have the same projection (longlat, datum=WGS84). Both rasters also have 'NA' values that correspond to landmasses.

Bathymetry data:

https://www.ngdc.noaa.gov/dem/squareCellGrid/download/4956

Sst data: June 25-July 2 (2016) 8 day composite, cropped to -131.8333, -127.6667, 50.125, 54.54166  (xmin, xmax, ymin, ymax)

https://oceancolor.gsfc.nasa.gov/cgi/l3?per=8D&prd=NSST_sst.nc&sen=A&res=9km&num=24&ctg=Standard&date=25Jun2016

I have manipulated the bathymetry raster (using R) to match the extent, origin, and resolution of the sea surface temperature raster.

bathy2<-projectRaster(bathymetry, sst, method="bilinear")

I made a data frame so that I could compare the sea surface temperature values and the bathymetry values at each pixel. I then noticed that I have a problem. Some of the pixels have 'NA' values for bathymetry, but numeric values for SST. How can this be, given that the landmasses should be in the exact same place? The original bathymetry raster has a much finer resolution than the SST raster, so its outline of the landmasses should be accurate.

Is this a result of using projectRaster(), or an issue with the data itself? Anyway I can fix it?


texts: Loading .qgs and .qpt files in standalone PyQGIS script?

texts:
I have a .qgs and .qpt files. And I would like to know if it is possible to use them in a standalone PyQGIS script?

I am using QGIS 2.18 and Python 2.7

This is the template that I would like to use.

And put this image inside that template using python if that is possible. I'm nearly finished with my thesis and this is the only part that is holding me back.


texts:
I mostly followed Creating animated route map in QGIS? for animating a route on QGIS. The main thing I did differently is that I started with a csv with timestamped lat/lon data, converted it to a GPX file using GPS Visualizer and skipped the step for adding timestamps but my animation only shows the current position, not the track.

How do I show the track?


texts: How to Copy single attribute into overlapping layer?

texts:
I have two polygon shapefiles, one of mine sites and one of parks. I've clipped one layer so that only mines within the park boundaries are shown. However, I want to take an attribute from the parks layer (park name) and add it into the mine layer so that when I export the table to excel it will show which park the mine is in. I can't merge the attribute tables because I need the park name to line up with the correct mine.

How do I copy the attributes from the Park layer to the Mine layer based on location?


texts: How to flatten spatially equal polygons in PostGIS?

texts:
I have the following dataset of polygons. All the pieces are broken based on their borders. I added some transparency to better understand what's going on. The pieces that have stronger colors, means that there're more "spatially equal"  polygons in that area (Note by spatially equal I mean ST_Within(A,B) == true and ST_Within(B,A) == true). Stronger the color is, the greater the number of polygons in this area.

(Here is my sample data)

I'm using the following code to flatten all the polygons (in other words, I want them all in the same transparency). Although, it works for some of the polygons, it does not flat all my polygons. Am I doing something wrong here?

SELECT
    ST_Union(geom)  geom,
    sum(val)        ct_values
FROM tbl_test
GROUP BY ST_Centroid(geom)

The expected output is to have something like the following image. There's no polygon over another.They just touch each other (the reason why there're still differences in the border line color).


texts:
I'm trying to match polyline layer features up to polygon edges in a different layer quickly.
I have a huge extent of polygons (river) and separate lines (bank features) that I need to get to follow/match the polygon outlines.

Spatial adjustment not working, reshape/modify features will be too time-consuming with snapping each vertex...

Is there a quicker easier way? I have ArcMap 9.2.


texts: How to overlay a map of Europe on a map of Australia?

texts: I'd like to overlay the cities of Europe on Australia for distance/size comparison purposes using QGIS and I'm wondering how best to achieve this while preserving distance between the European cities?


texts: Changing longitude latitude to latitude longitude for coordinate display in ArcMap?

texts: How do you change status bar display in arcmap 10.0 from longitude and latitude to latitude and longitude?


texts:
MapServerQuerySublayer (Version 10) does not seem to support setting a definition query, although according to documentation it should using ImapServerSublayer2.LayerDescription as IMapserverSublayer. Nor does it support selections, i.e. rendering objects matching a query or selection differently.

What are our options when making an Engine application which contains a MapControl consuming only arcgisserver layers? In the application we would like to change an individual layer's definition query to only show object that matches the user's interest. We would also like to render some objects as "selected" to reflect user selections.

It would also be nice to be able to get the selected objects to the client as features, which also does not seem to be supported out of the box in ArcGisServer.

Any ideas on how to accomplish this?

Edit:

I don't think that we can wait for ESRI to fix this bug, so I'm looking for an alternate solution, which may include fetching the selection as features and use a local layer. Engine seems to lack support for fetching features from a map service without having to connect to the context itself (not using the Layer's connection) or fetching them one by one using identify. Using the Feature service seems like a good idea, but Engine has no built in functionality to from a Layer fetch IFeatures from the FeatureService. I think that there would be a more graceful way than having to involve extra connections to the server, which implies more configuration in the application.


texts: Encoding issue in QGIS while loading CSV file

texts:
I am trying to create a value map on a field by loading a CSV file that was Exported from Google Sheets, but QGIS is loading unrecognized characters.

However, the Arabic characters do show up fine if I open up the CSV in notepad. As I understand, Google Sheets exports in UTF-8 which should be recognizable.

Thanks in advance.


texts: How to add a comment to a source field in QGIS 3?

texts:
How do I add a comment to a Field in the Layer Properties - Source Fields tab?
The documentation says "Comments can be added by clicking in the comment field of the column ...", but clicking in the Comment Field or the Alias does not work.

However, I may add an Alias using the Attributes Form.

Example: Adding "Dato" as an Alias.

But, how may I add a comment (no. "Kommentar")?
I see Comment below Alias in the Attributes Form window, but the "text field" is missing?


texts:
In ArcGIS Desktop I have a slope raster and calculated mean slope per watershed for a fish habitat model. There are zero ("0") values in some subcatchments where no information was available to calculate the mean.
This subcatchments have to get the value "NoData", otherwise my model will think zero is a valid input.

How can I change this zero values to NoData?

Before I used the "reclassify" tool. But this does not allow me to keep the decimal numbers of my many means. It only allows for categories.
I tried at the symbolgy tab, I can exclude the zeros there and then the maps looks like what I need. But I am afraid this is only the appearance and the zeros are still in the data.


texts:
New user to ArcGIS Pro here.
I was wondering if it was possible to export/convert the ArcGIS Pro .stylx file type to another format? Nothing specific, just another file type so that it can be used elsewhere.
I have seen a few links like the one below about converting it back to .style but its still a Esri format.
Converting styles from ArcGIS Pro (.stylx) to ArcMap (.style)?


texts: When you query the postgis database, you will get a number of row of record, how to write a sql query to get a valid geojson object instead of rows of record?


texts: How to get the BBOX of selected features in QGIS using the UI?

texts:
Is there an easy way to get the BBOX of the currently selected vector features in QGIS?

I've seen How to quickly calculate bounding box coordinates of a selected polygon with QGIS but it does not solve the problem if multiple features are selected.

I'm asking if QGIS has an out-of-the-box feature to do this.
I don't want to use the Python console, and I don't want to create a new layer just to get this information.

If it is not possible I'll just open a feature request.


texts: How to add a polygon coodinates field instead of drawing it

texts:
I am working on a simple web mapping app using Django framework, and I am trying to add a geometric field where to insert polygon geometry's coordinates instead of drawing it manually (these coordinates will be stored in a PostgreSQL database under the column geometry). How to make geom field a coordinates insertion field?
Here is my model class:
from __future__ import unicode_literals
from django.db import models
from django.contrib.gis.db import models as gis_models
from django.db.models import Manager as GeoManager
from django.contrib.gis.geos import Polygon
class projets(models.Model):
     titre = models.CharField(max_length=15)
     avancement = models.CharField(max_length=50,null=True,blank=True)
     geom = gis_models.MultiPolygonField(srid=4326,null=True,blank=True)
     class Meta:
         verbose_name_plural ="Projets"


texts: Connecting markers by following the roads in Google Maps

texts:
I have many markers in a kml file. All the markers are on the roads.

I wonder if there is a way to process the data and connect the markers by following the road? Of course, the markers should not be randomly connected with each other. The point is to mark the parts of the roads that have markers on them.

I could of course just manually create the lines myself. But there are so many markers, and the data is expanded almost daily.


texts: How to apply k-means algorithm on image in pyqgis?

texts:
I try to apply k-means clustering in qgis console like below shown in this site.Unsupervised KMeans image classification

processing.runalg('otb:unsupervisedkmeansimageclassification', 'D:\data\sub_66', 128,100, 5, 1000, 0.0001, 'D:\data\sub_66kmeans')

but it gives error like this.


texts: Creating Leaflet layer that can be toggled in layer control from CartoDB viz.json?

texts: How do I create a leaflet layer that can be toggled in the layer control from a cartodb viz.json? I have read that the cartodb.createlayer does not create a leaflet layer. It seems you can add the cartodb table, but of course that takes away the whole point of using cartodb for easy visualization. This seems like a basic function, I am a novice in terms of coding so maybe I am missing something.


texts: CSV to Vector Tiles: Encoding polygons in Tippecanoe (or some other software)

texts:
I want to generate .pdf vector tiles from polygon features. From what I've found out, tippecanoe seems to be the way to go. Like so (found here):
tippecanoe --no-feature-limit --no-tile-size-limit --no-tile-compression --output-to-directory directory layerName.json

I'm pondering the optimal workflow. My data is provided in a database from which I can easily build CSV. Tippecanoe can read CSV but apparently for point features only. However, I need polygons. I could generate the input as a geocsv as specified here with the geometries provided in some way, e. g. as WKT (polygon examples from here).
id   name      attr1     attr2    geom
1    garden1   green     round    POLYGON ((30 10, 40 40, 20 40, 10 20, 30 10))
2    garden2   emerald   edgy     POLYGON ((35 10, 45 45, 15 40, 10 20, 35 10), (20 30, 35 35, 30 20, 20 30))

There is a file geojson.cpp in the tippecanoe Github but it doesn't look like it's using this specification and I can't find any documentation for it. So I wonder if there's a way to build vector tiles in tippecanoe (or any other software) without taking the detour of generating GeoJSON first.


texts:
I am trying to create a Meshlayer in memory and tried this simple test in QGIS 3.6:

from qgis.core import QgsProject,QgsMeshLayer
proj = QgsProject.instance()
uri="1.0, 2.0 \n2.0, 2.0 \n3.0, 2.0 \n2.0, 3.0 \n1.0, 3.0 \n---0, 1, 3, 4 \n1, 2, 3 \n"
outLayer = QgsMeshLayer( uri, 'test', "memory_mesh")
proj.addMapLayer(outLayer)

Apparently dataprovider has no attribute name and QGIS crashes.
Is there something that I am doing wrong? What shoulb be the format of URI, a simple string?


texts:
I am trying to use the GDAL merge tool and receive the following error message:

  The process failed to start. Either the invoked program is missing, or
  you may have insufficient permissions to invoke the program.

On install I checked the "Path to the GDAL Executables" and it already had the /usr/bin added to the install. Here is what I have in this field:

  /Library/Frameworks/GDAL.framework/Versions/1.9/Programs/usr/bin

Thanks in advance for any help with this issue


texts: Layer disappears on OpenLayers after specific zoom in

texts:
I have added some layer(s) and a base layer using OpenLayers v 2.11, layer does not appear when it first opened the page, but after zoom out  it does appear. Why is behaving in this way, are there any parameteres for always visible?

var map = new OpenLayers.Map("map");
var wms = new OpenLayers.Layer.Google("Google Streets", {
                numZoomLevels: 12
            });
map.addLayer(wms);
var wms2 = new OpenLayers.Layer.WMS("Basic1", "http://localhost:8080/geoserver/wms",
           {
               layers: 'SomeLayers',
               transparent: true,
               numZoomLevels: 12
           }, {
               isBaseLayer: false
           });
map.addLayer(wms2);
var ls = new OpenLayers.Control.LayerSwitcher();
map.addControl(ls);
map.addControl(new OpenLayers.Control.ZoomBox());
ls.maximizeControl();

map.setCenter(lonlat, 6);


texts: find and Calculate centroid within polygon either in ArcMap or in MapInfo

texts: How can I calculate single/multiple centroid(point) within polygon either in ArcMap or in MapInfo. Both object(point, polygon) are in different layers.


texts: Making DEM from contours without elevation data?

texts: I got a shapefile with contours, but no elevation data is present. We can assume as known info, that the contours mark elevations at known steps (say, 20 m). Is it still possible to restore the elevation information and build the DEM model of relief, assigning arbitrarily the 0 elevation value to some starting contour? (either in ArcMap, or in QGis/Grass)


texts: How to configure tilestache.cfg "projection"?

texts:
I am using tilestache-1.51.6 and mapnik 3.0.9.
The tilestache.cfg is

{
  "cache":
  {
    "name": "Disk",
    "path": "/var/www/cachedata",
    "umask": "0000"
  },
  "layers":
  {
    "test":
    {
        "provider": {"name": "proxy", "provider": "OPENSTREETMAP"}
    },
    "kerala":
    {
        "provider": {"name": "mapnik", "mapfile": "district.xml"},
        "projection": "WGS84"
    }
  }
}

and the district.xml is:

<?xml version="1.0" encoding="UTF-8"?> <!--<!DOCTYPE Map SYSTEM "mapnik3.dtd">-->
<Map background-color="#fff" srs="+proj=utm +zone=43 +datum=WGS84 +units=m +no_defs" buffer-size="128" minimum-version="0.7.2">

<Style name="District_style_line">
<Rule>
     <MaxScaleDenominator>6600000000</MaxScaleDenominator>
     <PolygonSymbolizer fill="#ebb1af"/>
         <LineSymbolizer stroke="rgb(0,0,0)" stroke-width="0.98" stroke-linejoin="round"/>
</Rule>
</Style>
<Style name="District_style_text">
    <Rule>
    <MaxScaleDenominator>6600000</MaxScaleDenominator>
    <TextSymbolizer face-name="DejaVu Sans Bold" size="8" fill="black" wrap-width="10" avoid-edges="false" placement="interior">[name]</TextSymbolizer>
    <MinScaleDenominator>2000000</MinScaleDenominator>
    </Rule>
</Style>
<Layer name="District_Boundary" srs="+proj=utm +zone=43 +datum=WGS84 +units=m +no_defs">
    <StyleName>District_style_line</StyleName>
    <StyleName>District_style_text</StyleName>
     <Datasource>
        <Parameter name="dbname">testdb</Parameter>
        <Parameter name="extent">291627,904686,958569,1426831</Parameter>
        <Parameter name="geometry_field">the_geom</Parameter>
        <Parameter name="host">localhost</Parameter>
        <Parameter name="password">password</Parameter>
        <Parameter name="port">5432</Parameter>
        <Parameter name="srid">32643</Parameter>
        <Parameter name="table">(select * from district) as districts</Parameter>
        <Parameter name="type">postgis</Parameter>
        <Parameter name="user">test</Parameter>
    </Datasource>
</Layer>
</Map>

I am not getting the tile properly, I think the projection is wrong.
What projection should I give ?


texts: Represent huge amount of polygons with Leaflet

texts:
I'm actually aiming to represent a huge GeoJSON file using Leaflet (JS). The GeoJSON size would be around 300 Mb, so I'm looking to raster the polygons while being at a higher zoom, and vectorizing (on the bound) as we zoom into the region I'm interested in.

I did some research, I found this extension Here
but it looks like its not maintained and it doesn't raster then overtime or over zoom vectorize the polygons.


texts: Unattractive image artefact when exporting from QGIS's Print Composer

texts:
When I export my map from Print Composer in QGIS 1.7.4, an odd red band appears across the middle of the image. This happens when I export to both PNG and PDF. Its not visible on the screen before export, or present in the original raster image.
The band is to the left of the island halfway down the coast, below Pwani.

Any ideas what's causing this and how to fix it?


texts: Delorme GPS and Delorme X map alternatives

texts:
We have a group doing field work and are having issues using Delorme GPS and Delorme X Map and Delorme X Map.

Wondering if there is a way to connect the delorme GPS pn 40 to ArcReader / ArcExplorer in the field?

All shape files used in the field are polylines ArcGIS 10.


texts: Upload shapefile to ArcGIS Server, using Python API

texts:
I have shapefiles (.zip) that I want to upload to an ArcGIS Server programmatically using the ArcGIS Python API.

All the examples I've found assume the data is already on the server. Is there a tutorial on uploading data to the server using the Python library?


texts: Formula of mapnik bounding box (spherical mercator)

texts:
I want to implement a tile server and I'm following the code in the repo node-mapnik-sample-code.

When a request is sent with the parameters x, y, z (longitude, latitude, zoom), those parameters are used to create a bounding-box, with the following code:

/**
 * Convert tile xyz value to Mapnik envelope
 *
 * @param {Number} x latitude number.
 * @param {Number} y longitude number.
 * @param {Number} zoom zoom.
 * @param {Boolean} tms_style whether to compute a tms tile.
 * @return Object Mapnik envelope.
 */
SphericalMercator.prototype.xyz_to_envelope = function(x, y, zoom, TMS_SCHEME) {
    if (TMS_SCHEME) {
        y = (Math.pow(2, zoom) - 1) - y;
    }
    var ll = [x * this.size, (y + 1) * this.size];
    var ur = [(x + 1) * this.size, y * this.size];
    var bbox = this.px_to_ll(ll, zoom).concat(this.px_to_ll(ur, zoom));
    return mercator.forward(bbox);
};

(full code here: https://github.com/mapnik/node-mapnik-sample-code/blob/master/utils/sphericalmercator.js)

Why the longitude and latitude have to be multiplied by the tile size (256) here:

var ll = [x * this.size, (y + 1) * this.size];

I'm getting lost between pixels and geographic coordinates.


texts: CartoDB Torque not showing year in slider

texts:
I'm trying to set up a map using Torque with CartoDB's in browser wizard.  I've set the time column to year, but the this does not show up on the slider.  Anyone know how to make year display in the slider rather than a count/ID number?

This is my test map:
https://jamaps.cartodb.com/viz/e956c248-5cd0-11e5-ae08-0e6e1df11cbf/public_map


texts: Digitized polygon not lining up with others of same coordinate system

texts:
To start: I created a new polygon shapefile to represent an old site boundary. I set the spatial reference at the time the shapefile was created. I then started an editing session and digitized a polygon from a JPG image (I set the spatial reference for this image too).

I need to overlay the old site boundary (the polygon I just created) on top of an updated site boundary polygon (already existed). My problem is that when I add the polygon I just created to the map, it doesn't show up with everything else. It has the same coordinate system and everything as the rest of the map.

Any suggestions? I'm sure it's an easy answer but I've tried everything I can think of at the moment.

EDIT:
Software: ArcMap 10.7.1

Spatial reference of digitized polygon:

Projected Coordinate System:    NAD_1983_UTM_Zone_12N
Projection: Transverse_Mercator
False_Easting:  500000.00000000
False_Northing: 0.00000000
Central_Meridian:   -111.00000000
Scale_Factor:   0.99960000
Latitude_Of_Origin: 0.00000000
Linear Unit:    Meter

Geographic Coordinate System:   GCS_North_American_1983
Datum:  D_North_American_1983
Prime Meridian:     Greenwich
Angular Unit:   Degree

Spatial reference of preexisting polygon/map canvas:

Projected Coordinate System:    NAD_1983_UTM_Zone_12N
Projection: Transverse_Mercator
False_Easting:  500000.00000000
False_Northing: 0.00000000
Central_Meridian:   -111.00000000
Scale_Factor:   0.99960000
Latitude_Of_Origin: 0.00000000
Linear Unit:    Meter

Geographic Coordinate System:   GCS_North_American_1983
Datum:  D_North_American_1983
Prime Meridian:     Greenwich
Angular Unit:   Degree


texts: Using Baltimore City Survey Control System?

texts:
A client has provided some surveyed points referenced to the "BALTIMORE CITY
SURVEY CONTROL SYSTEM", they're all projected and given in easting/northing numbers. I need to reconcile them with some GPS coordinates I have in the NAD83 (2011).

I'm new to the GIS world, and really don't want to ask the client any stupid questions, but I can't figure out what the transformation is going to be.

Should I ask them for an EPSG code?


texts: Guidelines for organizing Python Toolboxes (.pyt) in ArcGIS

texts:
For one of my tasks I need to create a Toolbox with three tools (python scripts).

Instead of ordinary tbx I've decided to write it as Python Toolbox (pyt).

Everything is clear except the "good" way to organize the code.

As I have three tools it is not good way to store them in one file (pyt). So, I've decided to keep each tool in a separate .py file.

Here is a problem: what is the good way to organize several files with PYT for distribution or deployment on ArcGIS for Server? Should I keep them on the same level with PYT file or should I place them in some subdirectory (i.e. "Tools")?

Can you recommend any Esri guidelines or reference "big PYT toolbox" sample?

I have not find anything on this topic. In version 10.0 there was so called ToolShare folder structure which I've used.

This is more a question of style of coding. Because the idea to create one PYT file with 500-1000 or more lines of code does not look good to me and I believe it is not "pythonic" way.


texts: I'm sorry because my english isn't very good, I hope you're able to understand it. I'm using arcgis 9.3 right now, but when i try the Tools dropdown, then extension, there are no extensions listed. Do you have any idea how to solve this problem? I want to use the network analysis extension to finish my project, please help me.


texts: How to visualize shortest path that is calculated using Networkx?

texts: I have a graph G (made from a road network shapefile), and a source node S, and a destination node D. I have calculated the length of shortest path using single_source_dijkstra_path_length, but now I need to save the actual path in a shapefile. This will not help, because it returns the graph representation of my network, which looks nothing like my original network. Is there a way to save this path separately as a shapefile?


texts: How often are GSMaP and IMERG in Google Earth Engine updated?

texts:
I am trying to get a near-real-time data from Google Earth Engine for Hourly Precipitation. How often are the GSMaP and IMERG updated?
At what time (UTC) will new data become available? Below is an example of the latest available images from IMERG and GSMaP obtained on May 18th 06.00 UTC


texts:
I want to create  PYQGIS plugin and for my task i need my QGIS GUI form to have one more button except to original button ok.

the additional button to do something and original button to do something else for my script.

first i add new button in initGui :

def initGui(self):
    icon_path = ':/plugins/twoobutton/icon.png'
    self.add_action(
        icon_path,
        text=self.tr(u'twoobutton'),
        callback=self.run,
        parent=self.iface.mainWindow())
    self.dlg.pushButton.clicked.connect(self.run1)

and after I create new def run1:

def run1(self):
    self.dlg.show()
    result = self.dlg.exec_()
    if result:
        myvar = str(self.dlg.lineEdit_5.text())
        if myvar==1000:
            text1='hello'
            text2='world'
            self.dockWidget.dockWidgetContents.lineEdit.setText(str(text1))
            self.dockWidget.dockWidgetContents.lineEdit_2.setText(str(text2))
        pass

update :

class fff:
    """QGIS Plugin Implementation."""

    def __init__(self, iface):
        self.dlg = fffDialog()
        self.dlg.pushButton.clicked.connect(self.run1)

        """Constructor.

        :param iface: An interface instance that will be passed to this class
            which provides the hook by which you can manipulate the QGIS
            application at run time.
        :type iface: QgsInterface
        """
        # Save reference to the QGIS interface
        self.iface = iface
        # initialize plugin directory
        self.plugin_dir = os.path.dirname(__file__)
        # initialize locale
        locale = QSettings().value('locale/userLocale')[0:2]
        locale_path = os.path.join(
            self.plugin_dir,
            'i18n',
            'fff_{}.qm'.format(locale))

        if os.path.exists(locale_path):
            self.translator = QTranslator()
            self.translator.load(locale_path)

            if qVersion() > '4.3.3':
                QCoreApplication.installTranslator(self.translator)

        # Declare instance attributes
        self.actions = []
        self.menu = self.tr(u'&ffffff')
        # TODO: We are going to let the user set this up in a future iteration
        self.toolbar = self.iface.addToolBar(u'fff')
        self.toolbar.setObjectName(u'fff')

But if I click on pushbutton then don't show me anything.
How to fix that?


texts: How to cartographically represent active faults, without using geological symbols

texts: I would like to ask how to cartographically represent active faults, without using geological symbols. I don't want to make a geological map, that's why I don't use that kind of symbology. I don't have the ages of faults, but I do have the type. I saw how the United States Geological Survey (USGS) made similar maps, but they used "age" instead of "type" (of faults).


texts: Can't display ONLY my GeoServer layer using OpenLayers 6

texts:
I am having some trouble getting my OpenLayers 6 maps to display my own layer from GeoServer. The code below let me choose between OSM , Bing Map or CartoDb as a base layer and displays a default GeoServer layer, topp:states!
And all this works fine.

But in the end, as I attempted to display my own GeoServer layer "grupo" created from PostGIS DB, nothing appears. The GetMap request is working fine as shown here:

so it must be something wrong within the layer or projection (my guess)...

This layer "grupo" was created with EPSG:3763. There are no errors in the developer console logs, if I open the GeoServer log there is record of my own GeoServer layer "grupo" and topp:states being requested!

My GeoServer layer "Setup"

Layer Data

Tile Caching

LayerPreview

LayerPreview url -> http://localhost:8080/geoserver/sig/wms?service=WMS&version=1.1.0&request=GetMap&layers=sig%3Agrupo&bbox=-933103.256159661%2C5052848.3016269%2C-931781.706532744%2C5053702.88175803&width=768&height=496&srs=EPSG%3A3763&format=application/openlayers

Here's the code:

function init(){

    const map = new ol.Map({
            view : new ol.View({

               center: ol.proj.transform([-8.3770108, 41.2771817],'EPSG:4326','EPSG:3857'),
                    //center: [-932510.7343462185,5053142.307944268],
                zoom: 15,
                maxZoom: 20,
                minZoom: 4

            }),
            target : 'js-map'
    })

    // Base Maps
const openStreetMap = new ol.layer.Tile({
    source : new ol.source.OSM(),
    visible: true,
    title: 'OpenStreetMap'
})

const bingMap = new ol.layer.Tile({
    source: new ol.source.BingMaps({
            key: 'AmsC0864kukoBixOpup8lsv3yFO_qzXenrScHMmdHz5HVDSI0wZsGcP9DyW9zIU-',
            imagerySet: 'Aerial'
    }),
    visible: false,
    title:'BingMap'
})

const CartoDB = new ol.layer.Tile({
    source: new ol.source.OSM({
        "url" : "http://{a-c}.basemaps.cartocdn.com/dark_all/{z}/{x}/{y}.png"
    }),
    visible: false,
    title:'CartoDB'
})

//layer group
const baseLayerGroup = new ol.layer.Group({
    layers: [
        openStreetMap, bingMap, CartoDB
    ]
})

//Layer Switcher (BaseLayers)
const baseLayerElements = document.querySelectorAll('.baseSideBar > input[type=radio]');
for(let baseLayerElement of baseLayerElements){
    baseLayerElement.addEventListener('change', function(){
        let baseLayerElementValue = this.value;
        baseLayerGroup.getLayers().forEach(function(element, index, array){
        let baseLayerTitle = element.get('title');
        element.setVisible(baseLayerTitle === baseLayerElementValue);
        })
    })
}
map.addLayer(baseLayerGroup);

 // layers geoserver
 var wms_layer = new ol.layer.Tile({
    source: new ol.source.TileWMS({
    url: 'http://localhost:8080/geoserver/topp/wms',
    params: {'LAYERS': 'topp:states','VERSION':'1.1.1', 'TILED': true},
    serverType: 'geoserver'
    })
   });
     map.addLayer(wms_layer);

 var wms_layer1 = new ol.layer.Tile({
    source: new ol.source.TileWMS({
    url: 'http://localhost:8080/geoserver/sig/wms',
    params: {'LAYERS': 'sig:grupo','VERSION':'1.1.1', 'TILED': true},
    serverType: 'geoserver'
    })
   });
      map.addLayer(wms_layer1);

    }


texts: Auto-Styling QGIS vector layers symbology

texts:
Is there a way to get QGIS3 to automatically apply specific symbology based on the layers name for example?

My situation:
I digitize maps often based off satellite imagery. I draw the map using one set of symbology, when I wish to produce the final map I must manually go through each layer and change the symbology for the output map. This is extremely tedious and quite often there is small variations between maps, as each time I set the symbology it is not exactly the same as last time.

I have started saving styles and loading them up for each layer which is an improvement, but I still have to go through each individual layer and load the style. I often have hundreds of layers.

Each layer has different symbology.


texts: How to fix garbled fields when importing csv to Arcmap

texts:
I'm having trouble reading fields in Chinese when importing a csv data file to ArcMap. I imported the file by using 'add data' option in the tools section. I'm confused because some csvs work while some do not, meaning for some csvs the characters are clearly read in ArcMap while for the other csv all of them are garbled.

These files were all edited in R and Python and they were encoded as utf-8. They both do not get garbled characters when read in R, so I'm not sure why ArcMap would not display them correctly.

The garbled attribute table looks like this:

I tried turning the csv file to xlsx format and then tried importing it in ArcMap but I could not import the file. I'm using ArcMap 10.5.

Could it also have to do with long texts like one in the image? I have long descriptions of each data points (around 500+).


texts: How to get elevation value in CSV format in R from raster SRTM DEM file?

texts:
I am trying to get elevation value in meter using R programming.

I have loaded raster library in rstudio:

elevation <- getData("alt", country = "BD")

How I can get those value with lat lon?

Here is the documentation for the function I used: getData


texts: QGIS virtual layer path never saved as relative

texts:
When I try to use virtual layers in QGIS, I'm unable to use relative paths.

I've already set in the "Project Properties", in the tab "General", the "Save paths" as relative.

You can see a screenshot to illustrate the issue.

Removing the absolute path "breaks" the virtual layer.
For this reason, my *.qgs project is not portable...

How do you manage this use case?
It's due to my way of doing, is it a feature, a missing feature or a bug?

Any advice to bypass this issue is welcome


texts:
I have some multi-spectral data, where the individual bands are separate Raster files.

How do I layer stack these rasters so that I get one raster with 3 bands instead of 3 different rasters?

This is quite easy to do in ERDAS and ArcGIS, but I haven't figured out how to do it in QGIS


texts: mapbox gl expression to match multiple values

texts:
i have a layer that i want to filter by multiple values stored in array

let filter = ['in','id', 1,2,3]
map.setFilter('layer_name',filter)

this works fine and it will return only features with ID 1,2,3
but in Mapbox GL specs, they label this in expression as deprecated
https://www.mapbox.com/mapbox-gl-js/style-spec/#other-filter
and encourage using the new expressions, i think the expression match would work but don't know know it should look like, i tried to do like

['match', 'get', 'id', 1,2,3 ];

it's not really that clear in the Docs so if anyone can clarify it please


texts: Using Conditional Statement in ArcGIS Deasktop?

texts:
I'm having trouble performing calculations with single output map algebra ArcgIS 9.3 using conditional statements and applying it to several model outputs as they are generated: The statement is:

Con(dem%n%<= -0.02,1)

I have tried these:

Con('dem%n%' <= -0.02,1)

and

Con("dem%n%" <= -0.02,1)

but they seems unstable. What would be a better solution?


texts: Is there a good guide out there for setting MapServer up on an existing windows apache server?

texts: I have been interested in trying out MapServer and have no problem setting up a windows development box.  Seems like all the MapServer implementations on Windows include Apache, php, mapserver, and a variety of other software.  Unfortunately, if I wanted to use Mapserver in my production environment, I do not want to uninstall my functioning and stable Apache to replace with the packaged mapserver download.  I understand that I can take pieces out of the package and use with my existing Apache but I would really like to have a clear set of instructions rather than references from 4 different websites.  I am looking for clear instructions/how-to in 1 place.


texts: Query to return closest LineString given a Point

texts:
I have a road lines data in a table with field
roads(id SERIAL ,the_line geometry(LineString,4326))

suppose an accident occurs in a middle or at any location in that waterline i want to stop the entire waterline
so, given a location or coordinates(POINT(350205.805 148722.56)),  I need to find a roadline which is closest to that point. and return its start point and endpoint

 "SELECT id, ST_AsText(ST_StartPoint(the_line )) AS startpoint,ST_AsText(ST_EndPoint(the_line )) AS endpoint from roadlines where (some condition to get the closest roadline) close to 'POINT(350205.805 148722.56)'"

POINT(350205.805 148722.56) WILL BE THE INPUT


texts: ee.Algorithms.If throwing both cases in Earth Engine

texts:
I'm trying to work with the ee.Algorithms.If in the sense that:
if list = 0(i.e. not containing any points within a buffer) -> do a function()
else (buffer is containing points)-> do a function()

My problem is that both my true and false cases are thrown which leads to nothing is happening.
I know about the server side vs client side issue and that that is one of my problems. However, I don't know how to deal with it. Maybe I can get a helping hand on that? see code below.

var VMAsize = VMA.size();
ee.Algorithms.If(VMAsize= 0, size0 , sizeint);
  var sizeint = function(){
    var VMA_regionobj_list = VMApoints.toList(VMA.size());
    var VMA_regionname_list = tools.featureCollection.propertyValues(VMApoints, 'Title');
    var VMA_regionsummary_list = tools.featureCollection.propertyValues(VMApoints, 'Summary');
    var VMA_regionupdatetime_list = tools.featureCollection.propertyValues(VMApoints, 'UpdateTime');
    var VMA_regionauthor_link_list = tools.featureCollection.propertyValues(VMApoints, 'Author');
    VMA_regionpanel.clear();
    var VMA_regionselect = ui.Select({
      items: VMA_regionname_list.sort().getInfo(),
      onChange: function(key) {
        var index = VMA_regionname_list.indexOf(key);
        VMA_regionpanel_info.style().set('shown', true);
        VMA_regionname.setValue(key);
        VMA_regionsummary.setValue('Announcement: '+VMA_regionsummary_list.get(index).getInfo());
        VMA_regionauthor_link.setValue('Web page: ');
        VMA_regionupdatetime.setValue('Update time: '+VMA_regionupdatetime_list.get(index).getInfo());
        VMA_regionlink.setValue(VMA_regionauthor_link_list.get(index).getInfo());
        VMA_regionselect.setValue(null, false);
      }
    });
  VMA_regionselect.setPlaceholder('Choose an VMA...');
  VMA_regionselect.style().set({margin: '0 0 1% 0', width: '98%'});
  VMA_regionpanel.add(VMA_regionselect);
  return VMA_regionselect;
  };

  var size0 = function(){
    VMA_regionpanel.clear();
    var VMA_regionselect = ui.Select({
      onChange: function(key) {
        VMA_regionpanel_info.style().set('shown', true);
        VMA_regionsummary.setValue('Announcement: There are no announcements within this area');
        VMA_regionauthor_link.setValue('Web page: ');
        VMA_regionselect.setValue(null, false);
      }
    });
  VMA_regionselect.setPlaceholder('Choose an VMA...');
  VMA_regionselect.style().set({margin: '0 0 1% 0', width: '98%'});
  VMA_regionpanel.add(VMA_regionselect);
  return VMA_regionselect;
  }


texts:
I put a pourpoint at the end of a headwater stream, then another further downstream that is fed by that headwater stream plus two others.

Obviously the pourpoint further downstream is going to share a watershed with the headwater pourpoint.

Is there a way of identifying the individual watersheds so that the watershed of the pourpoint further downstream includes that headwater?

I want to be able to do many pourpoints together and get the real watersheds, but they always are cut off by the overlaps.

I DO NOT have Arc Hydro because I cannot get it to install properly and haven't been able to find instructions on how to do so - see https://gis.stackexchange.com/questions/137538/installing-arc-hydro-problems.

I'm using ArcMap 10.2.2


texts:
Given a set of 2D or 3D points:

How to find the center of geometry of an object?

According to the following figure, the center of geometry differs from the center of mass if it is calculated in the simplest form i.e., homogenous density of mass. The problem appears, indeed, in the computation of those.
Commonly, one approach is to average X coordinates and Y coordinates separately i.e., find an average position to the given points (here in 2D). This may be used as centroid for the set of points representing an object.
As shown, due to the extra vertex along the bottom edge, for a simple rectangle the resulting centroid is (0.5,0.4) while the correct answer is (0.5,0.5).
Note that the example given is too simple. The problem of interest however is for complex shapes in 2D and objects in 3D for which only coordinates of vertices are available.
BTW, an efficient computational way is of interest.

Just to mention that I have checked some web links such as Wikipedia's however my current problem is that there are group of 2D and 3D points wish to find a point as representative for those. Thus centroid became of interest. The points are given without any topological information. You may consider them as point cloud. The demonstration here provided to make it clear that the commonly known averaging of coordinates (see for example this Stack Overflow Q&A) may be incorrect as shown in the example.

Here are some implementations for comparison:

aa    = accepted answer below
chull = convex-hull of points i.e., the golden polygon
cent  = centroid proposed in Wikipedia and discussed in aa as the polygon centroid
centl = centroid of polyline as explained in aa

Visually, centl looks better representative for the geometry given compared to cent. Two others look promising here but usually they are too biased if the dispersion of the points were inhomogeneous as it is a usual case.
And also consider that although convex-hull makes problem reasonably simpler however it may generate too long and too short edges without any symmetrical positioning in the space, that is, awareness is necessary if you do simple averaging (i.e., without weighting) for both cases: entire points (green) or convex-hull polygon vertices (blue).

One application could be found in Finding minimum-area-rectangle for given points?.


texts:
I want to automatically count features in my opened layers, and add the result in the layer name, with pyqgis. I can do it with this code :

layers = QgsProject.instance().mapLayers().values()
for layer in layers:
   fc = layer.featureCount()
   fc = str(fc)
   layer.setName(layer.name()+" ["+fc+"]")

However, it does count features in each layer, but not in each renderer category when a layer have a categorized style. I was trying something like this just to print the result :

for layer in layers:
    renderer = layer.renderer()
    renderert = renderer.type()
    if renderert == "categorizedSymbol":
        for cat in renderer.categories():
            print(cat.label+cat.featureCount())

But i get an "AttributeError: 'QgsRendererCategory' object has no attribute 'featureCount'"

Do you know how to apply the featureCount() function to layers and categories?


texts:
I am trying to iterate over the values of a raster dataset (band1). I am able to do it in python with the following code snippet (sorry I can't provide the original raster)

import numpy as np
import gdal

path = "data/Isle_wight.tif"

ds = gdal.Open(path)
myarray = np.array(ds.GetRasterBand(1).ReadAsArray())
print(myarray.shape[0])
#print(columns)
for j in range(myarray.shape[0]-1):
    for i in range(myarray.shape[1]-1):
        print( myarray[j][i])

My goal is to emulate the same in C-GDAL,

here is my code snippet (which does not work)

#include "gdal/gdal.h"
#include "gdal/cpl_conv.h"
int main()
{

    GDALDatasetH raster;
    GDALAllRegister();
    raster = GDALOpen( "/home/roger/Documents/98_GDAL_C++_snippets/data/Isle_wight.tif", GA_ReadOnly);
    if(raster==NULL)
    {
    printf("Invalid Raster Dataset");
    }
    else
    {
    GDALRasterBandH raster_band;
    int XSize, YSize;

    raster_band = GDALGetRasterBand( raster, 1);
    GDALGetBlockSize( raster_band, &XSize, &YSize);
    float *pafScanline;
    int nXSize = GDALGetRasterBandXSize (raster_band);
    int nYSize = GDALGetRasterBandYSize (raster_band);
    pafScanline = (float*) CPLMalloc(sizeof(float)*nXSize);
    int address_of_pafScanline = *pafScanline;
    int band1[nXSize][nYSize];
    int band1 = GDALRasterIO(raster_band, 0, 0,GF_Read, nXSize, 1, pafScanline, nXSize, 1, GDT_Float32, 0, 0);
    for(int i = 0; i <= nXSize; i++) {
        for(int j = 0; j < nYSize; j++ ) {
            printf("%d", band1[i][j]);
        }
    }
}

}

Question: What does GDALRasterIO returns?, I did an array declaration when writing in line 25:

int band1[nXSize][nYSize],

And then I read the raster Data and assign it to the previous band1 array in the following line (26):

int band1 = GDALRasterIO(raster_band, 0, 0,GF_Read, nXSize, 1, pafScanline, nXSize, 1, GDT_Float32, 0, 0);

When compiling through

cc open_array.c -o open_array -lgdal

I got the following error:

open_array.c: In function ‘main’:
open_array.c:26:9: error: conflicting types for ‘band1’
     int band1 = GDALRasterIO(raster_band, 0, 0,GF_Read, nXSize, 1, pafScanline, nXSize, 1, GDT_Float32, 0, 0);
         ^~~~~
open_array.c:25:9: note: previous declaration of ‘band1’ was here
     int band1[nXSize][nYSize];
         ^~~~~
open_array.c:29:31: error: subscripted value is neither array nor pointer nor vector
             printf("%d", band1[i][j]);

Does that happen because of the declaration in line 25?


texts: What is the best public geoservice for getting data about buildings?

texts:
What is the best public geoservice for getting data about buildings, which is free? Or maybe it's better to ask "does exist such a public service?"

I want to add some 3d polygons to my 3d plane, which I render in WebGL project and I want to get such a service, which may help me with getting detailed information about buildings.

What protocol will respond me an answer on query (WFS or another) - doesn't matter.
I can't google such a service, so ask you in this topic about it.


texts: Get a specific value from a table with using a few conditions in Python

texts:
I have a table that each different values in the table is the result of a combination of different conditions. For example in the below image the condition is as follow: if coverType=fallow , Treatment=Crop Residue Cover, Condition/ImperviousArea=Good, SoilType=C then the value is equal to 83.

I want to have a tool which asks the user to choose a value from each column (e.g. choose CoverType; SoilType, ...) and then return the related number as output. Do you have any thoughts how should I do this?

So far I have just the first lines of the code as below:

import arcpy
path= r'H\Python\FinalProject\RunOff.gdb'
arcpy.env.workspace= path
arcpy.env.overwriteOutput = True

table=r'H\Python\FinalProject\RunOff.gdb\CN.csv'
cursor= arcpy.SearchCursor(table)
#CN= ?? # I don't know this part!
#S=(1000/CN) - 10
#I = 0.2(S)
#P = arcpy.GetParameterAsText(0)
#out_feature_class = arcpy.GetParameterAsText(1)


texts: Identifying the best OLS regression model

texts:
I have 10 independent variables to include in an OLS regression model.
I know how to interpret the results and that for comparing models, the AICc value can be used, where a smaller value indicates a better model (ArcGIS Help).

But to find the best OLS model I cannot try all different combinations of the independent variables myself. So if I start with all 10 independent variables for the first model, how should I define, which independent variables to keep/remove to get a better model? Removing those variables where the coefficient is small? Removing those variables that are not statistically significant?


texts:
My road network is broken up into segments. I want to produce a map that shows where the segments start and finish.

I don't really want to create a point layer that shows the location of the ends of each line, but will if I have to. And I know that I could categorise segments and use a colour ramp or random colours, but I was looking for something that would look a little more classy.

I was wondering if there is a way to use a colour ramp - or gradient fill - for lines in QGIS? What I envisage is every line segment starts with one colour and finishes with another colour. This would also serve to indicate line direction - another one of my requirements.

Any ideas?
Johanna


texts: Produce raster map using universal kriging as interpolation method

texts:
I want to do kriging with external drift(universal kriging) in R. The covariate I use is in WGS84 refrence system. Because of nature of kriging which calculates distance in planar surface I have to do transformation from WGS84 to ETRS89. what I want out of  kriging is a raster map. Because of the reference system transformation the prediction points are not regular grids and what I get as a result is a vector map.

I used 'gstat' package in and spplot() function to create plots. But When grids are regular (wich I can not define because of nature of the covariate data) with equal distance I have a raster map.


texts: Get latitude and longitude from just the name of a company, its city and country of location, no physical address

texts:
I have been trying to find the latitude and longitude of a list of company-names. For example, the addresses look like this-(1Mobility India Pvt Ltd Pune Maharashtra India), (21st Century Software Solutions Pvt Ltd Visakhapatnam Andhra Pradesh).

I am geocoding using python and I know how to geocode if I am given the complete set of addresses.

For that I use geocoders package from GoogleV3. Using this package however I am able to recover only the latitude and longitude of the city and not the exact location of the company within the city. In Google maps online , however, I can recover the exact lat long for the company.

Does anyone have any experience about which python package to use for batch geocoding when only company names are given or any other way to batch geocode with just names?


texts: Total newbie. Creating map from aerial photo of golf course in QGIS

texts: I have a aerial photo in tif with a tfw file of the golf course that I work at.  I would like to create vector data of different features on the course such as greens fairways tees, etc. to calculate areas and distances to start with.  So far I have been able to load the tif into QGIS but after that I am completely lost.  It seems that every tutorial is using data sets that are complete.  I imagine I am creating new data using the raster data? I'm lost.  Thanks for any help


texts: I want to "co-interpolate" temperature with DEM. I have SRTM 90m aggregated to 1km and I want to "co-interpolate" my DEM with points with temperature data. How do I do this? I am using ArcGIS 10.


texts: Create New field equal to existing field - convert to integer - PyQGIS 3.4

texts:
I'm using PyQGIS 3.4, I am trying to add a new field and set the values equal to an existing field. The existing field are population totals but are in string format, I want to convert these to integers or real numbers.

I can add the new field however the problem is adding values to the field.

I am getting the following error.

exec(code, self.locals)
  File "<input>", line 1, in <module>
  File "<string>", line 13, in <module>
TypeError: QgsExpression.prepare(): argument 1 has unexpected type 'int'

Script:

from qgis.core import *

from PyQt5.QtCore import *

output_layer = iface.activeLayer()

provider = output_layer.dataProvider()

provider.addAttributes([QgsField('new_pop', QVariant.Int)])

output_layer.updateFields()

output_layer.startEditing()

fieldIndex = output_layer.fields().indexFromName( 'new_pop' )

e = QgsExpression( 'toInt("exist_pop")' )

e.prepare(output_layer.fields().indexFromName( 'new_pop' ))

output_layer.commitChanges()


texts: Bounding box of set of MULTILINESTRINGs in PostGIS?

texts:
I have a PostGIS table with a geometry column (geom) and a type column. The geom column contains only MULTILINESTRINGs and the type column contains tags that tell me what each MULTILINESTRING represents. The tags are just letters. I’m trying to get the bounding box of all the MULTILINESTRINGs in the table that have type='a'.
I tried the following but it does not work:

SELECT *
FROM line_table
WHERE line_table.geom && SELECT ST_Envelope(line_table.geom) FROM line_table WHERE type='a';

I’m new to PostGIS and SQL in general so I’m not sure if I’m even close to the answer. Right now I think the problem is that I'm getting one independent envelope for each MULTILINESTRING where type='a' instead of the bounding box of all the MULTILINESTRINGs together. Perhaps there is a way of merging all the MULTILINESTRINGs of type='a' into one geometry but so far I have not been able to figure out how.


texts: Opening vector layer in QGIS Python standalone application?

texts:
I'm working on a project using the QGIS Python API.  I seem to be having some issues with opening a vector layer from a shapefile in my standalone application.  This is the Python script I'm using to test opening the layer:

#   Import statements
from qgis.core import *
import os

#   Instantiate the QGIS Application
GUIEnabled = True
app = QgsApplication([], GUIEnabled)
#   Update prefix path
app.setPrefixPath("C:\OSGeo4W64\apps\qgis", True)
app.initQgis()

#   File to be read by the application
vector_source = "absolute_path_to_file"

#   Read layer from file
layer = QgsVectorLayer(vector_source, "Vector_Layer", "ogr")

#   DEBUG: Test if the file path points to a valid file
print "File Path points to a file: ", os.path.isfile(vector_source)

#   DEBUG: Test if the layer initialized correctly
print "Output of layer.isValid(): ", layer.isValid()

layer.isValid() returns false when I run the script.  I checked to make sure the file was valid by opening it in the QGIS GUI, which worked.  I also checked out the answer to Creating QGIS layers in python console vs stand-alone application, but the answer there is specific to Mac users, and I am running on Windows 8.


texts: Transform an image with projection EPSG:4326 to EPSG:3857

texts:
I have an image with known boundary and projection EPSG:4326 and want to overlay it to openstreetmap using Leaflet.js.

I have tried using function imageoverlay, but the image is not aligned with the basemap.

L.imageOverlay(imgUrl, [[-15, 75], [45, 145]], {opacity: 0.6, autoZIndex: true});

So, I use gdalwarp to transform the image projection to EPSG:3857.

gdal_translate -of Gtiff -a_ullr -15 145 45 75 -a_srs EPSG:4326 test.png out_4326.tiff
gdalwarp -s_srs EPSG:4326 -t_srs EPSG:3857 -ts 850 728 out_4326.tiff out_3857.tiff
gdal_translate -of png out_3857.tiff final.png

However, the final image looks like distorted. Here is the original image.

Here is the final image


texts: Cannot pan to correct location using IEnvelope.CenterAt(IPoint)

texts:
It seems no matter what projection I try it pans me to the Atlantic Ocean. I have run across this before in an ArcPy script and it was due to being in a projected coordinate system and thus it was reading the coordinates as x meters away from the origin, but to my understanding a geographic coordinate system would read them differently. Whenever I use GoToXY with these coordinates it puts me in Texas. (10.1 SP1, C#.NET, VS2010, Desktop Addin)

IMxDocument mxdoc = ArcMap.Application.Document as IMxDocument;
IEnvelope pEnv = mxdoc.ActiveView.Extent;
ISpatialReferenceFactory srFactory = new SpatialReferenceEnvironmentClass();
int gcsCode = (int)esriSRGeoCSType.esriSRGeoCS_WGS1984;
ISpatialReference spatRef = srFactory.CreateGeographicCoordinateSystem(gcsCode);

try
{
    IPoint point = new PointClass();
    double xCoord = 33.217918;
    double yCoord = -98.158749;
    point.SpatialReference = spatRef;
    point.PutCoords(xCoord, yCoord);

    pEnv.CenterAt(point);
    mxdoc.ActiveView.Extent = pEnv;
    mxdoc.ActiveView.Refresh();
}
catch (Exception ex)
{
    MessageBox.Show(ex.Message);
}


texts: How to add GeoServer WMS layer to Mapbox GL API?

texts:
I have been using GeoServer WMS layers with OpenLayers 3 and it is working well.
Below is the code how I have used to add my GeoServer WMS layer to the Mapbox GL API:

   var mapOptions: MapboxOptions = {
  container: 'map',
  style: 'mapbox://styles/mapbox/light-v10',
  zoom: 17.5,
  center: [0.5500863948437313, 53.229982583617584],
  pitch: 60
}

this.mapInstance = new Map(options);
var layerConfig: mapboxgl.Layer = {
  'id': 'wms-test-layer',
  'type': 'raster',
  'source': {
    'type': 'raster',
    'tiles': [
      'https://geoserver.visur.io/geoserver/****/wms?BBOX={bbox-epsg-3857}&SERVICE=WMS&REQUEST=GetMap&CRS=EPSG:3857&WIDTH=256&HEIGHT=256&LAYERS=tourmalineqa:GasWell&FORMAT=image/png'
    ],
    'tileSize': 256
  },
  'paint': {}
}

this.mapInstance.addLayer(layerConfig)

same with OpenLayers is working fine but in Mapbox GL API I am getting below error

  ERROR Error: Style is not done a the loading
      at i._checkLoaded (mapbox-gl.js:33)
      at i.addLayer (mapbox-gl.js:33)
      at r.addLayer (mapbox-gl.js:33)
      at MapComponent.push../src/app/map/map.component.ts.MapComponent.addWmsLayer (map.component.ts:363)

LOOKS like the official sample itself is not loading as their services are no longer available Mapbox-gl WMS sample


texts: Quick way to get WKT from a GeoDataFrame's geometries in GeoPandas?

texts:
Is there a quick way to produce a column/series containing the WKTs of all the geometries in a GeoPandas GeoDataFrame?
I got around this using  the apply function (see example below), but I just wanted to check if there is a better\quicker way to do this.
import geopandas as gpd
from shapely.geometry import Point

mydf = gpd.GeoDataFrame({'x':[1,2,3],
                         'y':['a','b','c'],
                         'geometry':[Point([0,0]),Point([1,1]),Point([2,2])]},
                        geometry='geometry')

def get_wkt(geom):
    return geom.wkt

wkt_col = mydf.geometry.apply(get_wkt)

print(type(wkt_col))
print(wkt_col)

The output I'm looking for (which the code above produces) is copied below:
<class 'pandas.core.series.Series'>
0    POINT (0 0)
1    POINT (1 1)
2    POINT (2 2)
Name: geometry, dtype: object

Edit
Here are some performance measures from the main suggestions:
# My original method:
mydf.geometry.apply(get_wkt)
267 µs ± 9.31 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

# martinfleis' method:
gpd.array.to_wkt(mydf.geometry.values)
57.9 µs ± 629 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

# gene's first method:
mydf.apply(lambda row:row['geometry'].wkt, axis=1)
461 µs ± 4.46 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

# gene's second method
[geom.wkt for geom in mydf.geometry]
61.6 µs ± 1.41 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)


texts: Lat/Lon to UTM conversion across multiple zones

texts: I have a data set that represents the extents of 1:25000 map sheets over a whole state. Its stored in PostGIS in WGS84 lat/lon. Id like to retrieve the MinX, MaxX, MinY and MaxY for the bounding box for each record. Easy... but... Id like to display it in UTM. Problem is that this state crosses multiple UTM zones. The DB should be smart enough to calculate the central meridian and give me a full UTM value for each extent. I can do this via code or in excel (shudder) but Id rather do this in the database so I can create a view. Any ideas/examples I can follow?


texts: Rasterio error: Input shapes do not overlap raster but reprojection not working

texts:
I'm trying to use rasterio (v1.0.13) and fiona to perform a raster clip on a geotiff using a geojson polygon. For reference, the clip works perfectly from the command line using GDAL:

gdalwarp -cutline tmp_yard.geojson  -crop_to_cutline -of GTiff -r cubic leaf_index.tif tmp.tif

But the following does not (as per rasterio docs: https://rasterio.readthedocs.io/en/latest/topics/masking-by-shapefile.html#)

with fiona.open("tmp_yard.geojson", "r") as geojson:
    features = [feature["geometry"] for feature in geojson]
    pprint.pprint(features)
    with rasterio.open("leaf_index.tif") as src:
        out_image, out_transform = rasterio.mask.mask(src, features,
                                                           crop=True)
        out_meta = src.meta.copy()

Produces error:

  [{'coordinates': [[(-11680153.379982425, 4825380.211452511),
                     (-11680101.881582342, 4825359.649286965),
                     (-11680099.47141429, 4825337.825410188),
                     (-11680153.677325355, 4825348.074481825),
                     (-11680153.379982425, 4825380.211452511)]],   'type': 'Polygon'}]

  ValueError: Input shapes do not overlap raster.

I tried to apply the solution here: Masking GeoTIFF file after GeoJSON through rasterio - "Input shapes do not overlap raster" to reproject my geoson, using:

import shapely
import pyproj
def project_wsg_shape_to_csr(shape, csr):
    project = lambda x, y: pyproj.transform(
        pyproj.Proj(init='epsg:4326'),
        pyproj.Proj(init=csr),
        x,
        y
      )
    return shapely.ops.transform(project, shape)

yard=project_wsg_shape_to_csr(features, 'epsg:4326')

However, this produces the error:

  AttributeError: 'list' object has no attribute 'is_empty'

Anyone got any ideas?


texts: Don't allow null values in Collector for ArcGIS

texts:
I have created a geodatabase in ArcMap with four Feature Classes for a Survey using Collector for ArcGIS.

I have multiple Fields per Feature Class and various Domains for single and multiple use in the Feature Classes. Some fields in the Feature Classes should not allow null Values. This is correctly implemented in the different Feature classes. As well after loading the Layers to ArcGIS online the don't allow null values seems to work.

If I now publish the Service and use it with Collector for ArcGIS it does not work any more. Which means, I can collect new data without an error if the fields contain nulls.
Any idea?


texts: Create shapefile of raster outline using python

texts:
I want to create a single bounding shapefile polygon of valid pixels in raster data (GeoTIFF) using python.

I need something efficient as I need to do this for hundreds of rasters. I had thought of writing out a mask using simple numpy manipulation, then passing that to gdal_polygonize.py, but that will not be quick enough.  Does anyone know of a more direct way?


texts: Rasterizing vectors in SAGA via saga_cmd - how to set the exact cellsize, extent (basicly grid system) based on other grid?

texts:
I'm doing some batch analysis using SAGA's saga_cmd. Basicaly I have a lot of grids and shapes (mainly generated in SAGA by vectorization of raster classifiaction results).

Now I'm trying to make thing's the other way:

I want to run saga_cmd grid_gridding 0 to rasterize polygons but I need the resulting grid to be adjacent to every other raster I got right now (each one was used in some analysis and each one is in the same SAGA's grid system).

Because spatial extent of vector data isn't allways the same it is quite problematic to set it manually using [-USER_XMIN <str>] [-USER_XMAX <str>] [-USER_YMIN <str>] [-USER_YMAX <str>] [-USER_SIZE <str>] [-USER_FIT <str>] options.

When using SAGA's GUI there is a wizzard to do such thing (to choose Grid system when any grid is loaded in SAGA) but my question is how to do it using batch mode?


texts: How to programmatically Pause the drawing of the Map?

texts: In ArcMap, there is a button to Pause the map's drawing. How do I programmatically call it from ArcObjects Code, before and after my heavy code runs?


texts:
I have a Feature Layer that I would like to buffer using user input and a geometry service.

FeatureLayer:

    var texasPipeline = new FeatureLayer(pipeURL, {
        mode: FeatureLayer.MODE_ONDEMAND,
        outFields: ["*"],
        definitionExpression:texasPipeQuery
    });

BufferParameters:

        var params = new BufferParameters();
        params.distances = [distance];
        params.unit = units;
        params.outSpatialReference = map.spatialReference;
        params.geometries = texasPipeline;
        map.graphics.clear();
        geomSvc.buffer(params, showBuffer);

The server is returning an error saying that geometries must be supplied. My guess here is that i need to pass in the geometry of the FeatureLayer as opposed to the FeatureLayer itself. How do i get at the geometries of the FeatureLayer and appropriately pass that into the BufferParameters??

EDIT:
Additionally I have tried to loop through as you can see in the code bellow. passing the array of geometries into the BufferParameters still does not return successfully.

        var texasPipelineGeom = [];
        var graphics = texasPipeline.graphics;
        for (var G in graphics) {
            var g = graphics[G]["geometry"];
            console.log(g);
            texasPipelineGeom.push(g);
        }


texts: Creating Perpendicular lines between groundwater contours and flow direction

texts: I have groundwater contours and flow direction arrows, each in their own feature class. I'm wondering if there is a way to draw the flow direction arrows so that they cross over the groundwater contours at a perpendicular angle?


texts: GEE - Add image band (NDVI) from one ImageCollection to corresponding image in an other ImageCollection

texts:
I have two ImageCollections one from Sentinel-1 and the other from Sentinel-2. From S-2 i derive the NDVI. I want to add the NDVI band from the S2 Collection to the Sentinel-1 ImageCollection with the closest date. How can i do this?

GEE basic script


texts: Connecting ArcGIS Desktop to remote PostgreSQL database?

texts:
I have setup PostgreSQL 9.5 on my local machine. I am trying to create a database connection from a coworkers laptop on the office network (let's say the network IP is 192.168.1.X). My database server is currently just called localhost (on port 5432). I was able to connect it to ArcGIS Desktop10.5 on the same machine using the add Database Connection option in ArcCatalog, no problem. So, I want to do the same on my coworkers computer. When I try, I get an "Unable to make database connection..." error.

I have configured my pg_hba.conf file multiple times to make sure that any computer on the networks (either 192.168.1.1 - the wifi or 10.1.10.1 - the ethernet) can access it  (with database authentication of course). I followed these ESRI configuration instructions.

Is there something I am missing?

Is there something else I need to do to make the database connection to postgres?

Can anyone see a problem with my pg_hba.conf file or with using "localhost" as the database instance?


texts: How do I transfer data from a geopackage to MS SQL server

texts:
I am trying to upload a geopackage (created with ogr2ogr though I don't think that matters) to a MS SQL Server 2008 instance.

I can't access the instance with ogrinfo all I can get is

C:\>ogrinfo -so "MSSQL:server=localhost\EXACTRAKDEVA;database=VectorMapData_2020-02-17;uid=sa;pwd=password"
ERROR 1: Unable to initialize connection to the server for MSSQL:server=localhost\EXACTRAKDEVA;database=VectorMapData_2020-02-17;uid=sa;pwd=password,

Try specifying the driver in the connection string from the list of available drivers:
SQL Server
SQL Server Native Client 10.0
SQL Server Native Client 11.0
ODBC Driver 11 for SQL Server
ODBC Driver 13 for SQL Server
SQL Server Native Client RDA 11.0

ERROR 1: Unable to initialize connection to the server for MSSQL:server=localhost\EXACTRAKDEVA;database=VectorMapData_2020-02-17;uid=sa;pwd=password,

I can't find any clues in the OGR manual page on how to specify a driver in the connection string. I'm also guessing how to manage to connect without  a trusted connection, though I have tried user, username as well as uid. I also am not sure what my "server" name is in.

Is there any (easy) way to generate a working connection string using the information available in the SQLServer connection params window?

If that's too hard I'll take a simple explanation of what to put in the QGIS connection boxes instead.


texts: How to show one to many relationships in GML from GeoServer/PostGIS?

texts:
I am just figuring out how to use GeoServer and PostGIS.  I am using the States.shp example that came with GeoServer, and I have uploaded that into a table into the PostGIS database.  I have created a new layer, and can make WFS requests to get back GML.  For example:

http://localhost:8080/geoserver/mine/ows?service=WFS&version=1.0.0&request=GetFeature&typeName= mine:usa_view&maxFeatures=1&

This returns a GML of the top 1 state (which happens to be Illinois).  I use a view, and aside from the geometric data, I also get all the fields I have exposed in the view, which happen to be “name” and “abbreviation”.

If I add a second table, say a list of cities within each state, is it possible to return this information within the GML?  Is there an easy way to do this using the setup I have described?


texts: v.surf.bspline and Interpolation(cublic spline) persistent errors in QGIS 3.2 & 3.3

texts:
I'm new to QGIS and as such I'm pretty hesitant to post anything, but I cannot for the life of me figure out why two raster creation techniques are not working.

For the v.surf.bspline GRASS algorithm: I keep getting 'ERROR: Choose either vector or raster output, not both' as the first error in the log.  No matter what I do, skip outputs, change parameters, etc, this happens.

For the SAGA Interpolate(cubic spline) algorithm, I keep getting '###' is not a correct numeric value for option 'K'. 140.0 is the default, but no matter the number, I get this result.

Have people successfully used these algorithms in 3.2 and/or 3.3?  I've tried both of these in 3.2 and the 3.3 dev release.  Failure.  In 2.18.22 the SAGA algorithm works fine, but I can't get any GRASS tool to work, although that's not the point of this post.  Here's a LINK to the data I've primarily been using (UDawsThick is the column with values to interpolate), however I have tried both algorithms with other point data, same result.

Apologies if I've missed something obvious. My coding skills are limited and there's every possibility that the answer was staring me in the face and I didn't have the lexicon to understand it.  With the limited documentation available for these tools, however, I'm not seeing another option other than appealing to you fine folk.

EDIT:
The log for the GRASS algorithm should now be in the link.


texts: How can I stretch the colors in my vizParams using EarthEngine API?

texts:
I'm displaying images on a map :
var landsat = [];
//load all the landsat years
for (var year = 2005; year < 2016; year++) {
  var start = year.toString() + '-01-01';
  var end = year.toString() + '-12-31';

  var dataset = ee.ImageCollection('LANDSAT/LE07/C01/T1_SR')
                  .filterDate(start,end)
                  .map(cloudMaskL457);

  var clip = dataset.median().clip(buffGeo)

  landsat[year.toString()] = clip;

  //Map.addLayer(clip, visParams, 'landsat-' + year.toString())
}

Map.addLayer(landsat['2010'], visParams, 'landsat')

To do so I'm using the following vizParams :
var params = landsat['2010'].reduceRegion({
    reducer: ee.Reducer.percentile([0, 100]),
    geometry: buffer,
    scale: Map.getScale(),
  });
print(params)

//vizparams for landsat
var visParams = {
  bands: ['B3', 'B2', 'B1'],
  min: [
    params.get('B3_p0').getInfo(),
    params.get('B2_p0').getInfo(),
    params.get('B1_p0').getInfo()
  ],
  max: [
    params.get('B3_p100').getInfo(),
    params.get('B2_p100').getInfo(),
    params.get('B1_p100').getInfo()
  ],
};

Map.addLayer(landsat['2010'], visParams, 'landsat stretch')

But I obtain a different stretch than when I use the button on the interactive map of the JavaScript code editor. Is it possible to reproduce this behaviour ?


texts: Statistically calculating real road from set of GPS tracks?

texts:
I participate in a long-distance cycling club, and we started to collect GPS data routinely from our riders.

My interest is to calculate "the real trajectory" for future events based on accumulated GPS data over the same roads. Basically, this would mean to pass some pre-selected tracks to an algorithm, and the algorithm would generate points at an appropriate sample rate (an appropriate distance from one another depending on road curves). I will discard timestamps, taking only spatial track information into account.

Which algorithm/statistic methods could I use? I don't use any GIS package and I plan to implement this in Python.

Below, some sample trajectory sets:


texts: How can i Calculate Road Width in ArcMap?

texts:
I have Road shapefile which is in polygon feature. I need to calculate the width of each road.

Can it be done in ArcMap?


texts: Where can I get free New Zealand and Melbourne geodata?

texts:
I have just installed QGIS. Do any forum users know where I could download or access files and maps for 1) New Zealand, and/or 2) Melbourne/Victoria, Australia. Preferably it would be good to know where free material is.

Thanks.


texts:
With QGIS 2.18, I used to capture the composer closure with PyQGIS, the class QgsComposerView and its composerViewHide signal:
self.obj_compo.composerView.composerViewHide.connect(self.raiseModule)

Now, there is the QGIS 3 answer:

AttributeError: 'composerClass' object has no attribute 'composerView'

I guess this is because "composer related methods have been removed from the public API and Python bindings. These classes have been replaced with the new layouts engine, based on QgsLayout, QgsLayoutItem, and the other related classes."
But I can't find a signal in QgsLayout class or a parent to capture the composer closure. If someone has any clues...


texts: How can I create an inset map in DDP that zooms to a feature in a layer besides the index layer?

texts:
(ArcMap 10.6.1)

I have a map that is using DDP to index through and show ~1800 short line features (bridges), and I need the map to be pretty zoomed in because I'm wanting to show the immediate context (especially local road names) of where the bridge is. While this is easy enough to do, in this application, I think it is important to give the broader context of the location. These bridges cover a 21 county area, and I want to have an inset map that shows the current county with an extent indicator that shows where in the county the current extent is showing. I realize that I could do this by having a separate map document for each county, but that feels like it would kind of defeat the purpose of DDP.

From what I've found researching this issue, it isn't possible to have multiple indexes, which is what I would want. Is there a way in python to make the inset map "zoom to feature" based on the county name that is listed for the current feature in the index?


texts: Using Python if/then/else to return value from either of two fields based on value of third using ArcGIS Field Calculator?

texts:
I'm trying to avoid manually inputting Begin and End engineering station values already inputted into a seperate field.

I'm dealing with(4) fields:

Begin_Sta
End_Sta
Beg_End and
Station.

The Begin_Sta and End_Sta fields have engineering stations.  The Beg_End field has either "BEGIN" or "END" inputted. The Station field is the one I want to populate using python in the field calculator.

If Beg_End has "BEGIN" then the engineering station value in the Begin_Sta field is returned into the Station field, else the End_Sta value is used to populate the Station field.

Is this possible using python in the Field Calculator?


texts: Getting all services for different years in GeoNetwork

texts:
I have some referential data services like country.
I have different versions of my country service for each year (geometry or attributs update) (country_2020 ; country_2012 ; country_2010 ; ...). In french, the word is "Millésime".
Is there a way with GeoNetwork to arrange all my country version/Millésime service and get all version/Millésime available service with a request ?


texts: Losing features when merging vector layers in QGIS

texts:
Trying the v.clean command and checked the layer with Topology checker and it turned out every single geometry was invalid.

What to do now?

Original question:

I have two shapefiles with 3D building blocks that I want to merge (Both MultiPolygonZ-type, the same columns in attribute table with the same projection). The blocks represent houses in two different area of the city (no overlaps).

Vector > Data Management Tools > Merge vector layers works great except that some features in the merged file looses it's solidness, see example below:

Houses from one of the shapes before merge (all houses are solid):

The same canvas from the merged file (some features loose their solidness):

Has anyone experienced the same problem and been able to workaround it?

Links below to two of the shape-files I'm trying to merge if anyone want to try and re-create the problem:

Dataset A: https://drive.google.com/open?id=0BwFHJFrNX33hbGFHS2ZCeUVJNWc

Dataset B: https://drive.google.com/open?id=0BwFHJFrNX33hZGFoT2xiVXZSbGs


texts: Matching legend symbology with data driven symbology?

texts:
I am tightly bound to a client's CI and fortunately, the data I am working with comes with RGB values in the attribute table. Perfect for data-driven symbology. But how can I (easily) manage to make my legend match with those?

Now I know about rule-based symbology but as stated in this answer to Changing marker shape based on a rule or column value in QGIS? that requires a lot of tedious clicking.


texts: Extract mean values of climate variables by level 2 administrative area country by country using R?

texts:
[Working with R version 3.2.2 in a Mac computer.]

I am trying to calculate mean values for climatic variables (http://www.worldclim.org/) for spatial objects in R. The challenge is that I'm trying to calculate these means for every level 2 administrative area in the world (http://www.gadm.org/), and I need an efficient way of calculating the statistics.

At the moment, I am trying to get this information country by country. Notice that in order to save time, I am just trying to do it for Austria (later on I will create a loop to do it for every country). Before starting I had to download the SRTM data and part of the script from here: https://github.com/sikli/srtm_country.

I think I miss some command to tell R to join tiles, boundaries, and climate data when I try to extract mean values for the climate variables.

In what follows, I present the procedure I am following.

#Clear all
rm(list = ls())

library(raster)
library(rgeos)
library(rasterVis)

setwd('/.../srtm_country-master')

#------------SETTINGS--------------

#Specify target ISO country code and path to downloaded shapefile

country_name <- "AUT"                             #Austria
shp          <- shapefile("/.../srtm_country-master/srtm/tiles.shp")       #Path to SRTM Tiles (can be found in subfolder srtm)

#------------EXECUTE FROM HERE--------------

#Get country geometry first
country <- getData("GADM",
                   country = country_name,
                   level=2)

#Intersect country geometry with tile grid
intersects <- gIntersects(country, shp, byid=T)
tiles      <- shp[intersects[,1],]

#Download tiles
srtm_list  <- list()
for(i in 1:length(tiles)) {
  lon <- extent(tiles[i,])[1]  + (extent(tiles[i,])[2] - extent(tiles[i,])[1]) / 2
  lat <- extent(tiles[i,])[3]  + (extent(tiles[i,])[4] - extent(tiles[i,])[3]) / 2

  tile <- getData('worldclim', var='bio', res=10,
                  lon=lon,
                  lat=lat)

  srtm_list[[i]] <- tile
}

THIS IS WHAT IT IS NOT WORKING -->

#Create a list with all exported .bil files from the climate data folder
ras_lst <- list.files("/.../wc10/",full.names=TRUE, pattern=".bil$")

#Load raster mosaic into R
climatemosaic <- stack(ras_lst)

#Extract mean values
AFGmeans <- extract(climatemosaic, country, fun=mean, na.rm=TRUE, small=TRUE, layer=1, nl=19, sp=TRUE)
AFGmeansdf <- as.data.frame(AFGmeans)

setwd("/.../srtm_country-master/Exports")
write.csv(AFGmeansdf, file= "AUTmeansdf.csv")


texts: How to segment polygon in to line segment?

texts:
I want to segment a Polygon or PolyLine in segment of line and remove some segments. I have an polygon that i convert in polyline by the code below and it work well:

from shapely.geometry import mapping, LineString
import fiona as fn
import operator
from functools import reduce

with fn.open('multipolygon.shp', 'r') as source:
    source.profile['schema']['geometry'] = 'LineString'
    with fn.open('segment_line.shp', 'w', **source.meta) as sink:
        for feat in source:
            xcoords_poly = feat['geometry']['coordinates']
            liste_coords = reduce(operator.concat, xcoords_poly)
            lines = LineString(liste_coords)
            assert lines.is_valid
            assert lines.geom_type == 'LineString'
            feat['geometry'] = mapping(lines)
            sink.write(feat)

After this step i want to remove some segments (lines) of my polyline where they overlays with some polygons as this picture show it :

For doing that i make the code below :

# great function that i find here to get coords of geometry
def pair(list):
'''Iterate over pairs in a list -> pair of points '''
for i in range(1, len(list)):
    yield list[i-1], list[i]

to_del = []
fn_obs = fn.open('obs.shp')
with fn.open('segment_line.shp', 'r') as source:
    with fn.open('avoid_obstacle.shp', 'w', **source.meta) as sink:
        for f in source:
            geom_f = shape(f['geometry'])
            try:
                for ob in fn_obs:
                    xob = shape(ob['geometry'])
                    if(geom_f.intersects(xob)):
                        to_keep = []

                        for seg_start, seg_end in pair(geom_f.coords):
                            line_start = Point(seg_start)
                            line_end = Point(seg_end)
                            segment = LineString([line_start.coords[0], line_end.coords[0]])
                            if(segment.intersects(xob)):
                                to_del.append(segment)
                            else:
                                to_keep.append(segment)
                        new_seg = linemerge(to_keep)
                        f['geometry'] = mapping(new_seg)
                sink.write(f)

            except (Exception):
                logging.exception("Error type of file:", f)

And i got this results :

It's like good but what i want it to remove all the line as this picture show it :