texts: Reprojecting with ogr2ogr -t_srs does not write AUTHORITY nodes to .prj

texts:
(Note: This seems to be the same thing that Mohayemin ran into at  "GDAL does not recognize EPSG code of a shapefilefile reprojected by ogr2ogr", but the comment there that "GDAL is writing out EPSG:3857 incorrectly" doesn't apply to my case here.)

Is it me, or does gdalsrsinfo know and display the different EPSG codes for an SRS, but a shapefile reprojection using ogr2ogr -t_srs <targetSRS> <dst> <srs> does not write these EPSG codes to the shapefile's .prj file?

I have a shapefile shapes.shp that I presume to be in an EPSG:4647 projection.
gdalsrsinfo shapes.shp gives me

PROJ.4 : '+proj=tmerc +lat_0=0 +lon_0=9 +k=0.9996 +x_0=32500000 +y_0=0
          +ellps=GRS80 +units=m +no_defs '

OGC WKT :
PROJCS["ETRS89_UTM32",
    GEOGCS["GCS_ETRF_1989",
        DATUM["ETRF_1989",
            SPHEROID["GRS80",6378137.0,298.257222101]],
        PRIMEM["Greenwich",0.0],
        UNIT["Degree",0.0174532925199433]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["False_Easting",32500000.0],
    PARAMETER["False_Northing",0.0],
    PARAMETER["Central_Meridian",9.0],
    PARAMETER["Scale_Factor",0.9996],
    PARAMETER["Latitude_Of_Origin",0.0],
    UNIT["Meter",1.0]]

For comparison, gdalsrsinfo EPSG:4647 shows this definition:

PROJ.4 : '+proj=tmerc +lat_0=0 +lon_0=9 +k=0.9996 +x_0=32500000 +y_0=0
          +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs '

OGC WKT :
PROJCS["ETRS89 / UTM zone 32N (zE-N)",
    GEOGCS["ETRS89",
        DATUM["European_Terrestrial_Reference_System_1989",
            SPHEROID["GRS 1980",6378137,298.257222101,
                AUTHORITY["EPSG","7019"]],
            TOWGS84[0,0,0,0,0,0,0],
            AUTHORITY["EPSG","6258"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4258"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",9],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",32500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH],
    AUTHORITY["EPSG","4647"]]

So it's basically identical, apart from the AUTHORITY nodes and the TOWGS84[0,0,0,0,0,0,0], right?

But since I'd like to have these AUTHORITY nodes in there to make it easier to identify the projection later, I call ogr2ogr -t_srs EPSG:4647 shapes_4647.shp shapes.shp. This results in (a gazillion warnings about feature fields not being convertible due to type width, which I ignore as spurious, and ...) a new shapefile called shapes_4647.shp, with this projection:

PROJ.4 : '+proj=tmerc +lat_0=0 +lon_0=9 +k=0.9996 +x_0=32500000 +y_0=0
          +ellps=GRS80 +units=m +no_defs '

OGC WKT :
PROJCS["ETRS89_UTM_zone_32N_zE_N",
    GEOGCS["GCS_ETRS_1989",
        DATUM["European_Terrestrial_Reference_System_1989",
            SPHEROID["GRS_1980",6378137,298.257222101]],
        PRIMEM["Greenwich",0],
        UNIT["Degree",0.017453292519943295]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",9],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",32500000],
    PARAMETER["false_northing",0],
    UNIT["Meter",1]]

Note: no AUTHORITY nodes anywhere.

I thought that including the EPSG codes in the target SRS would be trivial, given that they are obviously on file somewhere...

Should I just manually overwrite shapes_4647.prj with the "official" SRS WKT returned by gdalsrsinfo EPSG:4647?


texts: Line Intersections tool missing from QGIS 2.16

texts: In QGIS 2.14 the Line Intersections tool was under Vector->Analysis. I can't find it in 2.16 and it doesn't appear to have been integrated as an option in the normal intersection tool. Google hasn't been helpful. Does anyone know where it's been moved too? Related question, does QGIS have a tool search function the same way esri does?


texts: Intersecting two polygons with same FID?

texts:
I need to perform an intersection between two layers, one is blocks and the other is buffers for each block which are of the same size of the road. Blocks and buffers share the same FID.
What I'm trying to obtain is the compactness of each block. I tried using an iterator but it either takes forever or after hours of processing gives me only the intersection for one pair. Also tried tabulate intersect but the output. A regular intersection doesn't work because some buffers go on neighbouring blocks and I only want the intersection between polygons sharing the same FID.
My python skills are limited but I'm learning.


texts: Field Calculator Only Working for Beginning and End of Data Set

texts:
I'm working with a dataset where every entry is a unique zip code in the United States (approx. 33,000 entries). I ran a simple field calculation:
def Rank(minrate):
    if minrate > .4:
      return 'HiMin'
    else:
      return 'LoMin'

and it worked perfectly for all 33,000 entries. I'm running a very similar calculation:
def Rank(zip,state):
    if zip > state:
      return 'HiLoan'
    else:
      return 'LoLoan'

This is an example of the data. "zip" refers to the "loanpen" field and "state" refers to the "statepen" field. Both are double fields, with the LoanRank field being text. You can see AK working correctly above, with CA failing below. (A few natural nulls in the AK data that are correct).
And I'm running into an issue I've never seen before. The calculation appears to be working fine, if you look at Arizona, Arkansas, Alaska, Alabama, or start from the bottom and look at Wisconsin, Wyoming, and West Virginia... But for every other state in the middle (that doesn't start with A or W), the calculation is failing. It's not giving any error, just populating as Null.
I've tried doing a definition query on one of the failed states and running the code just on that state - still fails. I tried running the code on a single cell within one of the failed states - failed. I tried creating new columns, still fails.
I tend to think the issue is not with the Python, but something with ArcGIS.
If anyone else has experienced this, do you have an idea on what the issue might be or how to fix it?

Here's the data sorted by UID in the natural order, showing the null values do NOT break the code.


texts: I want to upgrade my QGIS version to the newest one, Lisboa, but I'm not sure if I should uninstall Capiapo first. Is there a way to upgrade without uninstalling?


texts: Calculate Raster mean with gdal_calc - unexpected result

texts:
I have in input a single band 2 axis tif raster (float32 values). Using gdal_calc I need to produce an output raster (same dimensions of input) that contains in every pixel the mean value of the whole input raster (to be used in other calculation) - basically a "uniform" raster. To do this I use:

gdal_calc.py    --calc="mean(A,dtype=float64)"  --type="Float32"    --overwrite --NoDataValue=-32768    --outfile   output.tif  -A  input.tif

Based on

gdalinfo -stats input.tif

I get a mean of 0.45 approx. So I was expecting an output.tif with this value for every pixel, but what I get is a different image (constant horizontally but different values vertically). I tried several numpy.mean parameters  (axis, keepdims) but nothing changed.

Any suggestion or work around?


texts: Find the first occurrence of any letter in an alphanumeric string

texts:
Working in arcgis desktop 10.2, I use python scripts to split some strings, found in an attribute table. Most of them were quite easy to split, but here is one that i could't manage to break.

string_example = "123a456"

I need to find were the first letter appears. Knowing it's position I will test if everything left from that letter are numbers AND everything right from that letter are numbers. If so i will write each part of the string in a new field.

FieldA = "123a456"

find the first letter, test if it is surrounded by numbers

FieldB= "123"

FieldC= "a"

FieldD= "456"

There are variations of the string, plenty of them: "12a1","12A1", "123ab1234", "a123b456". I need to find the structure: number-letter-number.


texts: Adding an EPSG:2056 WMTS layer to Cesium

texts:
I can implement EPSG:4326 WMTS layers to my Cesium App like this:

const imageryProvider = new this.Cesium.WebMapTileServiceImageryProvider({
        url: "https://openwms.domain.ch/service/wmts",
        layer: "layerName",
        style: "",
        tileMatrixSetID: "EPSG:4326:256",
        tileMatrixLabels: this.getTileMatrixLabels(),
        maximumLevel: 19,
        tilingScheme: new this.Cesium.GeographicTilingScheme()
});

const imageryLayer = new this.Cesium.ImageryLayer(imageryProvider);
scene.imageryLayers.add(imageryLayer);

But my normal WMTS layers are EPSG:2056. I could only find the two TilingSchemes:

GeographicTilingScheme (EPSG:4326)
WebMercatorTilingScheme (EPSG:3857)

How can I get my EPSG:2056 WMTS layers in Cesium? The options I could think of:

Save all my layers (3k+) in EPSG:2056 and EPSG:4326
Create an own TilingScheme for EPSG:2056.

How could I create an own TilingScheme? Does this make sense? Or is there a better way to deal with this problem?


texts: Trying to measure point density along a line

texts: I have a dataset that consists of animal sightings along a road and I want to calculate the density of points (sightings) along that road.  My primary goal is to measure sighting hot spots, but I also want to subsequently look at the habitat characteristics associated with these hotspots.  Kernel and point density do not seem appropriate b/c I only have measurements along this road (a line) I have no measurements of occurrences at locations away from this line. I'm using ArcGIS 10.3  Suggestions?


texts:
I'm just starting with ArcObjects and I want to update coordinates of geometry vertices. Specifically, I want to swap X, Y coordinates for each vertex of polygons (which can be multipart and have inner rings).

I've found web help topic Updating geometry of existing features where ITransform2D interface is suggested. But it's methods does not seem solve my problem directly. I think I can use Move method, but it requires calculating dx and dy.

Is it possible to do coordinate replacement in ArcObjects like in Python vertex.X = vertex.Y but without recreating geometry? (this is quite tricky in Python for polygons with inner rings)

I adapted @vinayan code, but it does nothing. What do I miss?

protected override void OnClick()
{
    IMxDocument mxDocument = ((IMxDocument)(ArcMap.Application.Document)); // Explicit Cast
    IContentsView currentContentsView = mxDocument.CurrentContentsView;
    IFeatureLayer featureLayer = (IFeatureLayer)currentContentsView.SelectedItem; // Explicit Cast
    IFeatureClass featureClass = featureLayer.FeatureClass;

    //Get the Feature or FeatureCursor as you like..
    IFeatureClass pFC = featureClass;
    IFeature polygonFeature = pFC.GetFeature(1);

    IPointCollection pPtsColl = (IPointCollection)polygonFeature.ShapeCopy;

    //swap X and Y
    double tempValue;
    IPoint pPoint;

    for (int i = 0; i < pPtsColl.PointCount; i++)
    {
        pPoint = pPtsColl.get_Point(i);
        tempValue = pPoint.X;
        pPoint.X = pPoint.Y;
        pPoint.Y = tempValue;
    }

    //Update Feature
    polygonFeature.Shape = (IGeometry)pPtsColl;
    polygonFeature.Store();
}


texts: automated selection by attributes and location within another Layer

texts:
I have two Polygon Layers. Layer A are Districts of a City and Layer B are Blocks of a City. Currently i'm working with QGIS 2.6.1 and i don't have  access to any other GIS.

What i'm trying to do is: select Blocks in Layer B within (contain) a specific District in Layer A. I'm already able to select a certain District in Layer A by using "select by Attributes" followed by "select by Location" with Layer B by checking the "current selection" Option.

These Kind of Operation i've to do quite often and i'm looking for a way now where i can build a automated process with the "graphical modeler". So far the single steps are working but i'm not able to choose a specific Field and Value (District; Distict Name) i.e. a dropdown menue of Layer A. So even with that Modeler i've to fill out the Field and Value for each Process.

Is there any Way for an automated Process?


texts: How do you programmatically read the coordinate display format in ArcMap?

texts:
I'm looking for a property (or properties) that will give me the display units and the formatting.

I assume I need IMap.DistanceUnits but I can't find the formatting.  Is this exposed?

I am also interested in accessing the additional formatting in this window:

I'm hoping there's a INumberFormat property someplace that I can hook into.


texts:
I'm trying to make a filter function that I add to an URL.
I managed to make it work for polygon, but I'm struggling with filter function for point type.

Currently I am trying to use something like this:

<Filter xmlns="http://www.opengis.net/ogc/" xmlns:gml="http://www.opengis.net/gml"><Contains><PropertyName>geometry</PropertyName><gml:Point><gml:coordinates>520757.13197114784 120868.02826434454 </gml:coordinates></gml:Point></Contains></Filter>

I keep getting <ows:ExceptionText>java.util.NoSuchElementException</ows:ExceptionText>, I narrowed my problem to filter, since I get correct response without it.


texts: Processing Script: now gives error in QGIS 2.2

texts:
I recently updated to QGIS version 2.2, and am now getting an error when running Processing scripts that worked previously (in version 2.0.1).

Uncaught error while executing algorithm
Traceback (most recent call last):
Traceback (most recent call last):
File "C:/PROGRA~1/QGISVA~1/apps/qgis/./python/plugins\processing\core\GeoAlgorithm.py",    line 202, in execute
self.processAlgorithm(progress)
File "C:/PROGRA~1/QGISVA~1/apps/qgis/./python/plugins\processing\script\ScriptAlgorithm.py", line 252, in processAlgorithm
exec script in ns
File "<string>", line 25, in <module>
AttributeError: 'module' object has no attribute 'getobject'

I believe it's this line of code causing the error, but i'm not sure how to modify it to work with QGIS 2.2:

indivLayer = processing.getobject(Couche_compteurs_individuels)


texts: Using Adobe Illustrator files from CloudMade in ArcGIS for Desktop?

texts:
Anyone here able to enlighten me with your knowledge of Adobe Illustrator files?

I have spent several hours trying to find a way to view an Adobe Illustrator (.ai) file from CloudMade.

Is there a way (preferably free) to at least view these *.ai files in ArcGIS 10 (or any other GIS software)?

ArcGIS has had "export to .ai" for at least 10 years but there seems to be no way to use (or at least view) them in ArcGIS.


texts: Multispectral image segmentation for natural-resources applications using R

texts:
R has the capacity for image segmentation, although all of the examples I have come across utilize a single band for the segmentation (example).  I am interested in combining R's capacity for random forest image classification with an object-oriented segmentation approach.

What functionality does R have for multispectral image segmentation that is suitable for natural resources-based analysis?  Or how to link the results of a single band segmentation for further analysis.


texts: OpenGeo postgres installation not working from command line but works from PGadminIII on fresh install

texts:
Problem:

I just did a default install of opengeo suite (4.1) on Ubuntu 12.04 and followed instructions exactly.  The problem is that I can connect just fine to my postgresql through pgadminIII (1.18.1), but I can't connect on command line.  Has anyone encountered this and if so, how did you fix?

background:

As I said above, I am installing opengeo suite (4.1) on a 32bit ubuntu 12.04 installation.  In fact I did it twice, once at home and once at work and I get the same error so it is repeatable:

Basically, i follow all the instructions on the opengeo website. (http://suite.opengeo.org/opengeo-docs/installation/ubuntu/install.html)

I am sure I am using the 12.04 repository to download software.  I changed the postgres.conf and pg_hba.conf to listen on my machine and accept outside connections and use md5 authentication.

I created a password for the postgres user manually like it instructs.  However, when it comes time to test the command line, it doesn't authenticate. I get the error:

~$ psql -U postgres -W
Password for user postgres:
psql: FATAL:  Peer authentication failed for user "postgres"

The problem is, that I can connect just fine through pgadminIII (1.18.1).  I can run the postgis extension and it creates my postgis database.  I can even edit files I created in QGIS.  But I can't log in via command line so I can't run the shp2pgsql to load any preexisting data.  This is a fresh installation on a machine and, as I mentioned, I tried it on two separate machines (both running ubuntu 12.04 x32).

my pga_hb.conf file:

# TYPE  DATABASE        USER            ADDRESS                 METHOD

# "local" is for Unix domain socket connections only
local   all             all                                     peer
# IPv4 local connections:
host    all             all             192.168.0.0/16          md5
# IPv6 local connections:
host    all             all             ::1/128                 md5
# Allow replication connections from localhost, by a user with the
# replication privilege.
#local   replication     postgres                                peer
#host    replication     postgres        127.0.0.1/32            md5
#host    replication     postgres        ::1/128                 md5


texts:
I'm having a problem using tmaps to map out points in R. The data is taken from the Chicago data portal, and I was able to turn it into a sf object using the code below. However, when I try to use tmaps, it just returns one point in Missouri(?). I've checked the coordinates for each point online and they are correct (or at least in Chicago), and this same problem is happening for another similar sf object I created from affordable housing projects in Chicago. I've played around with a number of different tmap options, and plot() seems to work fine.

here is the code I used:

tif_data <- read.csv("~/Documents/code/ChicagoPackage/data/TifProjects.csv")
tif_data <- tif_data[!is.na(tif_data$LATITUDE), ]
tif_sf <- st_as_sf(tif_data, coords = c("LONGITUDE", "LATITUDE"), crs = 102671)
tif_sf <- st_transform(tif_sf, crs = 4326)
tm_shape(tif_sf) + tm_dots()

Here is the website where I got the data, and below is a screenshot of where the point is.


texts: How to alter labels in Mapbox Streets in Mapbox Studio

texts:
I am creating vector tiles in Mapbox Studio v5 based off of the MapBox Streets layer.  I want to alter the labels that appear over buildings, parks, schools, churches etc. but I cannot find the code that deals with these labels.  I found the housenum_label section of the code but that does not deal with building or area names.  How can I change the style of these labels?

Here's an example of the labels I want to alter:


texts: Measuring different values in latlong Geotiff whether OTFR enabled or not

texts:
I'm using a GeoTiff file in QGis that has no projection per-se (it uses lat-long coordinates) and WGS84 datum. Further it has an x pixel size of 2.40475e-07 degrees and a y pixel size of -1.99028e-07... so multiplying by the radius of the earth where it is standing according to the WGS84 geoid (6371.30 km) I get a pixel size of 2.6 cm in x and 2.2 cm in y.

Importing directly the file into QGis and measuring the size of a pixel with the Measure Line tool i get exactly those values... so far so good.

But if I enable OTFR in any CRS (inluding the same used by the file, namely lat-long WGS84) when i measure the size of the pixel y get a size of 2.2 in x and 2.2 in y.

Is the pixel size (in degrees) given by the Geotiff as if it were in the equator?


texts:
using QGIS 1.8.0 lisboa on windows 7 pc

adding vector layers and they are not adding properly they don't show up in the layers list but they are there when I happen to merge the the layers altogether, I can see them listed, they are not hidden.

it seems to happen when i try to add layers and have control layer render order on, because when i turn it off the layers adds fine, but i cannot delete the previous layer that didn't add correctly, i can't see them to select and remove them


texts:
I have a netCDF file and would like to convert it to LAT-LON coordinates.
1) I do not know the projection in the file but from the website of the data provider, the data is on a polar-stereographic (PS) grid covering North America and adjacent waters with a 10 km resolution at 60 degrees north.
2) How can I project to latlon using a suitable coordinate system for the north pole?

Here is my sample file:

https://www.dropbox.com/s/lgni1s5tgu1sbwh/data.nc?dl=0


texts:
I have a certain geographic region defined by the bottom left and top right coordinates. How can I divide this region into areas of 20x20km. I mean in practial the shape of the earth is not flat it's round. The bounding box is just an approximation. It's not even rectangular in actual sense. It's just an assumption. Lets say the bottomleft coordinate is given by x1,y1 and the topright coordinate is given by x2,y2, the length of x1 to x2 at y1 is different than that of the length between x1 to x2 at y2. How can I overcome this issue

Actually, I have to create a spatial meshgrid for this region using matlab's meshgrid function. So that the grids are of area 20x20km.

meshgrid(x1:deltaY:x2,y1:deltaX:y2)

As you can see I can have only one deltaX and one deltaY. I want to choose deltaX and deltaY such that the increments create grid of size 20x20km. However this deltaX and deltaY are supposed to vary based upon the location. Any suggestions?

I mean lets say deltaX=del1. Then distance between points (x1,y1) to (x1,y1+del1) is 20km. BUt when I measure the distance between points (x2,y1) to (x2, y1_del1) the distance is < 20km. The meshgrid function above does creates mesh. But the distances are not consistent. Any ideas how to overcome this issue?

For finding point p(1,0)

deg = km2deg(20)

deg =

    0.1799

>> reckon(25,-120,0.1799,-18)

ans =

   25.1711 -120.0614

Function to calculate distance between two points

function [ distance ] = calculateDistance( latitude1,longitude1,latitude2,longitude2 )
radius = 6371;
dLat = degtorad(latitude2-latitude1);
dLon = degtorad(longitude2-longitude1);
a = sin(dLat/2) * sin(dLat/2) + cos(degtorad(latitude1)) * cos(degtorad(latitude2)) * sin(dLon/2) * sin(dLon/2);
c = 2 * atan2(sqrt(a), sqrt(1-a));
distance = radius * c;
end

Referring to whuber's suggestions

syms p(1,1).lat p(1,1).lon
solve ( calculateDistance(p(0,1).lat,p(0,1).lon,p(1,1).lat,p(1,1).lon), calculateDistance(p(1,0).lat,p(1,0).lon,p(1,1).lat,p(1,1).lon))


texts: Retrieve OSM changeset information (pyosmium)

texts:
I would like to ask if there's an easy way to query OSM changesets' information like user, timestamp, comment and hashtags in that changeset.
In particular I would like to retrieve all the user names of the users which modified an area saving their changesets with a particular hashtag (similarly to default hashtags in HOT projects https://tasks.hotosm.org/)

My idea is to download the OSM data relative to my area of interest (in this case a country, taken from https://osm-internal.download.geofabrik.de/africa/senegal-and-gambia.html#)

So what I did until now is to download the following files from GeoFabrik.de:

senegal-and-gambia-latest-internal.osm.pbf
senegal-and-gambia-internal.osh.pbf

And I tried to run the following code (taken from https://oslandia.com/en/?s=timeline+handler) using Python and pyosmium library:

import osmium as osm
import pandas as pd

class TimelineHandler(osm.SimpleHandler):
    def __init__(self):
        osm.SimpleHandler.__init__(self)
        self.elements = []

    def add_elements(self, e, elem_type):
        self.elements.append([elem_type,
                              e.id,
                              e.version,
                              e.visible,
                              pd.Timestamp(e.timestamp),
                              e.uid,
                              e.changeset]) #only changeset id

    def node(self, n):
        self.add_elements(n, 'node')

    def way(self, w):
        self.add_elements(w, 'way')

    def relation(self, r):
        self.add_elements(r, 'relation')

tlhandler = TimelineHandler()
tlhandler.apply_file("senegal-and-gambia-latest-internal.osm.pbf")
colnames = ['type', 'id', 'version', 'visible', 'ts', 'uid', 'chgset']
elements = pd.DataFrame(tlhandler.elements, columns=colnames)

This doesn't solve my problem as I'm searching for changesets' information and not OSM data itself.

Maybe I'm missing a snippet to be added in the handler similar to:

    ...

    def changeset(self, c):
        self.add_elements(c, 'changeset')

    ...

I tried to add this snippet, but the handler gives no result.

Actually I'm wondering if I'm missing something in the data structure of osm/osh .pbf files, something on pyosmium or if the changeset information I'm searching for is contained in such files.

If the task cannot be performed with pyosmium and GeoFabrik, how can I do such operation?


texts: Splitting all polygon features in a shapefile into smaller equal width polygons

texts:
I have a shapefile of thousands of polygon features, that i would like to break up into smaller pieces (i.e. break a square into 8 equal pieces). I have read all over these boards for the right solution and haven't found one yet unfortunately. I've tried parcel editor but that only allows you to manually edit one feature at a time, and I've tried fishnet but that doesn't run through each individual feature it just processes across the entire extent of the shapefile. I have the full license for ArcGIS so all options are game. Ultimately I would like to be able to split this initial shapefile of features into a specified amount of equal width slots that each would be their own feature in the resulting shapefile.

Any thoughts?


texts: Why does Lidar to Raster tool in ArcGIS create shifted raster cells?

texts:
I am testing the results of Lidar to Raster with 0.5m resolution with ArcGIS 10.3.1 SAGA-GIS and Grass 7. My LAS File consists of regular spaced points with a 0,5m resolution. No coordinate transformation took place. The resulting rasters are shown in the screenshots below. I created the  LAS Dataset in ArcGIS  by right clicking in the Catalog window in ArcGIS >New>Las Dataset. The Tool LAS-Dataset in Raster was used to produce the raster. The ArcGIS Raster cells are shifted and do not center on the point data from the lidar. The rasters created in Grass-GIS by running r.in.lidar > Export to raster or with SAGA-GIS using the workflow lidar to Pointcloud - Pointcloud to Grid > Export Grid to GTIF both create cells whose center is the lidar point.

Can anyone explain this difference and whether it has an negative effect on the validity of the resulting dataset?

ArcGIS Settings:

ArcGIS Raster:

Grass Raster:


texts: Geocoding address involving non-ASCII characters in PyQGIS

texts:
I have installed the "GeoCoding" Plugin in QGIS and it worked well even though I have used non-ASCII characters to geocode (i.e. Hacettepe Üniversitesi). However, when I tried to do the same in PyQGIS with the following code I got the error:

"UnicodeEncodeError: 'ascii' codec can't encode character"

The problem is that the returned address contains non-ASCII characters. Any ideas to solve this?
from geopy.geocoders import Nominatim

#Next, we'll create our geocoder:
geocoder = Nominatim()

#geocode an address:
location = geocoder.geocode("Hacettepe Üniversitesi")

print location

Another relevant question: I rely on Joel Lawhead's book ('QGIS Python Programming Cookbook') and the import statement relevant to this example in the book is:
from GeoCoding.geopy.geocoders import Nominatim

However, this statement gave the following error

ImportError: No module named geopy.geocoders

Is there a source we could follow to make these import statements without a try-and-guess?


texts: Resampling gives higher MEAN when using zonal statistics

texts:
I am using ArcGIS 10.1 for Desktop.

I have a value raster and a polygon feature in GCS coordinates. My raster resolution is 0.5 decimal degrees which is much larger than the size of my polygons. My value raster has discrete floating points. So I resampled my raster to smaller resolution of 0.01 using the NEAREST option and used Zonal Statistics in table to get the mean value for each polygon. However when I compare the means of the bigger polygons that have zonal statistics values for both original and resampled raster, I find that the resampled values are higher. My raster values are in terms of mm/sq km. So the small difference during resampling is amplified when I scale the value up with the polygon areas. Any ideas how to prevent this?

I have also tried doing the same with projecting both my raster and polygons to WGS and then resampling and zonal statistics. But i get the same trends in the results. Am I missing something?

Following are 2 graphs showing:
1) comparison of regional zonal means evaluated using original GCS raster and evaluated using resampled/projected+resampled rasters
2) comparison of country level values (zonal mean*area) and sum of all district values

It seems that the country level and sum of district level values match each other. But both are generally higher for resampled and projected rasters. So resampling is giving me higher zonal means.


texts: Using Generate Near Table and keep the ID attributes

texts:
I try to use near table between two datasets to calculate the nearest distance. The original data is as the right one, and the near table generated is as the left one. The problem is the attribute of "IN_FID" and "NEAR_FID", none of them keep the correct values from the data. The ID number is important to me, how can I keep it? I have tried to join the attribute, but it turns out "Null" values.

Can anyone help me?


texts:
Is there any way to style the WMS layer using client side script
I use js to style WMS layer. Code is shown below....
but it's not work.

/*****************************************************************************/
//applying styles for unescap layer

var unescap = new ol.layer.Tile(
    {
        source: new ol.source.TileWMS({url: 'http://203.159.29.11:8200/geoserver/unescap/wms?',params: {LAYERS: 'unescap:rg_di_any_all'}}),
        style: new ol.style.Style({
                                    stroke: new ol.style.Stroke({color: '#319FD3',width: 1}),
                                    fill: new ol.style.Fill({color: 'rgba(17, 224, 36, 0.66)'}),
                                    text: new ol.style.Text({
                                                              font: '12px Calibri,sans-serif',
                                                              fill: new ol.style.Fill({color: '#ff0'}),
                                                              stroke: new ol.style.Stroke({ color: '#09C',width: 3})
                                                            })
                                })
    });


texts: Getting errors after about 30 mins of processing Flow Direction using ArcGIS Desktop?

texts:
I'm trying to process the flow direction from a filled DEM using ArcGIS 10.2.2 and I get these errors after about 30 mins of processing, every time, no matter what I try. I tried changing the raster's format (form GRID to TIFF), the projection, the Drive it's save on, etc.

ERROR 010429: Error in GRID IO: S_VXfer: Read error.
ERROR 010213: Error in reading raster [path].
ERROR 010429: Error in GRID IO: S_VXfer: Write error.
ERROR 010429: Error in GRID IO: S_VXfer: Write error. 2- Another person or application is accessing this directory
ERROR 010067: Error in executing grid expression.
Failed to executed (FlowDirection).

Can anyone help?

EDIT : It works when I reduce the processing extent to calculated the flow direction on a small part of my DEM, but it won't work for the whole layer.


texts: Deactivating InfoWindow on certain Elements using CartoDB?

texts:
In CartoDB: How can I disable/deactivate the InfoWindow on certain Elements? For example if I have a layer with multiple lines, but only want a on-click-Infowindow or on-hover-InfoWindow on line 1 and 3, but not 2?

Ditto for Labels.

example: show infowindow on the solid red lines, but not on the dashed line (all lines are on the same layer)


texts: Shared markers/paths in layer groups - toggle issues

texts:
I'm working on a map that displays the trips I've taken. I've created many markers and lines to designate locations and the paths between them.  I've created trips using the required markers/paths and added them as layerGroups to a layerControl.  The problem is several locations and drive paths have been used in multiple trips.  I'd like to have multiple layers render but the problem is unchecking one layer will turn off markers/paths that are common with other still active layers.

For example:

Trip 1 = Burlington, VT marker -> BurlingtonToAlbany path -> Albany, NY marker ->AlbanyToDC path -> Washington DC marker

Trip 2 = Boston, MA marker -> BostonToAlbany path -> Albany, NY marker -> AlbanyToSanDiego path

Both trip 1 and trip 2 display properly. However, when I turn off trip 1 all related markers/paths are removed. This includes the shared Albany, NY marker and Trip 2 is no longer rendered correctly.  How can I best resolve this?


texts: In a Label Tool is it possible to display value of a point, for example 3.2 (double type), with subscript after a dot, like that: 32 ?


texts: vb.net arcobjects validate topology function only works after I use arcmap topology tool validate topology first

texts:
I find a strange thing occasionally, I have a validate topology function developed by vb.net arcobjects. It work not properly. The case is:

open arcmap
start editing
use my validate topology function, nothing happen, no error find (there really several errors)
set the extent to global and use arcmap topology toolbar 'validate topology in current extent' tool to validate topology
open error inspector and get errors list (4 errors)
change the current extent and use arcmap topology toolbar 'validate topology in current extent' tool to validate topology
open error inspector and get errors list (2 errors)
remain the same extent as step 6, use validate topology function developed by vb.net arcobjects, all errors displayed in error inspector.

I want to know is it my code bug or operator issue?

My code:

    ' Step 1 - access a topology
Dim topoUiD As UID = New UIDClass
topoUiD.Value = "esriEditorExt.TopologyExtension"
Dim m_application As IApplication = My.ArcMap.Application
Dim topologyExt As ITopologyExtension = CType(m_application.FindExtensionByCLSID(topoUiD), ITopologyExtension)  'original code from esri
Dim topology As ITopology

topology = CType(topologyExt.CurrentTopology, ITopology)
MsgBox("successfully access current topology")

' Step 2 - Finding topology errors
Dim geoDS As IGeoDataset = CType(topology, IGeoDataset)
Dim errorContainer As IErrorFeatureContainer = CType(topology, IErrorFeatureContainer)
Dim eErrorFeat As IEnumTopologyErrorFeature
eErrorFeat = errorContainer.ErrorFeaturesByRuleType(geoDS.SpatialReference, esriTopologyRuleType.esriTRTLineCoveredByLineClass, Nothing, True, False)
Dim topoError As ITopologyErrorFeature

' check the eErrorFeat length to determine whether has error
topoError = eErrorFeat.Next

topologyExt.ClearActiveErrors(esriTEEventHint.esriTENone)
topologyExt.DelayEvents(True)

Do While Not (topoError Is Nothing)
    MsgBox("Topology error exists!")
    topologyExt.AddActiveError(topoError, esriTEEventHint.esriTENone)
    topoError = eErrorFeat.Next
Loop

topologyExt.DelayEvents(False)

Another thing is if I am not in editing session and use validate topology tool developed by vb.net arcobjects, the map will crash.


texts: How to use Qt3 Support widgets in QGIS plugin?

texts:
Similar to this post: My plugins crashes with “No module named qgsprojectionselectionwidget” in Windows.

Except I would like to also use Qt3 Support widgets (in particular, a Q3ButtonGroup). However, QGIS does not recognise it and provides the following error:

  NoSuchWidgetError: Unknown Qt widget: Q3GroupBox

To solve it with QGIS custom widgets, you can add qgis.gui to the ui file.

Is there a similar fix (or any fix) for Qt3 Support widgets?

Here is an example of my "dockwidget_base.ui" file:

  <customwidget>
   <class>Q3GroupBox</class>
   <extends>QGroupBox</extends>
   Qt3Support/Q3GroupBox
   <container>1</container>
  </customwidget>
  <customwidget>
   <class>Q3ButtonGroup</class>
   <extends>Q3GroupBox</extends>
   Qt3Support/Q3ButtonGroup
   <container>1</container>
  </customwidget>
  <customwidget>
   <class>QgsCollapsibleGroupBox</class>
   <extends>QGroupBox</extends>
   qgscollapsiblegroupbox.h
   <container>1</container>
  </customwidget>
  <customwidget>
   <class>QgsMapLayerComboBox</class>
   <extends>QComboBox</extends>
   qgsmaplayercombobox.h
  </customwidget>
 </customwidgets>
 <resources/>
 <connections/>
</ui>


texts: Filtering parameters in ArcPy Python script tool avoid repeating selecting same value in sequence parameters?

texts:
I am trying to make a tool in which there are many sequence parameters (based on the same field of a feature class), and each of them is multi value to make a user choose from a check a list one or more value from values in this field.

If the same field's values I need to use in sequence parameters, but I want to avoid the duplication in selecting values in these sequence parameters, so how to adjust the validation code, so that if the user select definite value(s) in a parameter, that or these values can not be selected in the following parameter, unless that or these values are unselected in the previous parameters?

For example, if in a parameter the user select the value "Canada" based on a field for countries' names, I want to avoid the repeating of selecting the value "Canada" in a following parameter, that is also based on the same field for countries' names, unless the user return to the previous parameter and unselect or uncheck  the value "Canada".

import arcpy
class ToolValidator(object):
    def __init__(self):
       """Setup arcpy and the list of tool parameters."""
       self.params = arcpy.GetParameterInfo()

    def initializeParameters(self):

        """Refine the properties of a tool's parameters.  This method is
        called when the tool is opened."""
        return

    def updateParameters(self):
        if not params[0].value:
            """Modify the values and properties of parameters before internal
            validation is performed.  This method is called whenever a parameter
            has been changed."""
            global previousProjectPath3
            previousProjectPath3 = ""
            if self.params[1].value and self.params[3].value:
                self.params[4].filter.list = sorted({row[0] for row in arcpy.da.SearchCursor(self.params[1].value, self.params[3].value.value) if row[0]})
            if self.params[1].value and self.params[3].value:
                self.params[5].filter.list = sorted({row[0] for row in arcpy.da.SearchCursor(self.params[1].value, self.params[3].value.value) if row[0]})

    def updateMessages(self):
        """Modify the messages created by internal validation for each tool
        parameter.  This method is called after internal validation."""
        if params[0].altered:
            path = params[0].valueAsText
            global previousProjectPath3
            if previousProjectPath3 != path:
                previousProjectPath3 = path


texts: Collecting a Boundary of an Area using ODK

texts:
I am trying to get a geo shape of a place here and I'm using ODK Collect to Map the area. I am able to create the Geo shape using ODK but the problem is, I am not able to export this data to QGIS for use and create shapefiles.

I do not know how I can do this or if there is another tool I can use.


texts: OSRM Nearest health facility to a village

texts:
I work in a health research project in rural West Africa. I want to measure distance by road from a village to closest health facility. I use Stata 16. Per now I only have distance "as the crow flies" using geonear. I have 258 villages.
I've downloaded a map from geofabric downloads (Guinea-Bissau) and followed instructions from https://journals.sagepub.com/doi/pdf/10.1177/1536867X1601600209 and https://sites.google.com/view/stephanhuber.
I am also able to select which health facility a village should link to. And I have latitude and longitude (DD) for all locations.

Is OSRM a good solution for my problem?
I also have trouble running it, and receives the code no route found, even though the map seems to be detailed enough.


texts: Add and change attribute / feature cascaded (QGIS)

texts:
I have a polygon layer with hierarchical polygons. The main polygon has a capital letter like 'A' or 'B'. All polygons inside the main polygon got a number and all inside them get a second, third, fourth... number. You can see it in the graphic below.
I want to create new polygons inside a other, they should get the letter and numbers of there parent features. And also if I change a parent, the child's (an maybe the "older" parents) should automatically update.
For example: in the second graphic, I changed the letter A to B. All related polygons should now update.

I did use an aggregate function to get the attributes of the parent feature. This works for adding new polygons but won't help if I change a parent. I think aggregate is not the perfect function since I'm working in just one layer:
aggregate('layer','concatenate_unique',to_string("parent"),within(centroid(geometry(@parent)),$geometry))

At the moment I use only one Geopackage-Layer for the polygons and would like to keep it that simple. But If there is a nice solution by using more layers or a "real" database, I'm okay with it! I'm running QGIS 3.10.


texts: QField - Display attributes table frame in pdf with editing offline mode

texts:
I'm using QField and QFieldSync to synchronize my data. In a project I build a composer with an attribute table frame, wich show attributes of my layer. When I run the project in QField and create a pdf, attribute table is empty.

I have noticed that attributes are well displayed in pdf table frame when QFieldSync configuration is set on "Copy" for my layer but not on "Editing offline".

QFieldSync rename layers with a "(offline)" suffix, and I think this could be the reason, so I tried to put an expression to deal with this but without success. I found options in dataDefinedProperties for this LayoutObject like dataDefinedAttributeTableSourceLayer but I don't know how to manipulate it.

How to show attributes table frame with QField and QFieldSync in editing offline mode ?

Is it posible to precise an expression in attribute table frame source layer properties (as specified in documentation) to prevent this ?


texts: Previous/next function for custom QGIS form

texts:
I'm currently working on a custom form in QGIS for archaeological application.

Basically, the custom form is used to read and write information about the walls and structures we found on field. Each item has two unique identification "numbers" : one is the "field number" we give to each find (it combines letters and numbers : zxxfxxx) and the other is a number automatically generated by the Postgresql database where all the data is stored.

I would like to implement previous and next buttons so that I can browse through my items. I have a rough idea of what I must do (i.e. attach a python script to my form and inside it work with iterations) but I don't know how to put all that together.

Could someone point me the way to go ? I had a look at the Selenext plugin but I don't know what is relevant in it in my case and how it could be attached to a form.


texts: How to use ArcPy to Check if Search Cursor Returns Anything

texts:
How can I use arcpy to check if a search cursor returns anything and then do something based on the answer.  So far I have this and it recognizes the positive - i.e. when i pass the query for 001 it gets the value, but when I substitute a bad value, the script does not print "No Search Cursor"

import arcpy

table = "A_Table"
query = "\"FeatureID\" = '001'"
theseRows = arcpy.SearchCursor(table, query)

for row in theseRows:
    if theseRows:
        print(row.getValue("FeatureID"))
    else:
        print "No Search Cursor"


texts:
I'm trying to use 'd.out.file' to save the displayed map to a PNG file. It failed and produced this message: "Command 'd.out.file' not yet implemented in the GUI. Try adding it as a command layer instead."

How can I resolve this? It also seems that I cannot run other d.* commands


texts: How to overwrite metadata in a tif file with GDAL?

texts:
I have a set of tif files generated from a process with multiple steps. Some metadata tags are added at different steps. I would like to update the metadata of those file before distributing the files. Using gdal_edit.py -mo, I can add some new lines in the metadata, but I would like to remove some the old lines (internal tags that we use for quality control).

Is there a way to clean the metadata and/or update some of the tags ?


texts: Polarity of Slope raster for use in Multi Criteria Evalution

texts:
I'm looking to use a slope raster in a MCE however i am struggling to ensure it has the correct polarity. The polarity of my other two layers is high values = good and low values = bad. The slope ranges from 0 (good) - 36 degrees (bad) and therefore currently the polarity is the opposite of what i require. Ive tried to multiply the grid by - 1 in raster calculator which gives me this.

When i attempt to standardise the grid so the values range from 1 (good) to 0 (bad) i get this.

Obviously i require all my layers to range from 0 - 1, and would like the areas of least slope to have high values once ive standardised them, does anyone know where i've gone wrong?


texts:
First of all I should say I am new to EPSG and QGIS world.Then
I want to Georeferced Raster image to EPSG:5235 in QGIS. But in Qgis I can't found above EPSG. How I add this EPSG to Qgis ?

Pls help me,

Thanks


texts: Finding average number of locations in each neighboring census tracts within a given distance

texts:
I have used the buffer and intersect tools to identify point locations within a certain distance of a census tract's boundaries.

I now want to calculate the average number of point locations in each of the neighboring census tract within a specified distance relative to a census tract of interest. I have included an image below for added clarity. In the census tract pictured, for example, I would like 7 different numbers that tell the sums of point locations in the seven neighboring census tracts within a specified distance.

I think that it might be possible to use the Spatial Join extract these averages in neighboring census tracts, but I am unsure of how exactly to specify it. In using the Spatial Join tool so far, I have only been able to calculate a sum of points within the entirety of each tract, which is not exactly my goal.

How would I specify this with the Spatial Join tool?

The layers used in my previous Spatial Join were a layer containing points within a distance of census tract boundaries and a layer of all census tracts.


texts: da.InsertCursor not adding correct point

texts:
I'm iterating a dictionary and inserting rows however, the shape is not inserted correct. All the rows ends up with the same shape.

My code:

    fc= r'C:\Connections\c\Osa\OSA@cB.sde\dvdb.Dv.E\Dv.Dv.Info'
    edit = arcpy.da.Editor(r'C:\Connections\c\Osa\OSA@cB.sde')
    IdShapeDict = dict()
    IdShapeDict['54255'] = (562122.9734, 6322034.7203)
    IdShapeDict['52367'] = (529588.5947, 6304689.1642)
    IdShapeDict['40965'] = (5563882.452, 6243252.3302)

    edit.startEditing(False, True)
    edit.startOperation()

    cursor = arcpy.da.InsertCursor(fc, ["Ke","SHAPE@XY"]):
        for strId, shape in IdShapeDict.iteritems():
            cursor.insertRow([strId,shape])

    del cursor
    edit.stopOperation()
    edit.stopEditing(True)

If a change the code to this:

    fc= r'C:\Connections\c\Osa\OSA@cB.sde\dvdb.Dv.E\Dv.Dv.Info'
    edit = arcpy.da.Editor(r'C:\Connections\c\Osa\OSA@cB.sde')
    IdShapeDict = dict()
    IdShapeDict['54255'] = (562122.9734, 6322034.7203)
    IdShapeDict['52367'] = (529588.5947, 6304689.1642)
    IdShapeDict['40965'] = (5563882.452, 6243252.3302)

    edit.startEditing(False, True)
    edit.startOperation()

    for strId, shape in IdShapeDict.iteritems():
        with cursor = arcpy.da.InsertCursor(fc, ["Ke","SHAPE@XY"]):
            cursor.insertRow([strId,shape])

    edit.stopOperation()
    edit.stopEditing(True)

It works. But keep creating and deleting cursors must have an impact on performance

To simplify the issue
This works as expected:

IdShapeDict = dict()
IdShapeDict['54255'] = (562122.9734, 6322034.7203)
IdShapeDict['52367'] = (529588.5947, 6304689.1642)
IdShapeDict['40965'] = (5563882.452, 6243252.3302)

edit = arcpy.da.Editor(r'C:\Connections\c\Osa\OSA@cB.sde')
edit.startEditing(False, True)
edit.startOperation()

for strId, shape in IdShapeDict.iteritems():
    with cursor = arcpy.da.InsertCursor(fc, ["Ke","SHAPE@XY"]):
        cursor.insertRow([strId,shape])

edit.stopOperation()
edit.stopEditing(True)

But it creates a new cursor for each iteration.

When moving the creation of the cursor out side the for loop all rows ends up with the same shape

IdShapeDict = dict()
IdShapeDict['54255'] = (562122.9734, 6322034.7203)
IdShapeDict['52367'] = (529588.5947, 6304689.1642)
IdShapeDict['40965'] = (5563882.452, 6243252.3302)

edit = arcpy.da.Editor(r'C:\Connections\c\Osa\OSA@cB.sde')
edit.startEditing(False, True)
edit.startOperation()

with cursor = arcpy.da.InsertCursor(fc, ["Ke","SHAPE@XY"]):
    for strId, shape in IdShapeDict.iteritems():
        cursor.insertRow([strId,shape])

 edit.stopOperation()
 edit.stopEditing(True)


texts: Display legend of WMS in layers panel in full extent

texts: When loading a WMS layer from e.g. http://www.wms.nrw.de/geobasis/wms_nw_dgk5 the legend of the layer with larger desciptions is displayed in the layer panel not to its full extent (see picture). When scrolling down, the layer collapse and so it is not possible to read the lower parts. Is there a possiblity to change?


texts: Batch task execution in google earth engine

texts: I need to make tens of clipped portions of an image, so I uploaded the boundaries in my assets and then loaded them all in code editor and did the clips with a for loop and exported them to the cloud. It will add a list of all the tasks in my task pane, but still I should go and run each of them individually. Is there any way to run all the tasks at once?


texts: Mapping habitat suitability using ArcGIS for Desktop?

texts:
I have been having somewhat of an issue recently with one of my projections. I am currently trying to map out the preferable habitat for a species in Arcgis.

To do this, I performed a GLMM in R and took the estimates (that go along with the chosen variables) and plugged them into the raster calculator as a means to get a projection of what it looks like.

The raster calculator formula is as follows:

-2.521 - (0.00006988 * "eucstream")
- (0.7524 * "Reclass_musi2")
+ (2.065 * "Reclass_nort5")
+ (1.325* "Reclass_stmt3")
- (0.1956 * "Reclass_summ2")
+ (0.02602 * "extrast")
- (1.180 * "Reclass_salt1")
+ (0.0002888 * "Extract_tif3")
+ (0.06971 * "northraster")
+ (0.01771 * "noreastrast")
+ (0.08177 * "northwestrast")
+ (0.4179 * "southrast")
+ (0.1964 * "southeastrast")
+ (0.5842 * "southwrast")
+ (0.3635 * "westrast")
+ (0.0001164 * "eucagri")
+ (0.00005264 * "eucurban")

As you can see, the red bits are the preferred (and my predicted) habitats and the green values represent areas that they seldom go (primarily agricultural areas).

http://imgur.com/dkuCjDA

The values for this projection are as follows: high 3.823017134 and low -3.139751899

Afterwards, I tried to apply a linear stretch to scale the predicted values between 0 - 1. I did this through this formula: 3.823017134 - (-3.139751899)

Thus, in the raster calculator: it looks like this: ("rastercalc112"--3.139751899) / 6.962769033 (with 6.962769033 being a result of the formula).

The pixel values have now changed to this: low = 3.35482e-011 and high = 1

The main issue that I am having now is that (and this applies to the original projection as well as the stretched) my pixel values are incredibly concentrated to one side of the histogram.

This is an issue because the preferable habitats appear in the middle of the histogram which somewhat skews up my analysis (as I am trying to eventually bin the data and run a Spearmans correlation).

One thing that I tried to test out involved going into properties and visually assessing exactly what pixels are causing this (through min-max):

here is 0.5-1.0 from the stretched data

 http://imgur.com/AEpxgzp

As you can see..from a cell range from 0.5 to 1.0, you pretty have all of the areas that would be considered 'perfect' (i.e. the red bits from the first projection, although, they are yellow in the picture).

That quickly changes when you switch the min-max to 0.75 to 1.0, where you only have a slight outline of the preferred areas. Thus, screwing up my data collection as most of the higher values are between 0.4-0.6.

I am new to GIS.

I am trying to simply distribute the raster pixels evenly so that the raster pixels that are causing an issue (anything above 6.0 which have very little to no data associated with them by the looks of it) are either removed, so that preferred habitat corresponds with the highest pixel values whereas the least preferred habitat type corresponds with the lowest (and of course, everything in between).

Can I evenly distribute these raster pixels so that areas of higher habitat suitability (which are now only located between 0.4-0.6, are redirected towards upper pixel values such as 1.0 rather than being literally in the middle of the histogram?


texts: I'm creating a web map in ArcGIS Online that references a MapServer from our internal server, which also has feature access enabled (FeatureServer) for editing.  When I add both services to Web AppBuilder I get duplicate layers in the Layer List (one for each layer individually, plus the "group" layer that normally shows for the MapServer alone).  Is there a way to display only the MapServer in the map and legend, but still allow feature editing through the FeatureServer?


texts: How to get, then export as .csv the lat lon centroids from shapefile layer

texts:
I have a shapefile with a whole lot of regions on it. I need to find the centroids for each of these regions in the shapefile (as latlons) and export this as a .csv

I have done extensive research, but no luck yet.

So far: I select the layer. I go to "Vector/Geometry Tools/Polygon centroids". hurrah, the centroids show up. BUT they are not yet added to my layer as lat lons. (although the info must be there because they are, you know, located.)

To get the lat lons added to my layer I go "Vector/Geometry Tools/ExportAdd Geometry Columns" BUT THEN when I look in my attribute table for this layer, some of the column headings have been changed (weird!)  :((  The good news is that there ARE the lat lons added.

To be more descriptive:

The first 16 column headings remain unchanged. That's good.
The 17th gets truncated from SA1_testData_SA1_AREA_SQM to SA1_testDa
The 18th column gets changed from SA1_testData_SA1_Dwellings to SA1_test_1
The 19th, 20th etc columns keep incrementing by one: SA1_test_2, SA1_test_3, SA1_test_4, etc
And then right at the end, the last 2 columns are XCOORD and YCOORD with the lat lons in them.

Yes, I have tried both "adding to this shapefile" AND also "saving as a new shapefile"... same outcome.  Why do the column headings get scrambled?


texts: Simultaneously habitat fragmentation metrics calculation within more buffers

texts:
Within several raster buffers I need to calculate habitat fragmentation metrics (i.e., number of patches, mean patch area, Euclidean distance among patches, landscape proportion) using the landscape statistics tool. Here an example:

However, because I have several buffers, I would like to know if there is an existing method in QGIS to calculate such metrics simultaneously and not considering each buffer one-to-one.


texts: Calculate area parts of polygon intersections

texts:
I have a problem of calculating area sizes of intersected polygons in postgis.
I would represent the problem in an example (please see the picture).
I have a gridsquare table with some attributes like color (green, blue, yellow)
I also have an affected_area table with some attributes (like area_reference, the affected area is not multipolygon, can contains areas with the same area_reference). I would like to create a query as an end solution which calculate the sum sizes of the green, blue, yellow part-grids where they are intersect with the affected_area. The sum of green areas highlighted in red in the picture example. One area can also be fully within a grid. Any idea where should I start? Probably I should go with something similar: http://www.bostongis.com/postgis_intersection_intersects.snippet


texts: Export entire Mosaic with "Equal-Arc" spatial reference

texts:
I have several mosaic datasets contain a total of about 1000 images, each about 500mb in size. Unfortunately our customer needs a file format different than the one we are using and needs the imagery all with the equal-arc spatial referencing.

Is there a way to export the imagery from the mosaic dataset so that the spatial referencing for each image is the same as the equal arc zone that corresponds to its location?

So far all I could think of was just choosing the Data frame as the spatial reference source when exporting the imagery as a different file type, but in that case I have to choose one of the arc zones for the data frame, meaning most of the data will be wrong (The mosaics cover entire countries).

Has anybody found any elegant solution to this problem.


texts: Get type of projection from pyproj CRS object

texts:
I'm used to using PROJ strings/dictionaries in python where I have access to a proj_dict['proj'] element. I have some calculations that I can only do with the PROJ.4 geos projection and otherwise fallback to a naive calculation. Since PROJ strings are not the preferred way of accessing this information anymore I'm wondering how I could do this same check with a pyproj CRS object.

Here's my example:

from pyproj import CRS
proj_dict = {'a': '6378137.0', 'b': '6356752.31414', 'h': '35786023.0', 'lon_0': '-89.5', 'proj': 'geos', 'sweep': 'x', 'units': 'm'}
c = CRS.from_dict(proj_dict)

The only thing I see related to "geos" is the coordinate operation method:

<Projected CRS: +proj=geos +a=6378137.0 +b=6356752.31414 +h=357860 ...>
Name: unknown
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- undefined
Coordinate Operation:
- name: unknown
- method: Geostationary Satellite (Sweep X)
Datum: unknown
- Ellipsoid: GRS 1980
- Prime Meridian: Greenwich

Is that the same thing as PROJ's proj parameter?


texts: Creating an interactive chloropleth map with Leaflet and Excel data

texts:
I am trying to create a chloropleth map with Leaflet.

I have an Excel spreadsheet with different columns such as "Project Name", "Project Description", "State" and "Province".

Is there a way to load the data from the Excel spreadsheet into a chloropleth map and have the details for each project display in a popup over each state/province?

I have tried using Leaflet before but I couldn't get the functionality out of it. I have also used ArcGIS which worked pretty well, but I cannot use it due to money constraints.


texts:
I am having trouble converting an ArcGIS Online JSON output to a more standard GeoJSON format. I have tried the solution in this thread, but the answer is either deprecated or no longer working. Every other solution I've come across seems to require an ArcGIS license - there must be an open source tool for this conversion, right?
And yes, I've tried modifying the 'f' query parameter to 'geojson' to no avail. The data is from 2014, so I guess this functionality is not supported.


texts:
I downloaded some road network data derived from OSM by GeoFabrik, see the initial state:

Near Székesfehérvár (Hungary), parallel lines are probably trunk or motorway roads.

I would like to simplify it like this (drawn in Paint):

I would like to use open source stuff for that (QGIS, GRASS), so I tried v.clean and v.simplify, v.generalize in some combination. Unfortunately, I cannot reach the desired result.

With v.clean 'snap' method, tolerance=50 I reached this:

But I don't need duplicated geometries, or identical zig-zag lines. Probably the v.generalize.network (Closeness, Betweeness threshold) would solve the problem, if it worked.

You can download the clipped geodata here:
https://github.com/ewirth/gistemp


texts:
I added a click event to the map.

But when I click on a layer on the map, I want to disable this click event.

I have to write a condition in the map click event.

How can I find if there is a layer when clicking on the map?

this._map.on('click', (e)=>{
 if(/* what should I write? */){
  // code
 }else{
  // code
 }
})


texts: PostGIS global Voronoi Thiessen Diagram Polygon

texts:
I am trying to create a global voronoi diagram for three points (cities): Miami, Boston and Honolulu. The geom column contains the points created by x and y.:

| index |   name   |      x       |      y      | geom |
|-------|----------|--------------|-------------|------|
|     0 | Miami    | -80.25551469 | 25.79555801 | .... |
|     1 | Boston   | -71.07082218 |  42.3585578 | .... |
|     2 | Honolulu | -157.8555857 | 21.30573449 | .... |

To create a Voronoi diagram I tried:

CREATE TABLE voronoiv07 AS
SELECT name,
    ST_VoronoiPolygons(geom,0,ST_Polygon(ST_GeomFromText('LINESTRING(-180 -90, 180 -90, 180  90, -180 90, -180-90 )'), 4326)) As geom
FROM cities

The result is a table with invalid geometries. Any ideas how to easily create a global voronoi diagram (with polygons that cross the -180 degree meridian)?


texts: Virtual layer with st_buffer from geopackage layer

texts:
I'm trying to make a virtual layer that has a buffer over a layer of lines.
If I run the following query on a layer in shapefile it works correctly

create view river as select st_buffer(geometry,50) from BCN0302L_RIO_line where etiqueta = 'Río Sarela';

but if I run the same query on the same data but in a layer in geopackage or spatialite, the query gives no error but does not show the geometry in the canvas

What's wrong?


texts:
I'm using GeoServer and want to add a WPS with default values.

How can I define an optional parameter as a WPS input?
I think, that I have to define this parameter with a multiplicity of 0, but how can I use a default value?

Do I have to add something in the ApplicationContext.xml?

Thanks for helping.


texts: Group location points that have a cumulative value 2000

texts:
I have a set of 400 locations (points) joined to a table having population attributes.
I want to form groups of nearby locations such that the cumulative population is 2000, group can have maximum 15 locations.
Is there a way to do it using QGIS and/or R?


texts:
Basically is it possible to copy georeferenced information from one raster (original) to another one (Photoshoped), not the tfw but actually header?

I looked hard but can't find any precise answer.


texts:
Below given is screen shot of GeoServer 2.13.0 installed on my machine.
Configuration of my machine is :
Intel(R) Core(TM) i7-3770 CPU @3.40 GHZ

Installed Memory 16.0 GB

OS Windows 10 64-bit

I am displaying WMS layer using leaflet to draw tiles on map.Database is SQL SERVER 2016. The problem is that mostly rendering of tiles is slow on map. When I investigated further, I found that Native JAI and Native JAI ImageIO is false and this was required for good performance : Geoserver Running in production environment

When I tried "Installing native JAI on Windows" in the following link :
Installing native JAI on Windows few links were not working. When I download JAI and JAI ImageIO from the given link they were throwing error that JRE 1.4 must be installed. However , on my machine JRE 1.8 is installed. Check the screenshot given below.

Any suggestion?


texts: Garmin GPX file, what units are the time values in?

texts:
I have a GPX file from a Garmin GPSMap 62s. The GPX file has a time column but the values are in the following format:

42079.98593

This point was collected at around 14:20 on 7/8/2015.

I need to convert the time to Pacific Standard Time so I can join the point locations with some other data.  If I convert the GPX file to feature class using ArcGIS 10.3 the time gets converted to 2015-03-15T 23:39:44Z.  Is that Zulu time on the wrong date?


texts:
I want to customize the legend box elements .Please check the screenshot

(I want to remove the highlighted information... it's the map service name and symbol in ArcMap)

How its possible? Is there any way or are there any parameters in the legend constructor??

I tried to modify following code but it was not successful...

Code:
dojo.connect(map,'onLayersAddResult',function(results){
          var layerInfo = dojo.map(results, function(layer,index){
            return {***layer: abc,title: 'new legend'***};
          });
          if(layerInfo.length > 0){
            var legendDijit = new esri.dijit.Legend({
              map:map,
              layerInfos:layerInfo
            },"legendDiv");
            legendDijit.startup();
          }
        });


texts: How to add a geometry to an existing feature that was created with NULL geometry?

texts:
I have a simple SpatiaLite geodb loaded in QGIS with one table. There is a POINT geometry field and some alphanumeric data. In some cases I can record data about a specific site, but the coordinates are not available at this time and will be retrieved and recorded later. In some other cases I have the coordinates and I want to insert them right into the geometry field, preferably by clicking on the map canvas in the same way available for adding new features (having a separate text field for the coordinates and generating the geometry from that is a no-go).

How do I add a geometry to an existing feature that was created with a NULL geometry field?


texts: Multiple Input Feature Classes with Different Names - Model Builder

texts:
I have created the model below to create "straight lines" from a route (Feature Class is the Route in the model). This model will save me a lot of time compared to doing it manual. The manual steps are as follows: File> Add Data> Add Route Data... New window appears, the route and route identifier are already selected from my TOC, but all the other information must be entered then you hit ok.. route event is created, then the event must be exported to a gdb, added to the map, then symbolized. My issue here is, I would like for the model to recognize my feature class (the route) from within my TOC like the manual window does. I do not want to run my model and every time have to keep adding the route.

I have about 30 maps this model will be run off of. Each route has a different name. Would it be possible to give them an alias name such as "Route" and set the model's "Feature Class" to recognize that instead? So as long as in each different map the route is named "Route" the model will recognize that and I will not have to manually add the feature class (route). Not opposed to python.


texts: Layer spatial reference system units EPSG::32633

texts:
I polygonized a raster with the following code:

from osgeo import gdal, ogr, osr
gdal.UseExceptions()

srs = osr.SpatialReference()
srs.ImportFromEPSG(32633)

ds = gdal.Open('H:/Threshold/savi_nsdvi66.tif')
srcband = ds.GetRasterBand(1)
dst_layername = 'H:/Threshold/c11111111'
drv = ogr.GetDriverByName('ESRI Shapefile')
dst_ds = drv.CreateDataSource(dst_layername + '.shp')

dst_layer = dst_ds.CreateLayer(dst_layername, srs = srs,   geom_type=ogr.wkbMultiPolygon )
polygone = gdal.Polygonize( srcband, None, dst_layer, -1, [], callback=None )

Everything works fine just one point is freaking me out. The shapefile shows the correct reference system (+proj=utm +zone=33 +datum=WGS84 +units=m +no_defs) in QGIS but I receive wrong layer spatial system units (xmin,xmax,ymin,ymax) dimensions from the input raster file not the usual EPSG coordinates.


texts:
I have points like this in a leaflet map.

The points come originally from a GeoJSON file and I used L.circleMarker to convert them to circles. How can I create a legend for these points? Based on the Interactive Choropleth Map Example  I was able to generate a legend like this:

The code I used is this:

var legend = L.control({position: 'bottomleft'});

legend.onAdd = function (map) {

  var div = L.DomUtil.create('div', 'info legend');

    categories = ['STX','HHX','HF','STX/HF'];

    for (var i = 0; i < categories.length; i++) {
        div.innerHTML +=
            ' ' +
             (categories[i] ? categories[i] + '
' : '+');
    }

    return div;

But I like to also match the form (circles) to the legend, and not only the color.


texts: My company currently uses MapInfo Professional to map our client's addresses into territories.  I have been asked to figure out a mechanism to do the same thing but real-time via a web service call.  Is there any way to store the data from teh Mapnfo TAB files in a database and query it with lat/long in order to determine which "object" the lat/long belongs in?  I've seen that there are multiple data formats for geospatial information, so if it isn't possible with the MapInfo TAB files, would it be possible with other formats?  Am I coming at this issue in the correct way or is there alternative mechanisms for solving the problem?


texts: I am running topology check "layer should have no gaps" and the result shows one error with feature id 0. Double clicking on the error shows a large Dodecagon that is not part of my layer. How can I further troubleshoot and identify the gap?


texts: FME PostGIS writer define column as serial/PrimaryKey with pre-existing id values

texts:
With FME I read a table that has an id column with pre-existing values, like

id       | some_attributes...
--------------------------------
56784645 | ...
56478687 | ...
...

I would like to write these into a PostgreSQL/PostGIS database and for further data input I define the id column as serial / PrimaryKey. Because of joins to other tables, the id values need to be preserved.

In the resulting table my id values are replaced by the serials:

id       | some_attributes...
--------------------------------
1        | ...
2        | ...
...

QUESTION: How can I set up FME PostGIS writer in a way that defines a column as serial, but preserves my pre-existing values for this column?


texts: Error when importing shapefile as vector layer into GRASS GIS

texts:
I'm tried to importing an Esri shapefile into GRASS GIS using  both v.in.ogr and v.in.import however methods fail with the following error message:

WARNING: Table <roadsSriLanka> linked to vector map <roadsSriLanka> does
         not exist
Creating attribute table for layer <edges>...
DBMI-SQLite driver error:
Error in sqlite3_prepare():
near "from": syntax error

DBMI-SQLite driver error:
Error in sqlite3_prepare():
near "from": syntax error

ERROR: Unable to create table: 'create table roadsSriLanka (cat integer,
       access varchar ( 80 ), bridge varchar ( 80 ), from varchar ( 80 ),
       highway varchar ( 80 ), junction varchar ( 80 ), key varchar ( 80 ),
       landuse varchar ( 80 ), lanes varchar ( 80 ), length varchar ( 80 ),
       maxspeed varchar ( 80 ), name varchar ( 161 ), oneway varchar ( 80
       ), osmid varchar ( 231 ), ref varchar ( 80 ), service varchar ( 80
       ), to varchar ( 80 ), tunnel varchar ( 80 ), width varchar ( 80 ))'

This shapefile was downloaded using the OSMnx Python package. I used another file which I downloaded using the same which imported without an issue


texts: Histogram Matching using Python to improve mosaicking process of multiple overlapping rasters?

texts:
I am trying to do histogram matching using Python to improve the mosaicking process of multiple overlapping rasters.  I am basing my code on that found at:

http://www.idlcoyote.com/ip_tips/histomatch.html

To date, I have managed to clip the overlapping area of two adjacent rasters and flatten the array.

so I have two 1 dimensional arrays of the same length.

I have then written the following code based on that found at the above website. In the code shown I have substituted two very small datasets for the gd and bd images.

import matplotlib.pyplot as plt
from scipy.interpolate import interp1d

bins = range(0,100, 10)

gd_hist = [1,2,3,4,5,4,3,2,1]

bd_hist = [2,4,6,8,10,8,6,4,2]

nPixels = len(gd_hist)

# here we are creating the cumulative distribution frequency for the bad image
cdf_bd = []
for k in range(0, len(bins)-1):
    b = sum(bd_hist[:k])
    cdf_bd.append(float(b)/nPixels)

# here we are creating the cumulative distribution frequency for the good image
cdf_gd = []
for l in range(0, len(bins)-1):
    g = sum(gd_hist[:l])
    cdf_gd.append(float(g)/nPixels)

# we plot a histogram of the number of
plt.plot(bins[1:], gd_hist, 'g')
plt.plot(bins[1:], bd_hist, 'r--')
plt.show()

# we plot the cumulative distribution frequencies of both images
plt.plot(bins[1:], cdf_gd, 'g')
plt.plot(bins[1:], cdf_bd, 'r--')
plt.show()

z = []
# loop through the bins
for m in range(0, len(bins)-1):

    p = [cdf_bd.index(b) for b in cdf_bd if b < cdf_gd[m]]
    if len(p) == 0:
        z.append(0)
    else:
        # if p is not empty, find the last value in the list p
        lastval = p[len(p)-1]

        # find the bin value at index 'lastval'
        z.append(bins[lastval])

plt.plot(bins[1:], z, 'g')
plt.show()

# look into the 'bounds_error'
fi = interp1d(bins[1:], z, bounds_error=False, kind='cubic')
plt.plot(bins[1:], gd_hist, 'g')
plt.show
plt.plot(bins[1:], fi(bd_hist), 'r--')
plt.show()

My program plots the histograms and cumulative frequency distributions successfully...and I thought that I had the part of getting the transformation function 'z' correct....but then when I use the distribution function 'fi' on the 'bd_hist' to try to match it to the gd dataset it all goes pear-shaped.

I am not a mathematician and it is highly likely I have overlooked something fairly obvious.


texts: Correct WFS URL Format (Gaia WFS-T Extender)

texts: I am trying to add a GeoServer WFS layer to Gaia WFS-T Extender for editing. The available layers are available alright but I am unable to load them into Gaia. Using the same URL I was able to add and edit the layer in both Quantum GIS and UDig. Any hints on where I could be missing it with Gaia? [ why Gaia: It is light-weight and packs functionality I'm looking for]


texts: Adding Earth Engine Raster to Mapbox map

texts:
I have a bunch of data created in Google Earth Engine and I want to add some raster layers to a map created with Mapbox GL JS.

One way I managed to do it is via the Mapbox Studio tileset upload. However, this seems quite complex, download image from EE, process to 8-bit, save as rgb, upload,...

Another answer is Adding layer made in Google Earth Engine Code Editor to Mapbox, which suggests to upload it to cloud storage bucket and add it by mapbox-gl using addSource. However, until now I didn't manage to get it running this way.

This example here is pretty much what I want to do, unfortunately I don't understand how to use it for my own data.


texts: How to produce a simple ZTV/ZVI in openWind?

texts:
I'm struggling to get to grips with openWind. I need to produce a ZTV showing a 5km radius around the chimney of a proposed industrial development. I'm not sure of the height of the stack yet but for argument's sake let's say it's 60m. I heard that I could pretend that this is a turbine by creating a new turbine type. I set the hub at 60m but it asks for the rotor diameter so not sure what to do there.

I'm using Ordnance Survey (UK) open data - have downloaded the raster and DTM for the area I need and added them as new layers in openWind.

Here's where I'm starting to get confused and just want to make sure I'm on the right track. I've taken the DTM and set it to terrain elevation. Then created a raster based on this which looks right to me. I've set my custom turbine position and now set a ZVI going. It's taking a very long time so I'll have to wait and see if it's right. Is there a way to set just the 5km radius though for the ZVI?

Any help greatly received.


texts: Possible to run "Merge shapefiles to one" algorithm from Python script?

texts:
Developing a python script in QGIS 2.2, I can run algorithms from the Processing plugin fine such as those from QGIS geoalgorithms and SAGA.

Can you do the same for the algorithm "Merge shapefiles to one" from the Vector menu and how can you call it from the Python console?

The "Merge shapefiles to one" algorithm is from the Ftools plugin.

Found a link which might be useful: Accessing QGIS plugins in Python?


texts: Why is ArcGIS Tabulate Area Better to Calculate Unique Landuse Values in Zone?

texts:
I am not a formally trained GIS professional nor is GIS a primary tool in my job. And therefore, without an incentive to optimize GIS workflows, I  may not understand why one tool, although seeming to accomplish the same goal as another, may be better for a given task .

Sometimes answers or comments are given that imply or clearly state that the best tool to accomplish a given task is tool x and perhaps not the one suggested in the accepted/upvoted answer. Missing is a reason the other tool is not as good or the other better. I think this is important, to me anyway, to educate the asker and the answerer who has provided a not so good answer.

In this case, I would like to understand why Zonal Statistics is not the best tool for Count raster values in each polygon?


texts: Alternative North Arrow using R GIStools package?

texts:
I'm really digging the R package GISTools for making some basic maps. The default north arrow, created with north.arrow looks ok, but I'm wondering if it would be possible to use a different north arrow. It seems like no other arrow is available in the package, so how could I create my own north arrow? The documentation for the package states

  Draws a north arrow on a map. The arrow itself is drawn using polygon and any extra parameters are passed to this call.

So I would need to create the arrow as a polygon? How could I implement this? And would it be possible to modify the source code of north.arrow to offer several different arrows?


texts: Where is Area and Volume statistics tool in ArcMap 10.1?

texts:
I want to measure surface volume of an area with different Plane Height above reference plane. There should be a tool for statistic calculation but I cannot find it in 3D analyst toolbox.
Is it exist in GIS 10.1 ?
Or which tool can help me to measure ?


texts:
I am running a Geoserver instance in Tomcat and want to know if it is possible to use Windows Authentication (Keberos) to authenticate users in Geoserver.

My idea is to setup Geoserver to use an Active Directory to find users and then let the client send a token like "Negotiate dxs24sdsedsx!dwwe" (Kerberos) to Geoserver.

Is it possible for Geoserver to authenticate that kind of header?


texts: How to display the 'Stop Editing' window of a layer using PyQGIS?

texts:
I'm trying to implement a function that would warn the user when he changes active layer and has just made changes on it that his changes have not been saved.
I found this forum: Prevent QGIS user to change layers without saving? Which inspired me. I tried to translate it to QGIS 3 (see the code below).
def check_edit():
    layers = iface.editableLayers()
    if layers:
        layers = [layer for layer in QgsProject.instance().mapLayers().values()]
        iface.setActiveLayer(layers[0])
        iface.activeLayer()
        iface.actionToggleEditing().trigger()

iface.currentLayerChanged.connect(check_edit)

But the result of it is not as expected because the changes are not saved. The changes appear, but the layer is not saved.
So I kept a little different code to which I added a pushMessage which warns the user which layer is concerned.
def check_edit():
    layers = iface.editableLayers()
    if layers:
        iface.messageBar().pushMessage("Attention : ", "Couche "+ iface.activeLayer().name() + " en édition et non sauvegardée !", level=Qgis.Warning, duration=7)
        iface.activeLayer()
        iface.actionToggleEditing().trigger()

iface.currentLayerChanged.connect(check_edit)

I therefore do not find my code very efficient in the sense that the user is warned but these modifications are not saved.
So I would like to know if it is possible to display this pop-up (screenshot) which would leave the choice to the user to save or not his modifications?


texts: Creating Ellipse in WGS84 with metric parameters

texts:
I have points in WGS84 (x,y), radiuses in meters (rx, ry), rotation angle in degrees. I would like to create a new WGS84 geometry.

I've tried this but the major/minor distances seem to be off by a few hundred meters in some cases.

Ellipse (x,y,rx,ry,rotation,#of segments in 1/4 of ellipse)

SELECT
transform(
 translate(
  rotate(
   scale(
    buffer(
     transform(setsrid(makepoint(0,0),4326),3395),
     1,
   $6),
  $3,$4),
 $5),
x(transform(setsrid(makepoint($1,$2),4326),3395)),
y(transform(setsrid(makepoint($1,$2),4326),3395))),4326)


texts:
I was able to "combine my shapefile and CSV data have made a lot of progress on my project.

Now I want to display the attribute data in a popup when I either click or hover over a country... I've done some goggling and looked at some examples but can't seem to make what I've seen work...

Can someone recommend a link or example I can look at?

I'm using MapServer and OpenLayers

Let me know if you'd like to see my script.


texts: TypeError when using gdal.Warp in Jupyter Notebook

texts:
I try to use the gdal python binding to project a raster in World Mollweide.
import gdal

map_raster_proj = pm.getTmpDir() + Path(map_raster).stem + '_proj.tif'
input_ = gdal.Open(map_raster)
gdal.Warp(input_, map_raster_proj, dstSRS='ESRI:54009')

but i get the following error :

TypeError: in method 'wrapper_GDALWarpDestDS', argument 4 of type 'GDALWarpAppOptions *'

If I use gdalwarp from a terminal it runs without bug. Am I doing something wrong?


texts: GeoServer 2.17 - added Shapefile with geometry of type "3D Polygon", rendering fails exception "Must have at least 2 spatial dimensions"

texts:
I've downloaded a public shapefile and added it to Geoserver 2.16 just fine. QGIS as well. I added it to a GeoServer 2.17 as well, but when trying to preview, it causes the following exception

2020-05-22 08:28:38,325 ERROR [renderer.lite] - Must have at least 2 spatial dimensions
java.lang.IllegalArgumentException: Must have at least 2 spatial dimensions
    at org.locationtech.jts.geom.impl.PackedCoordinateSequence.<init>(PackedCoordinateSequence.java:57)
    at org.geotools.geometry.jts.LiteCoordinateSequence.<init>(LiteCoordinateSequence.java:151)
    at org.geotools.geometry.jts.LiteCoordinateSequence.cloneGeometry(LiteCoordinateSequence.java:371)
    at org.geotools.geometry.jts.LiteCoordinateSequence.cloneGeometry(LiteCoordinateSequence.java:316)
    at org.geotools.geometry.jts.LiteCoordinateSequence.cloneGeometry(LiteCoordinateSequence.java:302)
    at org.geotools.geometry.jts.LiteCoordinateSequence.cloneGeometry(LiteCoordinateSequence.java:384)
    at org.geotools.geometry.jts.LiteCoordinateSequence.cloneGeometry(LiteCoordinateSequence.java:306)
    at org.geotools.geometry.jts.LiteCoordinateSequence.cloneGeometry(LiteCoordinateSequence.java:312)
    at org.geotools.geometry.jts.LiteShape2.<init>(LiteShape2.java:124)
    at org.geotools.geometry.jts.LiteShape2.<init>(LiteShape2.java:99)
    at org.geotools.renderer.lite.StreamingRenderer.processSymbolizers(StreamingRenderer.java:3000)
    at org.geotools.renderer.lite.StreamingRenderer.processFeature(StreamingRenderer.java:2807)
    at org.geotools.renderer.lite.StreamingRenderer.drawPlain(StreamingRenderer.java:2644)
    at org.geotools.renderer.lite.StreamingRenderer.processStylers(StreamingRenderer.java:2243)
    at org.geotools.renderer.lite.StreamingRenderer.paint(StreamingRenderer.java:900)
    at org.geoserver.wms.map.RenderedImageMapOutputFormat.produceMap(RenderedImageMapOutputFormat.java:601)
    at org.geoserver.wms.map.RenderedImageMapOutputFormat.produceMap(RenderedImageMapOutputFormat.java:275)
    at org.geoserver.wms.map.RenderedImageMapOutputFormat.produceMap(RenderedImageMapOutputFormat.java:135)
    at org.geoserver.wms.GetMap.executeInternal(GetMap.java:749)
    at org.geoserver.wms.GetMap.run(GetMap.java:300)

ogrinfo indicates the layer contains 3D polygons.  And I see each point does have three coordinates.  I don't know if the third value is in feet or meters.

To try this yourself, go here: https://openhartford-hartfordgis.opendata.arcgis.com/datasets/building and download the shapefile.  Add it to GeoServer 2.17, Preview it with Openlayers, and you'll get the exception.

I've worked around this by using QGIS and Exporting to a new Shapefile, but specifying that the geometries should be Polygons.  This 2 dimensional export renders fine in GeoServer 2.17.


texts: Method in QGIS to classify veg from HiRes geotiffs created in Pix4D?

texts:
I have 10 2D geoTIFF orthophotos (~2GB each, ~3km2 each, all 4.2cm/pixel) that I generated from Vis RGB 12MP images in Pix4D for a desert vegetation survey over 30km2 for a development. I need to calculate the veg cover under a polygon in each ortho. After a bit of searching I opted for the SCP plugin by Luca Congedo. It seemed to work a treat, I can understand the UI just enough to get it to work and can overlay my mask layer to only classify that part.

Trouble is for one polygon is has taken often days to process, often fails, and the resulting classified shapes have over 1M features, making area calculation laggy and frankly I haven't managed to get very far with it. I've even divided each mask into 4 to reduce the area to classify and it takes almost the same amount of time for a quarter of the area.

Eg. I completed one small area (probably 1/50 of the total area), took 22 hours, resulting .shp is 443MB with 1.5M features and is practically unworkable due to lag for any type of action. I removed the non-veg classes I'm not interested in and Save As'd and it has so far taken 1 hour and counting...

I'm using 2 identical Dell XPS 9550's: quad-core i7 6700HQ 2.6 GHz; 2133MHz 32GB RAM; 512GB SSD; NVIDIA GTX960M 2GB; so no room for change there. Although interestingly neither my CPU, RAM or GPU ever get close to maxing out, so I'm thinking the plugin just isn't optimized for large high res images.

Can anyone suggest a faster method using QGIS?  Or maybe there's some setting in SCP that I need to tweak?


texts:
Asked another question below and this will develop on from this.  The Bedrock 625K BGS mapping is all in UPPERCASE.  It looks clumsey and wrong in the legend.  I wish to change it to Sentence Case.  I know the long winded way is to do each entry manually...there is over 7000 of them...so no.

Working in either QGIS 1.8.4 or 2.0


texts:
I am experiencing this issue (Fig.1)
Basically, it seems that my data (red dots) are only partially falling on the world shapefile. Both the layers are projected in WGS 84 EPSG: 4326.
How can I solve the issue?
Shall I reproject the data or find a world shapefile that is showed in a different way?
Sorry for the non-technical language but I find difficult to explain the specific issue.


texts: How to use the 'to_dms' function in a CRS conversion expression in QGIS?

texts:
I need to create in a point layer attribute table new fields with the points coordinates converted in different geographic CRS through an expression like:

x(transform($geometry, 'EPSG:xxxxx', 'EPSG:xxxxx'))

Is it possible (and how) to integrate in this expression the function 'to_dms' in order to convert directly the decimal degrees output to DMS degrees?


texts: How to change standard button bitmap of esri tool/command in ArcEngine 10

texts: I'd like to use some of the standard 'out of the box' esri tools e.g. Measure Tool in an ArcEngine app, but I'd like to change the default bitmap of the tool button. Any suggestions would be gratefully received.


texts:
I'm working on OTP. I would like to create a graph for Munich (which requieres GTFS and OSM data).

However, there is no gtfs available, so do you know how to convert osm data to GTFS?


texts: Splitting line in proportional parts in QGIS

texts:
I have recorded flight paths of birds as lines. Now I want to transform the lines to points, which represent the bird's position in a minute interval. For this, I want to divide the lines by the time I have observed the flight.

So a flight track I could watch for 5 minutes I want to divide into 5 parts, a track I have watched for 3 minutes, I want to divide into 3 parts and so on. Afterwards, I take the middle point of these segments as the flight position in the respective minute.

My first approach was to calculate the length of the segments ($length/minutes) and store the information in a new field (=segment). Afterwards, I used the plugin 'Station lines' to set lines after each segment. Then I used the QGIS tool 'Divide lines with lines'.

First it seemed to work, but then I realized that the function also cut the line at every self- intersection. As the flight tracks contain many loops I ended up with far more segments as expected.

Does anybody know a simple and safe way to split lines into proportional parts?


texts: Creating one single coverage polygon feature class from XY point data?

texts:
I'm using ArcMap 10.2, and I need to map areas covered by London postcodes (E1 0xx, E1 1xx, E1 2xx, and so on). What I have is the coordinates for every single postcode (310,000 in total) in an Excel sheet, so I can create point feature classes, which I then want to group as explained above.

How can I create a polygon shapefile that will include a group of coordinates within E1 0xx postcode?

For example. I want that polygon shapefile to cover the are of those postcodes specifically. In the end I need to end up with multiple polygon shapefiles that will represent different Postcodes districts (E1 0xx, E1 1xx...).
Tried the Coverage tool, but it's not seem to be doing what I expected it to. Just installed XTools Pro, maybe that could help me.


texts: How to convert picture to annotation in ArcGIS 10.1?

texts:
I have been tasked to display analytical data tables (currently in MS Excel) in an ArcGIS (v 10.1) map document (.mxd) and would also like to convert them to map document annotation, if possible.  Ideally, I would like the embedded excel tables to stay the same size and in the same place when I am moving and zooming around the map document.  Each analytical table is linked to a specific monitoring well.

If I simply cut-and-paste the individual tables (e.g. below) into my map, they will appear as "pictures".  I know that pictures are neither graphics (point, line, polygon, text) nor a feature.  And when I right-click on the pasted table, I am not allowed to "Convert to Graphics" (as it is greyed out).  How can I solve this problem?  Is there a work-around?

Example table I want to add to my map document (which appears as a picture in ArcGIS):


texts: Creating a dummy variable for properties that intersect the buffer I created

texts: I created a 1000 foot buffer around waterways for a single county in the U.S. I need to create a dummy variable for the properties that intersect the buffer I created. It would be 1 if the property intersects and 0 if not. I need to add this dummy variable to the parcel shapefile attribute table.


texts: St_Azimuth : Difference between Google Earth and Postgis

texts:
I'm using Postgis to find angles between sites. End users will use Google Earth to play with these results.

According to st_azimuth, I have exactly 315.002 degrees between these two points :

SELECT degrees(ST_Azimuth( st_setsrid(ST_Point(48.718544, -3.554562),4326), st_setsrid(ST_Point(48.116875, -2.95284),4326) )) AS degA_B

But when I verify them on Google Earth Pro, the result is not just slightly different : the ruler give 326.65° (!)
AFAIK, Google Earth is WGS84 based too.

St_distance does'nt much agree with Google earh : the first  94.33333316013 km and the latter 80.37km

Frankly, I'm more confident on Postgis, but I don't know who to believe...

Have I made a mistake?


texts: Getting route segments from OSM data and assigning ETA for them

texts:
I have a problem of assigning possible ETA for a route.
So far, I have the following:

OSM data for a city.
I have open source data of taxi trails for the same city.

Now, I am planning to divide the OSM data into route segments with a unique identifier for each. Meaning each segment is a route without any other route intersection/branching off and a value will be assigned.

For each segment, I will look up in my taxi data to see if it has travelled that segment. If yes, how much time did it take.

I wanted to know if this has been done already. I looked at GraphHopper in which I can get routes, but not segments.

Before going into pgRouting, I wanted to ask if there is any off-the-shelf best possible open source solution for the same?


texts: How to find the center point which lies on the LineString geometry?

texts:
I have a LineString geometry which is U shaped like and I'd like to find the center point of the geometry which lies on the LineString. I've already used st_centroid but the returned geometry isn't on the LineString.

I would like to know which function should be used in PostGIS to achieve that?


texts: Python in ModelBuilder calculate field tool

texts:
Code works on ArcMap "Field Calculator..." when I click on attribute and calculate.

But if does not work in ModelBuilder calculate field tool.

Any ideas?

Pre-Logic Script Code:

  def ifBlock(Country):
    if Country == 'China':
     return "A"
    elif Country == 'Canada':
     return "B"
    elif Country == 'Russia':
     return "C"

result=

   ifBlock(!Country!)


texts: Extracting Text Objects from Vector Graphic Software?

texts:
Within a Canvas X 2017 GIS file, (an ACD Systems software) I have a number of vector graphics layers, all of which contain thousands of manually entered text objects only conceptually paired with line objects (within the same layers). I hear that Canvas is similar to illustrator, although I have had too little experience with illustrator to verify this. Using Canvas's GIS capabilities, I have given my vector data spatial coordinates and then exported the layers as shapefiles. (Each layer becomes its own shapefile). However, the text objects are not transported into Arcmap (ArcGIS 10.4) within the shapefile as best as I can tell.

(Image of the data in Canvas X)

I am trying to research the best way to export the text objects and then to open them up within geospatial software. My goal is to be able to spatially associate them with their corresponding vector objects and have the GIS program populate a vector object's attribute field with the text object's text. Other formats (excluding raster) I can export the data as are TXT, PDF, and KMZ, and KML files. GIS software that I have access to include ArcGIS, QGIS, Global Mapper, and MapInfo. I have included a link to a PDF and KML sample of the data at the bottom of this post.

I would rather do this than have to export a rendered geospatial raster image, open that up side-by-side with the same layer's shapefile, and manually type the numbers into the attribute table for thousands of objects.

I have found very little on this topic online or even in person at academic institutions, but hopefully someone out there has familiarity with my question?

https://drive.google.com/file/d/1YpDyCUZqvSgxxdrDvSacEsN4wfPw87ex/view?usp=sharing
https://drive.google.com/file/d/1xxZfzPEqLwuyiMQYaIwqGO2Fip5vOG_l/view?usp=sharing


texts: Unable to get Feature Layer to apply Unique Value Renderer in ArcGIS Javascript API

texts:
I am trying to make a unique value renderer on a map using the JavaScript API. I am using a layer that represents continents. I have used the samples and modified the code for my own layer, but I can not get the feature layer to populate or the renderer to work. I have changed 'Continent' to 'CONTINENT' and 'continent' in the script but it does not matter. The layer does not populate and therefor it does not generate the unique value renderer. The developer tools and inspect tool in Chrome does not throw me any errors.  I am using these resources:

http://services.arcgis.com/BG6nSlhZSAWtExvp/arcgis/rest/services/Continents_simple/FeatureServer
https://developers.arcgis.com/javascript/3/jssamples/renderer_unique_value.html

Here is my code:

<!DOCTYPE html>

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

<meta name="viewport" content="initial-scale=1, maximum-scale=1,user-
 scalable=no">
<title>Unique Value Renderer</title>
<link rel="stylesheet" href="https://js.arcgis.com/3.22/esri/css/esri.css">


texts:
I want to mask pixels based on the quality band information for image collection Landsat 5. this is a code but the image-collection is empty afterward. I would be thankful if anyone could help me.

var startDate = '1992-01-01'
var endDate =   '1993-01-01'
var cloudCover = 30
var collection1992 = ee.ImageCollection('LANDSAT/LT05/C01/T1_SR')
                        .filter(ee.Filter.eq('WRS_PATH', 170))
                        .filter(ee.Filter.eq('WRS_ROW', 32))
                        .filterDate(startDate,endDate)
                        .filterMetadata('CLOUD_COVER','less_than',cloudCover)
                        .sort('system:time_start')
// cloud mask 66,130
var collection1992_noclouds = collection1992.map(function(img) {
                var mask1 = img.select('pixel_qa').eq(66).eq(130)
                return img.UpdateMask(mask1)
                    })
Map.addLayer(collection1992_noclouds, {}, 'collection1992_noclouds')