texts:
I need help with WFS filtering for GML. I want to filter only by one property (name of administrative unit) - that is all.

Lets say I would like to get the second member of this GML:
http://services.cuzk.cz/wfs/inspire-au-wfs.asp?service=wfs&version=2.0.0&request=GetFeature&typeNames=au:AdministrativeUnit&count=2

The name of administrative unit is in this tag: <gn:text>Hlavní město Praha</gn:text> (path: au:name > gn:GeographicalName > gn:spelling > gn:SpellingOfName > gn:text)

Here are some of my tests which do not work:

http://services.cuzk.cz/wfs/inspire-au-wfs.asp?service=wfs&version=2.0.0&request=GetFeature&typeNames=au:AdministrativeUnit&FILTER=<Filter><PropertyIsEqualTo><ValueReference>gn:text</ValueReference><Literal>Hlavní město Praha</Literal></PropertyIsEqualTo></Filter>

When I use ogr2ogr to transform the GML to GeoJSON (the WFS returns only GML) then the name of the administrative unit is in "text" attribute.

http://services.cuzk.cz/wfs/inspire-au-wfs.asp?service=wfs&version=2.0.0&request=GetFeature&typeNames=au:AdministrativeUnit&FILTER=<Filter><PropertyIsEqualTo><PropertyName>text</PropertyName><Literal>Hlavní město Praha</Literal></PropertyIsEqualTo></Filter>

This one returns whole layer instead of the 1 feature/member:

http://services.cuzk.cz/wfs/inspire-au-wfs.asp?service=wfs&version=2.0.0&request=GetFeature&typeNames=au:AdministrativeUnit&CQL_FILTER=text='Hlavní město Praha'

Here is meta for the layer & service:

http://services.cuzk.cz/wfs/inspire-au-wfs.asp?service=wfs&version=2.0.0&request=DescribeFeatureType&typeNames=au:AdministrativeUnit

http://services.cuzk.cz/wfs/inspire-au-wfs.asp?service=wfs&version=2.0.0&request=GetCapabilities

Any idea how to filter the nested tags in GML?


texts: I have geography data stored in SQL Server and I am trying to merge polygons from two data sets together.  The problem is that they do not line up exactly, so unioning them results in narrow gaps along the boundaries between the two sets.  What's the easiest way to removed these seams?  I am currently doing this in SQL Server, so if I could do this with tools SQL Server provides, that would be ideal.  Otherwise I'm open to using another, preferably free tool.


texts: Write string to status bar that doesn't go away

texts:
I'm attempting to write the number of items matching a given query to the statusbar.  when i debug through the code, i am able to step through and see the message hit the statusbar, but then it immediately disappears.  Below is the function in question.

public class CollisionViewExtension : Extension
{
     private static IQueryable<Collision> AllCollisions
    {
        get { return _allCollisions ?? (_allCollisions = Database.Query<Collision>(Collision.Query).ToList().AsQueryable()); }
    }
    ...
    public static void DrawCollisionPoints(string whereClause)
    {
        var matches = new List<Collision>();
        try
        {
            matches = AllCollisions.Where(string.IsNullOrWhiteSpace(whereClause) ? "1=1" : whereClause, null).ToList();

            if (CurrentSpatialQueryShape != null)
            {
                matches = matches.Where(x => IsPointInPolygon(CurrentSpatialQueryShape, (IPoint)x.Point)).ToList();
            }
        }
        catch (Exception)
        {
            System.Windows.Forms.MessageBox.Show("Invalid query");
        }

        FeatureInserter.InsertFeatures(matches);

        var layer = new FeatureLayerLocator(ArcMap.Document.FocusMap).GetFeatureLayerByModelName("CollisionView-Collisions");
        if (layer == null)
        {
            ArcMap.Document.FocusMap.AddLayer(FeatureLayer);
        }

        ArcMap.Application.Caption = string.Format("{0} collisions", matches.Count);
        Application.DoEvents();
        ArcMap.Application.StatusBar.set_Message((int)esriStatusBarPanes.esriStatusMain, string.Format("{0} collisions", matches.Count));

        ArcMap.Document.ActiveView.PartialRefresh(esriViewDrawPhase.esriViewGeography, FeatureLayer, null);
        ArcMap.Application.StatusBar.set_Message((int)esriStatusBarPanes.esriStatusMain, string.Format("{0} collisions", matches.Count));

        CurrentQuery = whereClause;
    }
    ...
}

I feel like this is pretty much how I've seen it used in all the examples, but for some reason my message is not displayed.  Asked another, probably simpler way, is there any way to show the user a count of the features in a particular layer in the statusbar?


texts: Calculating Slope, Flow Accumulation and Topographic Wetness Index (TWI) using ArcGIS for Desktop?

texts:
I am working on a project that requires me to extract values from  rasters (viz the slope and TWI) for a set of point dataset that I have. I tried calculating the Slope, Flow accumulation and the TWI from the DEM (using ArcGIS 10.1) but somehow the slope and the TWI raster outputs seems erroneous!

For  slope, I get values ranging from 0 to 89 (file attached here) and for TWI, I get values between 19.58 and 9.40 (file attached), with no data data for most of the area!
However, while extracting values for both slope and TWI to my point shape file, I mostly get 0 as the RASTERVALUE for most of the point data.

Could anyone kindly review  where exactly I am going wrong and how do I correctly get the slope and the TWI?


texts: Performing Sine inverse calculation in FME?

texts:
I am trying to find the formula for sine inverse calculation using Safe software FME.

Is there any function in ExpressionEvaluator?.


texts: Return all results within a 30km radius of a specific lat/long point?


texts: Is there a way to clip areas that are intersecting in Cartodb?

texts: I am new to cartodb and wish to know how can one clip only the areas that are intersecting?


texts: ArcMap - Data Driven Pages: Group features where page name field is the same

texts:
I'm using DDP to show a series of farmland areas.  Some of these have the same farmland area name but different attributes (based on land cover), so are two different features (by necessity). Is there a way to group these so that I can group these features for the DDP without having to edit the featureclass to create a multi-part feature, and so lose the detail I need to display in the map??

I'm using ArcGIS 10.0.


texts: Determine when I have crossed state line?

texts:
We are currently using the Google Maps API to get route information and total mileage but we would like to determine when we are crossing into a new state. When we return the route information we can be handed 2000+ Latitude and Longitude points. Initial thought was to do a lookup on each point to determine what state it was in but that would be very intensive. Another option would be to try and find a way to determine if the Lat & Long where inside a polygon or KMZ. Has anyone seen code for this?

Otherwise what other options should I be investigating?


texts: In reference to How to determine elevation drop along stream with ArcGIS Desktop?, suppose I want to cut the stream every 100 m and find the elevation drop every 100 m. How can I do this?


texts: OSM Line 'way' from PostgreSQL (PostGIS) to Java Linestring (from JTS library)

texts:
I would like to read the data from PostGIS OSM database to Java object. I have 'planet_osm_line' table. I'd like to get 'osm_id' and 'way' column read by ResultSet, as here:

stmt = con.createStatement();
ResultSet rs = stmt.executeQuery(query);
List<OSMLineData> osmLineDataList = new ArrayList<OSMLineData>();
while (rs.next()) {
    LineString ls = (LineString) rs.getObject("way");
    OSMLineData osmLineData = new OSMLineData(rs.getInt("osm_id"), (LineString) rs.getObject("way"));
    osmLineDataList.add(osmLineData);
}
rs.close();
stmt.close();
return osmLineDataList;

However, I get an ereor:

Caused by: java.lang.ClassCastException: org.postgresql.util.PGobject cannot be cast to com.vividsolutions.jts.geom.LineString

When I change it to read String rs.getString(way), and I also change the class variable type from Linestring to String, I am getting some weird representation of 'way' data.

I thought that it would be ok to import data for 'way' column with 'st_astext(way)' to the variable of String type, and then somehow change it to Linestring type from JTS. Is this a correct way? Is there a different, more efficient way to do this?


texts: How to add polygon height in QGIS or PostGIS?

texts:
Is there an easy way to set the height of polygons in QGIS or PostGIS? I am trying to visualize polygons by turning them into volumes (same shape at the base and roof, but with walls). I would then like to visualize the polygons in qgis2threejs.

Preferably the height would be just a column in the attribute table. Or a plugin that opens a window when you click on a polygon and then enter the height, which then gets transferred to PostGIS.

Clarification: My data currently does not contain any height information. Just normal flat PostGIS polygons.


texts:
I have pretty large dataset (>10M records) with a geography points (id, latitude and longitude), and need to find closest points from this set. Data itself is not updated very often but I need quick response. And I'm looking for something fast. Unfortunately I'm not experienced GIS user.

I tried to use PostGIS with this setup. Made this schema:

CREATE TABLE place (
   id int,
   name varchar(100),
   latitude float,
   longitude float,
   gpoint geometry
);
SELECT AddGeometryColumn ('public','place','gpoint',4326,'POINT',2);
UPDATE place SET geometry = ST_PointFromText('POINT(' || longitude || ' ' || latitude || ')', 4326);
CREATE INDEX idx_geo ON place USING gist(gpoint);

And tried to call data from database using this query:

SELECT * FROM place
ORDER BY gpoint <-> ST_GeogFromText('POINT(44.509648 40.184822)')
LIMIT 1;

But speed is not "close to realtime", about 30 seconds on my work laptop.

I've tried to make some optimizations: VACUUM ANALYZE, updated to latest Postgres and Postgis, use clustered table, may be something else (my stack already overloaded with all this tries). Last what I did I moved to another server with much more memory and space, but it increased speed only 2 times. Didn't tried yet to put everything in memory, but honestly I don't know yet how to make it right. Have other databases with a lot of data inside on this sever.

After all I'm not sure that I'm moving right direction. Do I use right storage and database. For some reasons I think that it is should be something trivial. But after few days searching I didn't found solution in the internet.

I want something fast. Under 1 second response time. Is it possible to get this numbers using Postgris or can you recommend to switch to something else? Since this is my pet project, so I want to use something with cheap cost of ownership.


texts: Transformation datum from WGS84 to BTM (Bangladesh Transverse Mercator)

texts:
I'm changing raster layers in ArcMap from WGS84 to BTM (Bangladesh Transverse Mercator). I first set my data frame to be in BTM
[BTM
Projection: Transverse_Mercator
False_Easting: 500000.0
False_Northing: -2000000000.0
Central_Meridian: 90.0
Scale_Factor: 0.9996
Latitude_Of_Origin: 0.0
Linear Unit: Meter (1.0)]
When I add the layers I'm asked to do a datum transformation apparently. I found two transformation datum which differ in ΔZ. Which one is accurate?
First one:
ΔX  -283.729m
ΔY  -735.942m
ΔZ    -261.143m
Second:
ΔX  -283.729m
ΔY  -735.942m
ΔZ    -161.143m


texts: Tabulating average value through spatial join in QGIS

texts:
Previous ArcGIS user and have now switched over to QGIS for a project. I am performing a spatial join of two shapefiles - one of census tracts and the other of zctas. In ArcGIS you can set the merge rule to average the values of data in one of the original layers by the other in the joined layer. Is this possible in QGIS?
When I perform the join, I am seeing values for the field originally in the census tract data in the zcta (larger) units -- but I'm assuming this is the first value found, right? Any clarification?


texts: Postgis Distance calculation on huge datasets

texts:
We have 7,512,694 records in our table,

Columns:

id, name, lat (numeric), lon (numeric), location (geometry(Point,4326))

Indexes we created:
create index on table using gist(location);

and of course rest of the columns are also indexed.

Requirement:
We need to get places that are in 100m distance of a given point. The queries I am using right now are taking more than 100s (yes seconds:()

These are the queries I have tried so far:

Query 1:

EXPLAIN ANALYZE
SELECT name FROM businesses
WHERE ST_Distance(location::geography, 'SRID=4326;POINT(-121.696176 38.15399)'::geography) < 100 ;

Result:

Seq Scan on businesses  (cost=0.00..2997214.91 rows=2504231 width=543) (actual time=5982.917..120110.083 rows=6 loops=1)
    Filter: (_st_distance((location)::geography, '0101000020E6100000A510C8258E6C5EC09DF4BEF1B5134340'::geography, '0'::double precision, true) < '100'::double precision)
    Rows Removed by Filter: 7512688
Planning time: 0.138 ms
Execution time: 120110.129 ms
(5 rows)

Query 2:

explain analyze SELECT name FROM businesses
WHERE ST_DWithin( ST_Transform(ST_GeomFromText('POINT(-121.696176 38.15399)',4326),26986), ST_Transform(location,26986), 100);

Result:

Sort  (cost=3104972.93..3105223.35 rows=100169 width=543) (actual time=120216.612..120216.619 rows=6 loops=1)
Sort Key: name
Sort Method: quicksort  Memory: 27kB
->  Seq Scan on businesses  (cost=0.00..3072341.85 rows=100169 width=543) (actual time=39426.948..120216.532 rows=6 loops=1)
    Filter: ((st_transform(location, 26986) && '01030000206A6900000100000005000000E758D42C93364EC16028183E1ADA3941E758D42C93364EC16028183EE2DA3941E758D42C2F364EC16028183EE2DA3941E758D42C2F364EC16028183E1ADA3941E758D42C93364EC16028183E1ADA3941'::geometry) AND ('01010000206A690000E758D42C61364EC16028183E7EDA3941'::geometry && st_expand(st_transform(location, 26986), '100'::double precision)) AND _st_dwithin('01010000206A690000E758D42C61364EC16028183E7EDA3941'::geometry, st_transform(location, 26986), '100'::double precision))
    Rows Removed by Filter: 7512688

What method can I use to decrease the query time?


texts: Get all vertices of a polygon using OGR and Python

texts:
I'm having a little trouble with the Python OGR API.  What I am trying to do is get all the coordinates of each vertex of the outer ring of a polygon.

This is what I have so far:

import osgeo.ogr
import glob

path = "/home/woo/maps/"
out = path + 'output.txt'

file = open(out,'w')
for filename in glob.glob(path + "*.shp"):
    ds = osgeo.ogr.Open(filename)
    layer1 = ds.GetLayer(0)
    print layer1.GetExtent()
    for feat in layer1:
        geom = feat.GetGeometryRef()
        ring = geom.GetGeometryRef(0)
        points = ring.GetPointCount()
        #Not sure what to do here

file.close()

I have heard that you can just for over the region but that only returns the rings in the polygon, not the nodes.

Anyone able to help.


texts: Copying point features along line features by predetermined distance

texts:
I have points in a feature class. Each point has its centroid on a unique line feature in a polyline feature class ( 62 points / 62 lines). Furthermore, points and their corresponding lines share unique identifiers in each of their attribute tables. My aim is create a new point feature along each of these lines at a distance (from old point to new point) also determined by a field in the line attributes.

Does this need to be scripted with python or is there a tool I can use?

I'm using ArcMap 10.2.2.


texts: Very simple Latitude/Longitude to WGS84 conversion

texts: As from here, I'm using the OpenTopography API to download a geoTiff heightmap, and it uses 4 WGS84 Coordinates, but the map I've supplied to the user gives back longitude/latitude coordinates. I'm very new to GIS, so I was wondering if there was some kind of equation for converting them.


texts:
I have a county-level map that is displaying nicely ... except for around 300 counties, all in 5-6 states, which are blank. All of the FIPs, county names, lats/longs are consistent in the data file, so I'm puzzled. Any clues?

Map: https://jpaulwright.carto.com/viz/5dcdc786-6183-11e6-ad04-0e3ff518bd15/public_map

Data file: http://www.economicmodeling.com/wp-content/uploads/2016/08/Talent-Map-v2d.xls


texts: GeoDjango - Use GDAL python binding in the shell

texts:
I would like to use the GDAL python binding in my GeoDjango site. The aim of this is to use ogr2ogr capabilities in a script. I am on windows 7 so I have installed GDAL-1.9.2.win32-py2.7 downloaded here

The installation went correctly. I can import GDAL in python with: from osgeo import GDAL

However, in my geodjango app python shell, I have an error when I try to import osgeo. The error is the following:

>>> import osgeo
Traceback (most recent call last):
  File "<console>", line 1, in <module>
  File "C:\Python27\lib\site-packages\osgeo\__init__.py", line 21, in
_gdal = swig_import_helper()
  File "C:\Python27\lib\site-packages\osgeo\__init__.py", line 17, in
t_helper_mod = imp.load_module('_gdal', fp, pathname, description)
 ImportError: DLL load failed: La procÚdure spÚcifiÚe est introuvable.
>>>

Can someone help me with this?

Im using GeoDjango v. 1.4.3 | python 2.7.3 | GDAL 1.9.2 win32


texts: Installing ArcSDE 10 on SQL Server 2008 R2 Cluster 64bit or 32bit?

texts:
I am looking for Tips an Tricks on installing ArcSDE on SQL Server 2008 R2 Cluster. Our IT division would like to have all the databases on their new 2008 cluster.

I've only done standard install before and I don't know how to go about it.

Does SDE have to be install on each of the cluster servers?

Does it have to be installed outside the cluster , on a seperate machine?


texts: Reset canvas to match print layout/ editing existing map layout QGIS 3

texts:
I have multiple map layouts from a large QGIS project which need editing and updating. The layouts are locked so they can be viewed as printed without updating to the current layers in the main QGIS window. Is there a way to reopen or load the layer from a map layout to view in the main QGIS map canvas window? The same question was poised here Reset canvas match print composer QGIS 3, however as Brian Fisher states the solution of using map themes works for new projects not existing ones where the different map themes have not been set up.
I have reposted the question as there may be new developments since 2 years ago. I am currently using QGIS 3.12.1.


texts: How to overlay a buffered layer to another layer in Google Earth Engine?

texts:
So i have this feature collection data, an .shp. I have buffered the polyline for a 20kms radii and tried to overlay (add) the buffered layer to another layer. but it can't be done. and anyway, I'm sorry if I perhaps missused the term 'overlay'.

Here's my attempt

var sesar = ee.FeatureCollection('users/sdwijati/Sesar_Jawa')
var areas = ee.FeatureCollection('users/sdwijati/IDN_adm2')
var sand = ee.Image("OpenLandMap/SOL/SOL_SAND-WFRACTION_USDA-3A1A1A_M/v02");

var bufferBy = function(size) {
  return function(feature) {
    return feature.buffer(size);
  };
};

var sarbuf = sesar.map(bufferBy(20000))

var nama_prov = ['Jawa Barat', 'Jakarta Raya', 'Jawa Tengah',
                  'Jawa Timur', 'Yogyakarta', 'Banten']
var area = areas.filter(ee.Filter.inList('NAME_1', nama_prov));

var sandavg = sand.reduce(ee.Reducer.mean());
var gtsandavg = sandavg.updateMask(sandavg.gt(50));

var xx = sarbuf.add(gtsandavg)

var sandvis= {min: 0, max: 100, palette: ['fbff00','b4ff00','00c4ff','0005ad']};
Map.addLayer(xx.clip(area), sandvis)
Map.centerObject(area, 7)

the code can't be run, the error message is 'sarbuf.add is not a function', my assumption is that because the variable sarbuf itself have no value. Any idea how to overlay the buffered layer?


texts: Import GeoTIFF into MapInfo - why "control points", and how to batch-create TAB files?

texts:
I have 20 GeoTIFFs (produced by gdal-merge of JPG/JGW-filepair hi-res aerial photographs) in total for this job; several hundred, overall for current and pending projects.

These files open without incident in QGIS and ArcGIS: drag'n'drop the file into a workspace, and both suites even detect the idiosyncratic deprecated projection (EPSG 3785 - not my choice) and put the images in the right spot on a map.

To try and open the identical file in MapInfo Pro 12.5 I amend MAPINFOW.PRJ to take account of it ignoring the (deprecated) existence of EPSG 3785, but still the 'Open Table' command (since when is a TIFF a 'table'?) acts as if I'm opening a file with no geodata at all. It wants me to pick three points on each image, type in the co-ordinates manually for each point, and then it will consider my request to add the image.

If MapInfo wants its 'TAB' file in order to perform the function of adding a geoTIFF raster layer to a map, how do I go about generating such a file (preferably in a batch)?


texts: Web request to upload shapefile to GeoServer

texts:
I have been trying to publish a zipped up shapefile to GeoServer.

For background, I am using DotSpatial to create a Line shapefile based on uploaded points. Once I have created the shapefile, I compress all files generated (shp, prj, dbf, shx). Then I want to publish it to GeoServer

I am able to create a Workspace and DataStore using a web request. When I tried adding the zipped up shapefile, all shapefiles are placed in the data\ folder.

No featuretype.xml or or layer.xml file is created under the workspace folder. The datastore.xml does get created.

The error I get back is:

  500 Server Error

It worked at some point but is now not working and I haven't changed the code.

  Dim localShapeFile() As Byte = readLocalShapeFile(zipUri) ' read file bytes

  Dim sUrl As String = "http://localhost:8080/geoserver/rest/workspaces/" & workspace & "/datastores/" & dsname & "/file.shp"

  Dim request As WebRequest = WebRequest.Create(sUrl)
  request.ContentType = "application/zip"
  request.Method = "PUT"
  request.Credentials = New NetworkCredential("admin", "geoserver")
  Dim requestStream As Stream = request.GetRequestStream
  requestStream.Write(localShapeFile, 0, localShapeFile.Length)
  requestStream.Close()
  Dim response As WebResponse = request.GetResponse ' error is produced here


texts: Find Shapefiles for world map, country borders and ISO3

texts:
I am trying to look for a source of clean Shapefiles in .shp format with country borders, for every country around the globe.

Would you be able to point me to an online free source?


texts: Is there a way to add x , y, z coordinate to the attribute table of polygon in ArcGIS 10.1?

texts: I have a large No. of converted DXF files, that represent the pitched roof. I want to calculate the roof area for each single building (that may consist of 4 to 10 polygon). The adjacent polygon must be identified, for that I need a x,y,z for each set of points to specified polygon.


texts: Appending multiple shapefiles into existing blank feature classes using ModelBuilder?

texts:
I have a set of shapefiles that I need to append into an existing blank database. The database I have contains blank copies of the shapefiles, only as feature classes in a feature dataset. Each feature class has the same geometry and schema as it's corresponding shapefile.

What I am trying to accomplish is to essentially load the data from my shapefiles into their corresponding feature class based on the feature class's name. My workflow thus far has been to manually load each feature using the loading wizard in catalog.

Is there a way to accomplish this with ModelBuilder?


texts: How to prevent small polygon WMS features to be hidden when zoomed out (GeoServer + OpenLayers)

texts:
Currently I'm developing a Web-GIS app where I use GeoServer 2.61 and OpenLayers 2. The map data is served from PostgreSQL+PostGIS via GeoServer (WMS).

I'm having a problem when I'm viewing polygon map. Some maps have a very little polygon feature, which do not show at certain zoom level when I zoom out (because the feature will be too small). Is there any workaround to prevent this to happen? I need to see ALL the polygon features at the lowest zoom level.


texts: Adding button zoom to map using OpenLayers 3?

texts:
I am a beginner in web programming using OpenLayers 3.

In my project i want to add a button to zoom in the map by clicking or selecting a part in the map using OpenLayers 3.


texts: Using continuous color band for raster legend in QGIS composer legend?

texts:
I'm looking for a way to insert continuous colour legend for a raster layer.
I managed to create discrete symbols with values along the bar, but I would like to have a continuous bar with some values on the bar, similar to the old plugin "One-Band Raster", that it doesn't work with the QGIS version higher than 1.7

I tried to set in properties->style, where I selected the  Render Type "Singleband pseudocolor" and see the "generate new color ramp" box.  I tried both "Continuous" or "Equal interval". However, in the Composer the legend always is displayed with discrete symbols.

Is there any way that can make QGIS Composer to display a nice continuous band for a raster layer?


texts: Find GPX files which contain tracks located in a specified bounding box

texts:
I have a bunch of GPX files and I want to find out which files contain tracks located in certain geographic extent, e.g. files containing tracks in Pennsylvania. Also, the inverse is also useful too, e.g. any file containing tracks not in USA. How to do this using any command line tool?

Gpsbabel's polygon filter looks like what I want. But for some reason I could not make it work. The filter does not filter anything. I am not sure whether the polygon filter works with track data. This is how I run it:

gpsbabel -i gpx -f infile.gpx -x polygon,file=test_bb.txt -o gpx -F outfile.gpx

Adding exclude suboption to polygon filter does not make any different to the output as well.


texts: How do I perform a proximity search with Postgis?

texts:
I downloaded the Geonames database (cities1000) and wrote a small Ruby program to import them into my table (geo_cities).  I then added geography column called geog.

Then I converted all of the lat/lon numbers into the geog column using:

update public.geo_cities set geog = st_GeogFromText('SRID=4326;POINT(' || longitude || ' ' || latitude || ')');

Things are looking good.  Now what I want to do is find all cities within 100 miles of Prague.

So I can get Prague like:

select * from geo_cities where asciiname = 'Prague' and countrycode = 'CZ';

I'm still learning GIS and Postgres so could someone help me with the simple query?


texts: QGIS plug-in equivalent of ArcGIS's Spatial Analyst Tools -> Overlay/ Weighted Overlay

texts: Specifically I have 12 .bil raster files that show monthly global mean minimum temperature. I file for each month. I want to combine these to create an ANNUAL global mean minimum temperature. In Arc I think that Spatial Analyst Tools-> Overlay -> Weighted overlay would be the way to go. Is there an equivalent plug-in for QGIS?


texts: Reading out raster data from multiple regions into xyz format in ArcGIS Desktop?

texts:
I‘m trying to export maize yield data from a raster to some sort of text file in order to analyze the data statistically. So far I‘ve managed to do so for an entire country (France), so that I ended up with a table containing the coordinates of about 8000 10x10km squares with the corresponding yield data (see pictures)

Now, data on country level is not precise enough for my purpose, I need to know for each square at least the region (department) it is located in. Extracting each department individually is way too tedious, especially considering that I will have to do the same for multiple other European countries.

Is there a way to easily extract the yield data from the raster so that I know in which department each square is located?


texts: Can I convert this exponential function in QGIS with Raster Calculator ?


texts: ArcGIS Surface Length Tool equivalent in FME?

texts:
I am trying to calculate the 3D length of a line.

I currently have a rather complex and cumbersome model built in ArcGIS ModelBuilder that uses a tool called "Surface Length".

I was wondering if there is an equivalent tool in FME?

Ultimately I want to calculate the 3D length of a shapefile. I have DEM data for the same area.


texts: Retrieve polygon attribute when positionning a point - openlayers / leaflet

texts:
Using a client-side webmapping library (openlayers3 or leaflet), I want my user to be able to add one marker on the map. I will then be able to retrieve long/lat coordinates.

From that, I want to know if it is located inside a polygon (from another layer) and retrieve polygon attributes. Is that possible ?


texts: Equal length dashed lines while using 'segments_to_lines' in QGIS

texts:
I want to create tapered lines with a dashed line stroke (to style intermittent streams).
I followed this tutorial to create the taperered lines, basically using Geometry generator with the expression segments_to_lines($geometry) and added the following stroke width expression:
CASE
  WHEN @geometry_part_num < 15 THEN  (0.222750 * (@geometry_part_num / 15))
  ELSE 0.22750
END

The result is as it is supposed to be:
But when I style the tapered lines with a dashed line stroke the dashes (and spaces) have very varied lengths (I assume due to segments_to_lines($geometry)):
For reference, this is how the dashed lines are supposed to look:
Does anybody have an idea how I could achieve this effect while having the line tapered?


texts: How to join multiple records from a table to a point feature

texts:
I work in ArcMap 10.4.
I have a xls table of 320 records. Some of the records are of the same name (municipality name).
I have a feature point layer of all municipalities in the country. Both the xls table and feature point (.shp) layer have a common ID field - so they are "joinable" and "relatable".
What I need to do is to create a feature point layer for all the 320 records, so they have a spatial attribute too (a point).

So, to illustrate it better:

ID, Name, POI (.xls)

1, Town_A, Castle

1, Town_A, Church

1, Town_A, Theatre

2, Town_B, Castle

2, Town_B, Monument

3, Town_C, Castle

...

ID, Feature point (.shp)

1, Town_A

2, Town_B

3, Town_C

...

I tried "Join" which didn't work. It kept only one record for one point. "Relate" would work by principle (one-to-many) but I can't export it, it won't append the data.
I hope I making myself clear. I need something which would multiple (copy) the point features by how many records there are in the .xls table. So in the result, there would be many overlapping points (for the cases where they are multiple records for one municipality).


texts:
I have a table containing thousands of data points around London.

I'm trying to write a query to find the closest data points to a given latitude and longitude.

I'm getting different results based on whether I order by the <-> distance operator or the result of ST_DistanceSphere. Here are the two queries and their results:

Order By <->
------------
SELECT
uid,
ST_DistanceSphere(
  ST_PointFromText('POINT(-0.186243362115208 51.5128590051459)', 4326),
  the_geom
) as d
FROM points
ORDER BY the_geom <-> ST_SetSRID(ST_Point(-0.186243362115208, 51.5128590051459), 4326) LIMIT 5;

Results:

  uid   |         st_astext          |       d
---------+----------------------------+---------------
6997685 | POINT(-0.186723 51.500052) | 1424.46281677
6997684 | POINT(-0.187011 51.500056) | 1424.62219563
6997686 | POINT(-0.186435 51.500047) | 1424.69367394
6997683 | POINT(-0.187299 51.500061) | 1424.94947768
6997687 | POINT(-0.186147 51.500043) |    1425.09232

Order By ST_DistanceSphere
--------------------------
SELECT
uid,
ST_DistanceSphere(
  ST_PointFromText('POINT(-0.186243362115208 51.5128590051459)', 4326),
  the_geom
) as d
FROM points
ORDER BY d LIMIT 5;

Results:

   uid   |         st_astext          |       d
---------+----------------------------+---------------
6997686 | POINT(-0.186435 51.500047) | 1424.69367394
6997685 | POINT(-0.186723 51.500052) | 1424.46281677
6997687 | POINT(-0.186147 51.500043) |    1425.09232
6997688 | POINT(-0.185859 51.500039) | 1425.76968591
6997684 | POINT(-0.187011 51.500056) | 1424.62219563

As you can see, the first query returns 6997683 in the results, but the second query does not (and returns 6997688 instead). 6997683 is the closer value of the two, so 6997688 should not be returned by either query.

The type of the_geom is ST_Point and the SRID (checked with Find_SRID) is 4326.

What might be causing the difference, and what am I doing wrong?


texts: New to ArcGIS and want to compute distance between two points "as the fish swims"

texts:
We have a layer named Hydrography which is the water area for a given assignment.
I have a station point and a target point.  I need to compute the distance between those two points but inside of the hydrography.  (boats cannot go over land)

The DefenseService distance computation looks passable but it seems to be in a straight line since it does not take any Raster as input.  The Euclidean seems to almost be there but I cannot find anything that is straight forward.

I am using ArcGIS Explorer 2500, ArcGIS for Server 10.1 SP1 Enterprise, and we are writing our services and tools in C#.

I have the loading of the hydrography, the creation of the IPoints, but at a loss as to what to do next!

Any help appreciated, please.


texts: Select features of a shape which are not in another shape

texts:
I have a shape with some features (let's say 100) and a second shape with some of this features (50) buffered (5 meters inside buffer). I want to select and write in a new shape every feature of the first shape which doesn't contain the buffered feature of the second shape. I tried two different approaches with no success:

1) Select by location. Select the features of shape 1 which doesn't intersect with shape2 (ST_difference):

ogr2ogr diff.shp parcels1.shp -dialect SQLite -sql "select st_difference(a.geometry, b.geometry) as geometry from parcels1 as a, """parcels1_buffer5m.shp""".parcels1_buffer5m as b"

No error, but every feature is selected.

2) Select ID (in my case DN_OID) which is not included in the list of DN_OID of second shape:

import fiona
with fiona.open("parcels1.shp") as input:
    meta = input.meta
    with fiona.open('diff.shp', 'w',**meta) as output:
        for feature in input:
             if feature['properties']['DN_OID'] not in ("SELECT DN_OID FROM parcels1_buffer5m.geom"):
                 output.write(feature)

Any advice is welcomed.


texts: Why NQPD is better than inverse of MAD/SAD in closeness measurement?

texts:
sDNA recommend use NQPD to measure closeness rather than inverse SAD/SMD/SCD/MAD/MMD/MCD.  (Please refer the 3.b Closeness centrality sDNA measure description documentation )

In my understanding, their computation are as the following(for example with angular):
- NQPDA: sum of (each link 1/distance)
- Inverse SAD: 1 / (sum of each link distance)
- Inverse MAD: link_count / (sum of each link distance)

The major difference between them are get the inverse first or get the sum first.
But why the NQPD is better? Is there any explanation about this? maybe from mathematics?


texts: Layout view graph display not filling frame in ArcMap

texts:
I'm using ArcMap 10.6.1 and when I create a graph and add it to Layout view it doesn't fill the frame, but only a portion of it. The rest of the graph frame is black, as in the screenshot below.

Is this a setting issue, or is there another way to fix this?


texts: Does ora2pg support the SDO_GEOMETRY spatial type?

texts:
I've been searching the net up and down for days now and can't find a definite answer to this question, neither on google nor on the official page (http://ora2pg.darold.net/).

Does ora2pg support Oracle's SDO_GEOMETRY spatial type?

If I run ora2pg on an Oracle test table I get an output like

CREATE TABLE platzerl2 (
        objectid bigint NOT NULL,
        name varchar(40),
        class varchar(25),
        st varchar(2),
        stfips varchar(2),
        placefip varchar(5),
        houseunits bigint,
        pop2000 bigint,
        pop_class integer,
        arealand double precision,
        areawater double precision,
        shape SDO_GEOMETRY,
        se_anno_cad_data bytea
);

Clearly SDO_GEOMETRY should read something like GEOMETRY, Postgres/PostGIS's spatial type.

So, is SDO_GEOMETRY unsupported or am I doing something wrong?


texts: Calculating Topographic Wetness Index (choosing from different algorithms)

texts:
Topographic Wetness Index can be expressed as

 Ln(a/tanB) based on the idea of Beven and Kirkby (1979)

where

  a is the specific catchment area (a=A/L, catchment area (A)divided by contour length(L))

and

  tanB is the slope

The basic idea here is simple, but as there are multiple ways to calculate both a and tanB, the results of a TWI can vary widely (Qin et al. 2011).

Flow accumulation and catchment area can be calculated, for example by:

 D8 (O'Callaghan, J.F. / Mark, D.M. (1984))
 D-infinity  (Tarboton, D.G. (1997)
 Triangular Multiple flow direction (Seibert, J. / McGlynn, B. (2007)

algorithms, and there are many other algorithms available too.

Slope is usually calculated as local slope around the pixel (Sorensen et al. 2005).
Local slope can also be calculated as minimum, mean and maximum slope around the pixel. Another way to calculate slope is presented by Hjerdt et al. 2004 where slope is calculated to a point d meters below cell center.

Slope is a basic tool in most of the GIS softwares, however the calculation can differ.
Here are few examples:
ESRI: http://webhelp.esri.com/arcgisdesktop/9.2/index.cfm?TopicName=Calculating_slope
SAGA: http://sourceforge.net/apps/trac/saga-gis/wiki/Terrain%20Analysis%20-%20Morphometry%20module%20library

As you can see there are many options available to calculate both a and tanB.
So, the question is, in practice, which is the proper (best) way to calculate TWI using different these algorithms? Or is there any?

I personally like working in SAGA, mainly because there are a large selection of open source hydrology tools.

P.s. I'm having a hard time to find out exactly how catchment slope is calculated in Saga GIS, and exactly what does it mean here.  (Terrain analysis -hydrology: catchment area parallel).

EDITED:
Answered by Volker Wichmann from SAGA Forums:
"The catchment slope output grid of the Catchment Area (Parallel) module is computed like this: for each cell, the local slope is calculated using the approach of Zevenbergen & Thorne. These slope values are accumulated downslope. Finally, for each cell the accumulated slope values are divided by the derived catchment area of the cell. The unit of the grid are radians."

"The Topographic Wetness Index (TWI) module requires a normal slope grid as input. "

References:

Beven and Kirkby 1979. A physically based variable contributing area
   model of basin hydrology. Hydrological Sciences Bulletin, 24, pp.
   43–69.

Hjerdt et al. 2004. A new topographic index to quantify downslope
   controls on local drainage. Water Resources Research, 40, W05602,
   doi:10.1029/2004WR003130.

O'Callaghan, J.F. and Mark, D.M. 1984. The extraction of drainage
   networks from digital elevation data.Computer Vision, Graphics and
   Image Processing, 28:323-344

Qin et al. 2011. An approach to computing topographic wetness index
   based on maximum downslope gradient. Precision Agric  12:32–43.

Seibert, J. and McGlynn, B. 2007. A new triangular multiple flow
   direction algorithm for computing upslope areas from gridded digital
   elevation models,
   Water Ressources Research, Vol. 43, W04501

Sorensen et al. 2005. On the calculation of the topographic wetness
   index: evaluation of different methods based on field observations.
   Hydrol. Earth Sys. Sci. Discuss., 2, 1807–1834

Tarboton, D.G. 1997. A new method for the determination of flow
   directions and upslope areas in grid digital elevation models, Water
   Ressources Research, Vol.33, No.2, p.309-319


texts:
I have a Cartosat 2.5 m Imagery and Block files along with DSM Data. I wanted to generate DEM using Leica Photogrammetry Suite Software. Does anyone know the steps of generating same?

I  have to edit the existing DSM by augmenting with more breaklines / masspoints and editing of existing points to the exact ground surface as visualized on the stereo data supplied.

Also is DSM also in black/white color or is there a way to change bands and update the colors?

How can I perform the same on LPS Software?


texts: Getting error while running Longest Flow Path for catchment

texts:
I am getting the following error while executing the Longest Flow Path for Catchments in Arc Hydro Tools. The tool requires a catchment shapefile which is in float and other file is flow direction (not in float). Both of these files were created using the HEC and ArcHydro tools.

Any way to solve this problem?

  System.Runtime.InteropServices.COMException (ox80041098): ERROR 010236: Error in building VAT.
  ERROR 010067: Error in executing grid expression
  at
  ESRI.ARCGIS>GeoAnalyst.RasterCOnversionOpClass.ToRasterDataset(IGeoDataset dataset, String rasterFormat, IWorkspace pWorkspace, String name)
  at
  ESRI.APWR.ApHydro.ArcHydroOp.LongestFlowPathforCatchment(ApLayers aplayers, Int32 debug, ITrackCancel& trackcancel, IGPMessages&messages, Strings& exMessage) in
  CBuilds\9s]ArcHydroTools\src\ApHydroOp.vb:line 68183


texts: Points to shapefile, with curved lines given center and radius (QGIS)

texts:
I'm a newbie in GIS and for a project I need to create a shapefile given some points in lat/lon. I have been following this tutorial: https://community.pyramidanalytics.com/t/q5tk8w/creating-a-shapefile-using-latitude-and-longitude-coordinates and it's OK for polygons but I need some lines to be curved.

My data is described given vertices and circles, like this: 3400N 02500W; 3000N 02000W; 3139N 01725W from this point following an arc of a circle 100 NM centred on 3304N 01621W up to 3130N 01545W... and then it continue giving the points of the vertices of the polygon.

(There are other curved parts as well, described like the above part, but I think with one example I'd be able to do the others).

How can I solve this?


texts:
I have a CSV-file with "ID"s. I would like to identify a feature in my database table that matches those "ID"s.

How can I use a tester/feature merger or something similar to get the records I am looking for?


texts: I clipped a big DEM to a smaller one to reduce the DEM to the area I need. Now i want to show the different heights in different colors. I did the contour lines. Due to the clipping, the contour lines will give me, using the process of "Line to Polygon", some polygons which have a flat side. I know this is because the contour lines are not closed. So how do i close them without doing it by hand following the boundaries of my clipping mask?


texts:
We are using the metes and bounds in old pipeline alignment sheets in an edit session in ArcMap to plug in the "direction/length" values to generate the line segments.

Does anyone have any idea how I could explore further on how would this be done in Python?

For example if I used an insert cursor and some how mangled that with some trigonometry... (grasping for straws here).  All I have is an alignment sheet giving a string of these: S-34-47-00-E 567'.


texts: QGIS: Export labels from layer to shapefile / Geojson

texts: I am having trouble finding a way to export only the labels of a layer in QGIS to either shapefile or geojson format.  I also tried to find a way to export both the geographic feature plus the label to shapefile or geojson format, but every time I get only the feature itself.  My ultimate goal is to have the labels in geojson format for web mapping purposes, since I already have a basemap but not the corresponding labels for some of the features.


texts:
I have a number of DEMs from the National Elevation Dataset (NED), that come standard in the NAD83 coordinate system.  I need to project them in order to define linear units for analysis.  Every time I project the DEM, the projected result contains a systematic pattern of error, where elevation values are distored up to 250 meters in either direction (higher/lower).  The errors appear to occur on a horizontal scale of 200 meters or so.  Not only are the errors apparent when calculating the difference in elevation between the projected and unprojected DEM, but they are obvious even just on visual inspection of the DEM and particularly in the hillshade.  I tried two different projections (Transverse Mercator and UTM Zone 12), both based on NAD83, and got similar errors but in different patterns.  I have attached a few images in hopes of illustrating what I may not be describing very well in words.  These patterns of error occur regardless of the resampling technique used when projecting the raster.

I know that prjecting a DEM inherently introduces some error by "stretching" but I have not noticed errors like these in projecting other DEMs. I am hesitatnt to use these projected DEMs in terrain analysis as the results don't look right even to the naked eye.


texts:
I have a script that looks at feature class inside a .gdb and prints some simple information about the attributes field. Note- these is actually only one feature class in the .gdb. I want use the arcpy.Describe method with fieldInfo. Every time the script gets to the line field_info = desc.fieldInfo, it throws the error AttributeError: DescribeData: Method fieldInfo does not exist. My feature class definitely has an attribute table with fields.

What is the problem here; why is this error being thrown and how can I fix it?

My end goal is to just get a list of the field names in my feature class attributes table:

import arcpy, sys, os

arcpy.env.workspace = r"C:/Workspace/Sandbox/MapChangeProject/data/Alabama.gdb"

fcList = arcpy.ListFeatureClasses()
for fc in fcList:
    print("FeatureClass: {}").format(fc)
    desc = arcpy.Describe(fc)

    field_info = desc.fieldInfo
    for index in range(0, field_info.count):
        print("Field Name: {0}".format(field_info.getFieldName(index)))


texts: Personal geodatabase items have an exclamation (!) in ArcMap

texts:
A user in our department is having a problem accessing personal geodatabase data in ArcCatalog and ArcMap.  Receiving error message:

Failed to connect to database.  An underlying database error occurred.
. . . you need permission to view its data.

Some background info:  This error is only occurring for a new user that has recently been given GIS rights by out IT department.  All other users can access this personal geodatabase data that the new user cannot.  I've consulted with our IT department and have been assured that the new user has the same rights as all other users.

So, all that said I am wondering if there is a permissions lock in the gis environment that I have to address.


texts:
I am interested to download ESRI Bump Map Tools for ArcGIS referred in the below link,

https://blogs.esri.com/esri/arcgis/2010/01/21/introducing-the-arcgis-bump-map-tools/

I have searched for these tools but cant find the working link to download and utilize these tools in ArcGIS 10.5.

If these tools are compatible with ArcGIS 10.5, where can I download them?


texts: Create automatically temperature or isotherms map with PostGIS?

texts:
Is there any possibility to create automatically temperature or isotherms map through a points table stored in PostGIS?
I need to serve this map on a web (Geoserver). The information is updated each hour and I need to have the map updated as well.


texts: How can I increase limits of OSRM table function?

texts:
I want to generate a distance matrix for more than 50 points with my own OSRM server. Because I want to minimize the number of calls, I decided to increase the size of the query from 50 to 100 or even more.

This is the answer of OSRM with 100 loc

{
  status_message: "Number of bearings does not match number of coordinates."
}

Query:

http:/.../table?loc=40.29864,-3.64438&loc=40.3986148,-3.6975552&...

I imagine there are limits, but is it possible to increase this limit?


texts:
I need an step-by-step Introduction to build in an GeoCoder by ESRI. I have downloaded the Story Map - Text and Legend - Template and now I am trying to build in a simple geocoder Services. But Anything works. The readme.html have an false Introduction for that. So I Need a simple step-by-step introduction please.

In German:

  Also ich brauche eine Anleitung um einen Geocoder-Service in meine
  WebGIS-Seite einzubauen. Ich habe von ESRI folgendes Template
  runtergeladen und im Aussehen geändert: Story Map - Text and Legend
  Ich habe bisher versucht mit anderen Tutorials das widget einzubauen,
  aber bisher ohne Erfolg. Ich wäre für eine Schritt-für-Schritt
  Anleitung sehr dankbar. In der ReadMe-Datei ist die Anleitung, schätze
  ich mal , veraltet. Und daher nutzlos


texts: How to convert shapefile of India into a system with longitude and latitude system?

texts:
I have two shapefiles of India. One is without PoK and other is with Pok. I want to convert the shapefile of India with PoK into a co-ordinate system with longitude and latitude (same as without PoK). I have tried in R, but could not crack the problem.

Here are the files for Without PoK and With PoK India?


texts:
I format my notebook because this problem, but it wasn't solved. Any suggestion to fix it?

The following packages have unmet dependencies:
 python-qgis : Depends: python-qgis-common (= 1:2.14.4-1+20trusty-ubuntugis) but it is not going to be installed
               Depends: libqgispython2.14.4 but it is not going to be installed
               Depends: libqgis-analysis2.14.4 but it is not going to be installed
               Depends: libqgis-core2.14.4 but it is not going to be installed
               Depends: libqgis-gui2.14.4 but it is not going to be installed
               Depends: libqgis-networkanalysis2.14.4 but it is not going to be installed
               Depends: libqgis-server2.14.4 but it is not going to be installed
 qgis : Depends: gdal-abi-2-1-0 but it is not installable
        Depends: libgdal20 (>= 2.0.1) but it is not installable
        Depends: libgeos-c1v5 (>= 3.4.2) but it is not installable
        Depends: libqgis-analysis2.14.4 but it is not going to be installed
        Depends: libqgis-app2.14.4 but it is not going to be installed
        Depends: libqgis-core2.14.4 but it is not going to be installed
        Depends: libqgis-gui2.14.4 but it is not going to be installed
        Depends: libqgis-networkanalysis2.14.4 but it is not going to be installed
        Depends: qgis-providers (= 1:2.14.4-1+20trusty-ubuntugis) but it is not going to be installed
        Recommends: qgis-plugin-globe but it is not going to be installed
        Recommends: qgis-provider-grass but it is not going to be installed
 qgis-plugin-grass : Depends: qgis-provider-grass (= 1:2.14.4-1+20trusty-ubuntugis) but it is not going to be installed
                     Depends: libgdal20 (>= 1.8.0) but it is not installable
                     Depends: libqgis-app2.14.4 but it is not going to be installed
                     Depends: libqgis-core2.14.4 but it is not going to be installed
                     Depends: libqgis-gui2.14.4 but it is not going to be installed
                     Depends: libqgisgrass7-2.14.4 but it is not going to be installed
E: Unable to correct problems, you have held broken packages.


texts: Symbolizing a line feature based on polygon location

texts:
using 10.2.2 desktop to be imported as a service on to AGOL

I am looking for a way to change symbology of my line features, water mains, based on a polygon, pressure zone. Field operations would like to see the color of water mains change as it flows through different pressure zones. I was hoping for a way where I could symbolize water mains to match the color of the known pressure zone without cutting up my line feature classes per polygon.


texts:
I'm not able to install pgrouting on Centos 6.4 with Postgresql 9,3 and postgis 2.0.
I tried from the rpm, but I didn't found the repository.
I tried building from the source but I got the following error:

# cmake -DPOSTGRESQL_EXECUTABLE:FILEPATH=/usr/pgsql-9.3/bin/postgres -DPOSTGRESQL_PG_CONFIG:FILEPATH=/usr/pgsql-9.3/bin/pg_config -L ..

[...]

    Boost  found.
Found Boost components:
   thread
-- CGAL_INCLUDE_DIR=CGAL_INCLUDE_DIR-NOTFOUND
-- CGAL_LIBRARIES=CGAL_LIBRARIES-NOTFOUND
-- BOOST_THREAD_LIBRARIES=optimized;boost_thread-mt-shared;debug;boost_thread-mt-shared-debug
-- GMP_LIBRARIES=GMP_LIBRARIES-NOTFOUND
-- CGAL not found.
CMake Error at CMakeLists.txt:294 (message):
   Please check your CGAL installation, or set correct path to CGAL_INCLUDE_DIR and CGAL_LIBRARIES.

-- Configuring incomplete, errors occurred!

I don't need CGAL, so How can I proceed without CGAL or How can I install it in order to proceed further?


texts: Distance matrix batch process for different area types in one layer

texts: I have a layer of the centroids of all neighbourhoods in a city, with neighbourhood ID and borough ID as attributes of each point. The city is made from boroughs, which are themselves made from neighbourhoods. I want to create a distance matrix between all neighbourhoods within each borough, but there are so many boroughs that selecting each and finding the distance matrix would take a long time. Is there a way to run this as a batch process? The batch process option seems to act between different layers.


texts: QGIS Print Composer Grid Draw Coordinates line offset

texts: I have had a good google about this and can't find anything on it. is it possible to offset  grid labels from the grid basically trying to move the label to where the red lines are on this crude example.


texts: GeoJSON - Resymbolize without redownloading?

texts:
GeoJSON served via GeoServer to OpenLayers.

Simple question - Is it possible to apply a new styleMap without redownloading the data from the server?

I am currently using the following, but am redownloading.

layer.styleMap = new_style
layer.refresh({force:true})


texts:
I have created a sqlite point layer in a QGIS (ver 3.8.1) project with both the absolute and relative path to photos stored in a subfolder of the project. I've added an action to the layer, which works great. In the image below, "Path" is the fieldname pointing to the absolute path

I have also created an attribute form, with 2 text fields and the relative path field using the drag-and-drop designer. The widget type is set to "Attachment", and the integrated document viewer is set to image (you can't see it in the screenshot). This results in an attribute form within QGIS with the image embedded within it.

All of this works great in QGIS, and mostly in QFIELD. The difference being that when editing in QField I am able to take a photo and see it in the attribute form, and when I hit save it claims to save the record but in fact it doesn't. The image is saved (in a subfolder of the project called DCIM, which is not configurable) but no record is created.

I found a reference from a couple years ago which stated that the project had to be stored in the Android/data/ch.opengis.qfield/ folder, which I have tried on the SD card and on the device memory. Device is a Samsung S8 with lots of available storage space.

There are several bug reports on GitHub which are the exact opposite, stating that the record is created but the image isn't saved (QField issue 571).

The attribute form designer is goofy, in that I can't seem to edit the fields after creating them. I have to delete and recreate, which makes it hard to confirm that I've tried every option, but I really think that I have. Given that others are having a different problem I feel that maybe my problem lies there though...


texts: How to separate one column with latitude and longitude into two columns. Lat/long are not separated by anything

texts:
I'm a GIS novice.  I have a CSV file that contains the latitude and longitude data for a small region in the northeast U.S. The data appear in one column appearing as one uninterrupted string of text in each field. I need to separate this data into two columns.
Each field is formatted like this, with all fields in this column having the same number of digits:
   xy4285021070865890

Two other issues: each field begins with "xy," and the longitude (or "y") should have a minus sign in front of it because the location for this map is within the western hemisphere.
How can I separate this column into two columns so that 1) the "xy" disappears and 2) the longitude includes a minus sign? OR will QGIS3 be able to automatically read and correctly interpret this column because of placement of "xy" before each number?
I've searched answers to this "one-column to two-column" issue, but all the questions referred to lat-long figures separated by a comma in their respective fields. There were some references to use of the field calculator - but I don't know how to use that, either.


texts: Shortening line geometry using expression dialog in QGIS?

texts:
I use the "Offset Point Symbol" tool and a geometry generator symbology to create dynamic leader lines. A sample file is available at http://hub.qgis.org/issues/15201 (Add snapping for 'Offset Point Symbol' tool).

Because I need a circle symbology without fill in some projects, I want to shorten the leader lines by the radius of the cirle symbology.

Can anyone help me with the expression? Symbology units are in meter ('Map units').

make_line (
make_point(
$x + regexp_substr("offset", '([^,]*)'),
$y - regexp_substr("offset", '([^,]*$)')
),
$geometry
)

Based on the answer posted by underdark I've created an expression that works properly:

-- radius = 2m

difference(
    make_line (
        make_point(
            $x + regexp_substr("offset", '([^,]*)'),
            $y - regexp_substr("offset", '([^,]*$)')
        ),
        $geometry
    ),
    buffer(make_point(
            $x + regexp_substr("offset", '([^,]*)'),
            $y - regexp_substr("offset", '([^,]*$)')
        ), 2)
  )


texts: Arc Hydro: Assign Stream Slope Producing _Mostly_ Correct Elevation

texts:
I'm pretty new to the Arc Hydro Toolbox.

I used assign stream slope in the Arc Hydro toolbox. It produced a mostly correct to and from elevation. However, if you look at the purple arrow, it shows a mistake. That river goes from low elevation to high and then back down to low (the highlighted main river channel is low again). It does it again at other locations above the purple circle.

Any suggestions on how to fix this? Are there any issues with a river line file that could trip up the tool?


texts:
I use OpenLayers 4.6 and I want to read a feature from GML string. The GML looks like:

<?xml version="1.0" encoding="UTF-8"?>
<gml:Polygon xmlns:gml="http://www.opengis.net/gml"
    xmlns:sch="http://www.ascc.net/xml/schematron"
    xmlns:xlink="http://www.w3.org/1999/xlink"
    xmlns:xs="http://www.w3.org/2001/XMLSchema" srsDimension="2">
   <gml:exterior>
      <gml:LinearRing srsDimension="2">
         <gml:posList>37.15 36.15 37.15 36.3 37.5 36.3 37.5 36.15 37.15 36.15</gml:posList>
      </gml:LinearRing>
   </gml:exterior>
</gml:Polygon>

I know that I need to specify the featureNS and featureType to read the feature using ol.format.GML3, but I don't know how to set them properly:

let feature = new ol.format.GML3({
  srsName: 'EPSG:4326',
  featureType: '???',
  featureNS: '???'
}).readFeatures(data);


texts: GeoServer: Best way to publish 2500 TIFF or 71 ECW files?

texts:
I have ortophotos of an area that I need to serve as a background map for my vectors. I got it as 2500 files (71,5MB each) in raw TIFF format with corresponding TFW word file - 180GB of data. Coordinate system is local and it matches with my vectors (there is no EPSG code for it but I named it as “32805” and put a proper definition).

For desktop use in MapInfo I converted them to ECW (with some tool that comes with MapInfo) and retile as much bigger to have only 71 file because opening 2500 tiff files is overkill. I just merged 49 TIFF’s (7x7) into one ECW - 35000x35000pixels – largest is around 200MB) It works great and its really fast in MapInfo.

Now I am quite confused - HOW to serve them in GeoServer ?

I have publish one TIFF and one ECW for comparison. ECW is much faster in browser preview (I am aware of ECW server licensing issue but that shouldn’t be a problem). I have found a presentation “GeoServer on steroids” and read about ImageMosaic, ImagePyramid, retiling, adding overview etc it was quite informative but still don’t know what to do.

My question is : how should I do it ? Mosaic or Pyramid, and if is one of the answers positive, I need your advice or some hints. I would really like to be ECW because of disk space so there will be no need to keep 180GB of tiff’s on a server.

Data will be served through LAN with maximum 20 users connected at peak hour.SQL server have not so huge amount of data. Sorry if I miss some other info, but I will send it if needed.

  Geoserver 2.1.4,  Windows 7 32bit,  2GB System Memory,  (1.7.0_09
  (Java HotSpot(TM) Server VM), Native JAI+Native JAI ImageIO = true

Original TIFF
gdalinfo D:\75720-47970.tif
Driver: GTiff/GeoTIFF
Files: D:\75720-47970.tif
       D:\75720-47970.tfw
Size is 5000, 5000
Coordinate System is `'
Origin = (7572000.000000000000000,4797500.000000000000000)
Pixel Size = (0.100000000000000,-0.100000000000000)
Metadata:
  TIFFTAG_SOFTWARE=Adobe Photoshop 7.0
  TIFFTAG_DATETIME=2006:10:09 13:02:57
  TIFFTAG_XRESOLUTION=72
  TIFFTAG_YRESOLUTION=72
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  ( 7572000.000, 4797500.000)
Lower Left  ( 7572000.000, 4797000.000)
Upper Right ( 7572500.000, 4797500.000)
Lower Right ( 7572500.000, 4797000.000)
Center      ( 7572250.000, 4797250.000)
Band 1 Block=5000x1 Type=Byte, ColorInterp=Red
Band 2 Block=5000x1 Type=Byte, ColorInterp=Green
Band 3 Block=5000x1 Type=Byte, ColorInterp=Blue
-------------

ECW file which is retiled from 7x7  original tiffs

gdalinfo D:\OF-45.ecw
Driver: ECW/ERDAS Compressed Wavelets (SDK 3.x)
Files: D:\OF-45.ecw
Size is 35000, 35000
Coordinate System is:
LOCAL_CS["LOCAL - (unsupported)",
    UNIT["Meter",1]]
Origin = (7571500.000000000000000,4798500.000000000000000)
Pixel Size = (0.100000000000000,-0.100000000000000)
Corner Coordinates:
Upper Left  ( 7571500.000, 4798500.000)
Lower Left  ( 7571500.000, 4795000.000)
Upper Right ( 7575000.000, 4798500.000)
Lower Right ( 7575000.000, 4795000.000)
Center      ( 7573250.000, 4796750.000)
Band 1 Block=35000x1 Type=Byte, ColorInterp=Red

  Overviews: 17500x17500, 8750x8750, 4375x4375, 2187x2187, 1093x1093, 546x546, 273x273, 136x136
Band 2 Block=35000x1 Type=Byte, ColorInterp=Green
  Overviews: 17500x17500, 8750x8750, 4375x4375, 2187x2187, 1093x1093, 546x546, 273x273, 136x136
Band 3 Block=35000x1 Type=Byte, ColorInterp=Blue
  Overviews: 17500x17500, 8750x8750, 4375x4375, 2187x2187, 1093x1093, 546x546, 273x273, 136x136


texts: Securing ArcGIS Server 10.x On .NET - User Store Point At LDAP

texts:
I have the fun job of securing ArcGIS Server, and one of the requirements is to use LDAP for authentication and authorization for services...

For ArcGIS Server installed on Java backend it looks to be a very simple task - LDAP is an option under User Store section.  It does not appear to be as straight-forward for .NET backend.

I know that I need to set up a custom ASP.NET membership provider in the ArcGIS Manager applications config file, but thats about all I know.  I've thoroughly searched ESRI's site for any type of "how-to" with no luck.

I am hoping someone can point me to a good link or give me any other useful info on how to accomplish this.


texts: Where can I find the SQL queries for use with OpenStreetMap-Carto?

texts: I found the OpenStreetMap-carto stylesheets (.mssfiles) here. I also found an issue discussing a problem with a SQL query here.  Does anyone know where I can find the SQL scripts to build the layers in TileMill? Maybe a TileMill project?


texts: How to measure quality of coordinate conversion?

texts:
I'm using PROJ.4 in my work and now I want to know how it convert coordinate well, especially longitude/latitude <-> UTM.

I can compare the result from proj command with values from some websites which provide coordinate conversion, but I don't think I can convince my clients or my supervisors with such unofficial data.

Is there any official or authorized reference to test coordinate conversion?


texts:
Yesterday I set up a Python script using a cursor on a layer.

Today I've tried running again and no matter what I run it in, the cursor is empty. I've tried removing the query on the layer, tried deleting the cursors and restarted my computer, but iterating through the rows in my table still does nothing. I've tried it directly on the FC (in SDE) and on the layer.

testcursor = arcpy.SearchCursor("testlayer")

Surely this should contain as many fows as the attribute table? Or have I got it wrong?

for row in testcursor:
    print "row exists"
    print row.getValue("SURVEY_NAME")
    survey = row.getValue("SURVEY_NAME")
    print survey
    if survey not in layerlist:
         layerlist.append(row.getValue("SURVEY_NAME"))

Yesterday this did something, but now it does nothing.

layerlist is empty

What have I done wrong?


texts: What is the fastest way to identify the two nearest linestrings using PostGIS?

texts:
I want to be able to identify the closest roads to two points that are generated through geocoding (a start and end address). The points will likely be near a linestring. Two approaches that I have been experimenting with are:

-- Method 1:
SELECT
    id
FROM
    routing
WHERE
    ST_DWithin(ST_Transform(ST_Setsrid(ST_Makepoint(-6.17706680000, 53.300791),4326), 900913), geom, 50.0)
ORDER BY ,
    ST_distance(geom_way, ST_Transform(ST_Setsrid(ST_Makepoint(-6.17706680000, 53.300791),4326), 900913))
LIMIT 1;
--Time: 228  ms

-- Method 2:
SELECT
        id
FROM
        routing
WHERE
        ST_Intersects(geom, ST_buffer(ST_Transform(ST_Setsrid(ST_Makepoint(-6.17706680000, 53.300791),4326), 900913), 50))
ORDER BY ,
        ST_Distance(geom, ST_Transform(ST_Setsrid(ST_Makepoint(-6.17706680000, 53.300791),4326), 900913))
LIMIT 1;
--Time: 248 ms

While both of these are quite fast, I'm unsure if this is the best way of performing this calculation. Both methods assume that the geo-coded point will be <50m away from a road. I also want to figure out the best way of finding two roads and cannot think of an alternative, other than to run to separate queries sequentially, which would mean that it takes ~0.5 seconds in total using the current set-up.


texts: How to do point-to-point near analysis conditional on common attributes of the two layers in Arcgis 10.2 (advanced)

texts:
I am measuring the closest straight distance from cities of interest (in a layer) to their capital cities (in another layer).  Basically, distance between a city and a capital city is calculated if only if their country codes are identical.  Arcgis license is advanced.

Solution:

The following solution is based on @dklassen.  Modification has been made so that I can loop a list of numerical number.

import arcpy
ctryID = [row.ID_0 for row in arcpy.SearchCursor(countryID_feature_class)]
arcpy.MakeFeatureLayer_management('PATH_TO_CITIES_FC','cityfc')
arcpy.MakeFeatureLayer_management('PATH_TO_CAPITALS_FC','capitalfc')

for grd in ctryID:
    arcpy.SelectLayerByAttribute_management('cityfc','NEW_SELECTION',""""ID_0" = {}""".format(grd))
    arcpy.SelectLayerByAttribute_management('capitalfc','NEW_SELECTION',""""ID_0" = {}""".format(grd))
    arcpy.Near_analysis('cityfc','capitalfc',"#","NO_LOCATION","NO_ANGLE")
    arcpy.CopyFeatures_management("cityfc", "c:/temp/ctry{}".format(grd))

arcpy.CreateFeatureclass_management("c:/temp/", "ctryall.shp", "POINT", "c:/temp/ctry1.shp", "DISABLED", "DISABLED")

for grd in ctryID:
    arcpy.Append_management("c:/temp/ctry{}.shp".format(grd), "c:/temp/ctryall.shp")
    arcpy.Delete_management("c:/temp/ctry{}.shp".format(grd))


texts: Creating Toolbox for ArcGIS that calculates route for maximum of 6 given addresses?

texts:
I like to build a Toolbox for the ArcGIS Network Analyst that takes a network and given addresses of a shapefile to calculate the best route.

Right now I am using the network-data of the [ArcGIS Network Analyst extension tutorial][1].

I have a shapefile with 1000 addresses, but I want to limit one route at one calculation with a maximum of 6 addresses.

Do I have to extract the addresses from the shapefile first before using the "Solve"-Tool or is there a better way?

In the end, I would like create an ArcGIS-Toolbox to have the possibility to choose the addresses with a drop-down menu for various address-lists and networks.

Right now I have the output of the Model Builder, Stores are my Addresses:

# -*- coding: utf-8 -*-
# ---------------------------------------------------------------------------
# Python.py
# Created on: 2017-02-16 10:11:54.00000
#   (generated by ArcGIS/ModelBuilder)
# Usage: Python <Routes>
# Description:
# ---------------------------------------------------------------------------

# Import arcpy module
import arcpy

# Script arguments
Routes = arcpy.GetParameterAsText(0)
if Routes == '#' or not Routes:
    Routes = "Route1\\Routes" # provide a default value if unspecified

# Local variables:
Streets_ND = "F:\\Streets_ND"
Route1 = "Route1"
Route1__3_ = Route1
Pizzaladen = "Pizzaladen"
Route1__4_ = Route1__3_
Stores = "F:\\Stores"
Routes__3_ = Stores
Route__3_ = Route1__4_
Solve_Succeeded = "true"
output5_shp = "F:\\output5.shp"

# Process: Make Route Layer
arcpy.MakeRouteLayer_na(Streets_ND, "Route1", "TravelTime", "FIND_BEST_ORDER", "PRESERVE_BOTH", "USE_TIMEWINDOWS", "", "ALLOW_UTURNS", "Oneway;RestrictedTurns", "USE_HIERARCHY", "", "TRUE_LINES_WITH_MEASURES", "")

# Process: Add Locations (2)
arcpy.AddLocations_na(Route1, "Stops", Pizzaladen, "", "5000 Meters", "", "Streets SHAPE;Streets_ND_Junctions NONE", "MATCH_TO_CLOSEST", "APPEND", "NO_SNAP", "5 Meters", "INCLUDE", "Streets #;Streets_ND_Junctions #")

# Process: Add Locations
arcpy.AddLocations_na(Route1__3_, "Stops", Stores, "Name Name #;Attr_Minutes ServiceTime 0", "5000 Meters", "", "Streets SHAPE;Streets_ND_Junctions NONE", "MATCH_TO_CLOSEST", "APPEND", "NO_SNAP", "5 Meters", "INCLUDE", "Streets #;Streets_ND_Junctions #")

# Process: Solve
arcpy.Solve_na(Route1__4_, "SKIP", "TERMINATE", "")

# Process: Select Data
arcpy.SelectData_management(Route__3_, "Routes")

# Process: Copy Features
arcpy.CopyFeatures_management(Routes, output5_shp, "", "0", "0", "0")

# Process: Select Data (2)
arcpy.SelectData_management(Stores, "Routes")

Is it possible to make a selection of 6 Stores/addresses among these Tools?


texts: Display maps from QGIS in my GUI created using PyQt5

texts: I wish to create a GUI using PyQt5 for an application that displays a map created using QGIS. How do I integrate QGIS and PyQt so that the map I wish to display can be shown in my app?


texts:
Leaflet has a 'click' event:

map.on('click', function(e) {
  alert(e.latlng);
});

This does not appear to work when map layers are displayed on the area clicked. Is there a way to use map.on('click') so that it ignores any layers on the map?


texts: Get approximately n points inside a polygon sf R

texts:
I need to get approximately n points in a pattern as regular as possible inside a polygon.

Current code:

library(tidyverse)
library(sf)
library(raster)
library(ggplot2)

n_points_wanted = 200

get_n_plants_grid <- function(shp, n_points) {

  area_box <- st_bbox(shp) %>% st_as_sfc() %>% st_area()
  area_shp <- shp %>% st_area() %>% sum()

  n_grid <- ((area_box / (area_box - area_shp  ) )*n_points) %>% sqrt()

  n_grid <- n_grid %>% as.numeric() %>% round()

  return (  n_grid )

}

shp <- getData('GADM', country = 'aut', level = 0) %>% st_as_sf()

n_grid <- get_n_plants_grid(shp, n_points_wanted)

grid <- st_make_grid(shp, n=c(n_grid, n_grid),
                     what = "centers",
                     square=TRUE) %>% st_intersection(shp)

ggplot() +
  geom_sf(data = shp) +
  geom_sf(data = grid, col ="red")

print(length(grid))

In this example I got 393 points inside the polygon but wanted 200.

Code corrected after @Spacedman’s advice.


texts: Use ogr2ogr to delete content of table and load new data from shapefile

texts:
Is it possible to load a shapefile to a PostGIS database by replacing the content of the table with the new content of the shapefile WITHOUT DROPPING the table?

I want to do this using ogr2ogr.

In the past I have used shp2pgsql, which gives you the option of four flags (-drop|-append|-create|-prepare) but none of these options actually maintains the table and deletes the data.

Is there something like this in ogr2ogr?


texts: Surface Volume: output table given all in one column?

texts:
The model iterates through a Point Feature Class to calculate different volumes. As an output I get two different styles of tables, both in one and the same iteration. One is the actual correct output table, the other consists of one column only, having all the values written one after the other separated by commas (see screenshots below). With the second table, there is no further access to the listed values possible.

It seems like a random thing to me when to get the actual table or the error one. Does anyone know about that problem?

The model is used to get the possible volume of retention basins. It is iterating through intersection points of a transport input (line) and a water body input (line), then creating a raster extraction that lies on the upper stream part of the buffer polygon only. With the z_values of the intersection as plane hight the surface volume below that level is calculated for every single intersection.

A "collect values" and "merge" model follows as the next step, to have all output tables in one.


texts: Where can get different satellite images for the same scene?

texts:
In a research paper, the authors make an experiment using 4 different images which are SPOT band 4, PalSAR HH polarization, Landsat TM band 5, DEM data. I attach one figure as example http://prntscr.com/fv6j9t.

I want to get images from a different source for the same scene. where can I find this for free?


texts: How to create a map from separate ascii files without merging them?

texts: I have three ascii files including longitudes, latitudes and pixel values. I just wonder whether there is anyway to display a map from the three files and save it in geotiff without merging them into a single file. I am using QGIS and ENVI-Idl.


texts: Can I upload a video file to embed in a custom infowindow in Cartodb?

texts: I want to embed video clips into custom infowindows in Cartodb, and I know how to do this via a youtube link, but is there a way I can embed a video from a file I have saved on my computer? Or does the video need to come from a website?


texts: ODP.NET invoked member is not supported in a dynamic assembly

texts:
I am using ODP.NET to get access to a geometry column in an Oracle database.
My code (derived from Oracle's sample) often (but not always) crashes on the last line where the geometry column value is converted to the SdoGeometry type.

Any ideas?

OracleConnection oraConnection = new OracleConnection("User Id=scott;Password=tiger;Data Source=oracle");
oraConnection.Open();
OracleCommand oraCmd = new OracleCommand();
oraCmd.Connection = oraConnection;
oraCmd.CommandText = "select geometrie from gw_pro where pro_id=3357";
OracleDataReader oraReader = oraCmd.ExecuteReader();
oraReader.Read();
Geo.SdoGeometry geom = (Geo.SdoGeometry)oraReader.GetValue(0);


texts: Getting attribute text error after running tool in QGIS?

texts:
I'm Korean, and I used Korean language text in attribute table.

I ran Add raster values to points. (input shapefile and raster(aspect))

The tool works well, but result file's attribute table has an error. Every Korean text character is changed to a question mark (?).

How can I maintain Korean text in an attribute table after running the tool add raster values to points?


texts:
I cannot open dashboard. I check ports and found that suite geoserver is working. (I click several times on dashboard.exe)

Please, tell me how to open dasboard.exe

UPD: when i enter ...localhost:8080/geoserver/web/ in browser i go to the geoserver. It is works.


texts: CartoCSS: Don't show border between adjacent polygons

texts:
In OSM, riverbeds are created as multiple adjacent or even overlapping polygons.

When styling these polygons in CartoCSS with a border, the overlapping can be seen:

I'd like to hide that line through the river. I guess that can be done either by merging the polygons, or by hiding the line between adjacent polygons.

Is that possible with CartoCSS / TileMill?


texts: I don't understand how v.voronoi GRASS module works with flag -a

texts:
I'm trying to use v.voronoi with flag -a. An example in the manual shows Voronoi diagram for urban areas but each try returns:

v.voronoi -a --overwrite --verbose input=FdT_fromSIE_int_parcelle_culturale_L93_nettoyee@PERMANENT output=test_voonoi_flaga

Reading features... Maximum segment length: 3.26421
ERROR: Found 0 vertices in <FdT_fromSIE_int_parcelle_culturale_L93_nettoyee@PERMANENT>, but at least 2 are needed

My input shapefile is type polygon.


texts: Snapping points to lines in ArcGIS Desktop and automate using VBA?

texts:
I use ArcGIS Desktop and I want create a network with a point feature and a line feature. I need my point features to be on the line features so I should do Snap point to line.

However, I do not know how do it in ArcGIS, and then automate it with code for this in VBA.


texts: Loading GEOJSON in Leaflet and separating attributes

texts:
I'm making a cloropleth map with geoJSON and Leaflet. Right now each year has a separate geoJSON file and I load each of them as below:

function loadGeoJSON(data) {
        var json = null;
        $.ajax({
            async: false,
            global: false,
            url: data,
            dataType: "json",
            success: function (data) {
                json = data;
            }
        });
        return json;
    }

  var ccc2 = loadGeoJSON("t2015.geojson");
    var t2015 = L.geoJson(ccc2, {
        style: style,
  onEachFeature: onEachFeature

    });
map.addLayer(t2015);

var ccc3 = loadGeoJSON("t2014.geojson");
    var t2014 = L.geoJson(ccc3, {
        style: style,
  onEachFeature: onEachFeature

    });

The geoJSON are of countries and are the same each year. I would like to combine the geoJson files into one with each feature having the attributes: 2014_1, 2014_2, 2014_3, 2015_1, 2015_2, 2015_3 etc. Is it possible to load one file with all years and break it into variables which only has the attributes from a single year (i.e. 2014_1, 2014_2, 2014_3) so they can be separate layers on the map?

EDIT: SAMPLE COMBINED GEOJSON (note: geographies are not in full)

{
"type": "FeatureCollection",
"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } },

"features": [
{ "type": "Feature", "properties": { "FIPS_CNTRY": "AG", "CNTRY_NAME": "Algeria", "2013": 9698, "2013_2": 9545, "2014": 266, "2014_2": 52918, "2015": 91, "2015_2": 842 }, "geometry": { "type": "Polygon", "coordinates": [ [ -5.152134895324707, 30.180469512939453 ], [ -5.139167308807373, 30.192359924316406 ], [ -5.080972671508789, 30.262359619140625 ], [ -5.048055648803711, 30.316665649414063 ], [ -5.029167175292969, 30.359165191650391 ], [ -5.010833740234375, 30.393886566162109 ]]}},
{ "type": "Feature", "properties": { "FIPS_CNTRY": "AL", "CNTRY_NAME": "Albania", "2013": 7326, "2013_2": 1245, "2014": 834, "2014_2": 135, "2015": 785, "2015_2": 134 }, "geometry": { "type": "Polygon", "coordinates":  [ [ 20.791923522949219, 40.431541442871094 ], [ 20.787220001220703, 40.394721984863281 ], [ 20.758609771728516, 40.311943054199219 ], [ 20.73680305480957, 40.307220458984375 ], [ 20.714302062988281, 40.270099639892578 ], [ 20.721246719360352, 40.225135803222656 ]] }},
{ "type": "Feature", "properties": { "FIPS_CNTRY": "AN", "CNTRY_NAME": "Andorra", "2013": 6698, "2013_2": 9578, "2014": 24266, "2014_2": 546, "2015": 425, "2015_2": 3145 }, "geometry": { "type": "Polygon", "coordinates":  [ [ 1.445833206176758, 42.601943969726563 ], [ 1.486527681350708, 42.650413513183594 ], [ 1.559722185134888, 42.655967712402344 ], [ 1.698333263397217, 42.626106262207031 ], [ 1.738610982894898, 42.616386413574219 ], [ 1.78171968460083, 42.569961547851563 ]]}},
{ "type": "Feature", "properties": { "FIPS_CNTRY": "AU", "CNTRY_NAME": "Austria", "2013": 2592, "2013_2": 311, "2014": 66, "2014_2": 511, "2015": 242, "2015_2": 8664 }, "geometry": { "type": "Polygon", "coordinates":  [ [ 10.471235275268555, 46.871353149414063 ], [ 10.488205909729004, 46.935993194580078 ], [ 10.429998397827148, 46.984161376953125 ], [ 10.390763282775879, 47.002567291259766 ], [ 10.350624084472656, 46.991241455078125 ], [ 10.323331832885742, 46.955551147460938 ]] }},
]
}


texts: bad results NDVI in rasterio with jp2 sentinel-2 bands

texts:
I am using 2 Sentinel-2A bands to compute the NDVI. Here is my code:

import os, rasterio

outfile = r'some\path\ndvi.tif'
#url to the bands
b4 = 'http://sentinel-s2-l1c.s3.amazonaws.com/tiles/30/T/TK/2017/4/12/0/B04.jp2'
b8 = 'http://sentinel-s2-l1c.s3.amazonaws.com/tiles/30/T/TK/2017/4/12/0/B08.jp2'

#open the bands (I can't believe how easy is this with rasterio!)
with rasterio.open(b4) as red:
    RED = red.read()
with rasterio.open(b8) as nir:
    NIR = nir.read()

#compute the ndvi
ndvi = (NIR-RED)/(NIR+RED)
#print(ndvi.min(), ndvi.max()) The problem is alredy here

profile = red.meta
profile.update(driver='GTiff')
profile.update(dtype=rasterio.float32)

with rasterio.open(outfile, 'w', **profile) as dst:
    dst.write(ndvi.astype(rasterio.float32))

The problem is that, this is working fine, except for negative values. It seems that the num (NIR-RED) gets positives values no matter what. So vegetation and all the ndvi values between 0 and 1 keeps good values, but ndvi values between 0 and -1 gets weird values.
I can solve this saving RED and NIR to 2 GTiff rasters in local, but I guess that there must be one solution in order to avoid to do that.


texts: What is the difference between ST_Union and ST_UnaryUnion?

texts:
Reading the PostGIS manual, it appears both St_Union and ST_UnaryUnion should dissolve boundaries of overlapping polygons. However, a simple test shows the following:

Start data:

insert into test (name, geom) values
('Top Left', st_geomfromtext('POLYGON((0 1, 0 3, 2 3, 2 1, 0 1))', 4326)),
('Bottom Right', st_geomfromtext('POLYGON((1 0, 1 2, 3 2, 3 0, 1 0))', 4326));

ST_Union:

create table test_union as select
st_union(geom) as geom
from test;

ST_UnaryUnion

create table test_unary as select
st_unaryunion(geom) as geom
from test;

What explains this confusing behavior?


texts: How to get output value from Feature Compare (schema) and append to existing csv?

texts:
I am working on a data integration project and am trying to single out feature classes which have a schema that varies from two base feature classes (one for point, one for polyline features).

I've played around with the Feature Compare tool to compare schemas, but find it tedious to use this tool, even in batch mode, because it outputs a single text file with a few values for each comparison performed.

I started writing a script in Python to loop through a list of features in a list and compare their schema based on geometry type. However, the problem of either integrating the output into a single text file persists. Is there a way to obtain this value and append it to an existing text file?

(Messy) code so far as follows:

# Import arcpy module
import arcpy

# Import environments; set workspace
from arcpy import env

# Set workspace
env.workspace = "Data_Comparison.gdb"

# List features to be compared
fclist = arcpy.ListFeatureClasses()
print fclist

# Loop through list to separate point and line feature classes
for fc in fclist:
    desc = arcpy.Describe(fc)
    name = desc.file
    geo = desc.shapeType

    print ("Feature name: " + name)
    print ("Feature geometry: " + geo)

    if geo == "Point":
        arcpy.FeatureCompare_management("swr_net_nodal", fc, "ENTID", "SCHEMA_ONLY", "", "", "", "", "", "", "CONTINUE_COMPARE", "fc + '_compare'")
        output = arcpy.GetMessages()
        print output


texts:
QGIS 2.8.2 on windows 8.1 produces the below error every time it is closed.

Anyone have any ideas as to why or possible fixes ?


texts: Automating classification of symbology in PyQGIS

texts:
I have created this code to categorize my layer symbology based on the attribute plant_allo. However I still have to press the classify button which I also want to automate in PyQGIS. Can anyone help me with this last step?

Here is my code:
year = 2030
allocated_regions = QgsVectorLayer(output_folder + "regions_w_plant_alloc_" + str(year) + ".shp", "allocated_regions_" + str(year), "ogr")
QgsProject.instance().addMapLayer(allocated_regions)

# Field name
expression = 'plant_allo'
# Set the categorized renderer
renderer = QgsCategorizedSymbolRenderer(expression)
allocated_regions.setRenderer(renderer)
allocated_regions.triggerRepaint()


texts: Append table in PostGIS - Geometry Violation

texts:
In PostGIS I have imported a shapefile - some of this is Invalid according to PostGIS.

Following online help I have first split the invalid polygons into a new table:

CREATE TABLE s_invalid AS
SELECT * FROM sdata
WHERE NOT ST_IsValid(the_geom);

DELETE FROM sdata
WHERE NOT ST_IsValid(the_geom);

I then buffered the invalid polygons by 0 to sort out the issue:

UPDATE s_invalid
SET the_geom = ST_Buffer(the_geom,0);

SELECT count(*), ST_IsValid(the_geom)
FROM s_invalid
GROUP BY ST_IsValid;

All seems ok, now I want to append the cleaned-up polygons back into the original table and this is where I have an issue. I have tried:

INSERT INTO sdata
SELECT * FROM s_invalid;

But there seems to be a mismatch with the geometry type from the best I can see - I get the message:

ERROR:  new row for relation "sdata" violates check constraint "enforce_geotype_the_geom"

********** Error **********

ERROR: new row for relation "sdata" violates check constraint "enforce_geotype_the_geom"
SQL state: 23514

How can I sort this... where is the issue? Is there an easier way to sort this out?
Thanks, I really want to get to grips with this!


texts: Request long/lat precision and request too slow

texts:
I've a PostGIS DB of 18Gb. The tab node (id, tags, geom) is very huge and I need to request long/lat on this table but it's very long when I'm using this sql request :

SELECT * FROM nodes WHERE geom = St_SetSRID(ST_MakePoint(longitude,latitude),4326)

like 300s :/ How can I improve that (i need to do more than 100). And i have another problem with the precision of long/lat because they are note the same (mine and in the Db).

SOLUTION : i just need to use the right request !!!

SELECT id
  FROM nodes ORDER BY geom <-> St_SetSRID(ST_MakePoint(lng,lat),4326)
    asc limit 10;

for 52M nodes it's take MAX(200ms) !!!


texts: Generate essential grid covering polygon with OGR & Python

texts:
I'm trying to come up with a method to recreate the type of grid shown below.
The grid is using a custom projection measured in meters. each cell is 256 x 256 meters.
Ideally I'd like to easily use the same method for a grid with cells of 512 x 512 meters etc.

It must be a grid with cells that are either inside or intersect with the polygon.. no cells that are completely outside of the polygon.
My thought process so far as been to get the extents of the polygon (shown above), fill the extents with the desired grid spacing (not implemented yet), then use ogr intersect functions to export the intersections to a new shapefile... but there must be a more elegant way to do this..
The code that I've started with is shown below which so far only creates a shapefile based on extents.
#!/usr/bin/env python2
import sys
from osgeo import ogr
import os

OutBaseName = os.path.splitext(sys.argv[1])[0]
resolution = "test"
# Get a Layer's Extent
inDriver = ogr.GetDriverByName("ESRI Shapefile")
inDataSource = inDriver.Open(sys.argv[1], 0)
inLayer = inDataSource.GetLayer()
extent = inLayer.GetExtent()
print extent

# Create a Polygon from the extent tuple
ring = ogr.Geometry(ogr.wkbLinearRing)
ring.AddPoint(extent[0],extent[2])
ring.AddPoint(extent[1], extent[2])
ring.AddPoint(extent[1], extent[3])
ring.AddPoint(extent[0], extent[3])
ring.AddPoint(extent[0],extent[2])
poly = ogr.Geometry(ogr.wkbPolygon)
poly.AddGeometry(ring)

# Save extent to a new Shapefile
outShapefile = OutBaseName + "_LMM_Blocks_" + resolution  + "m.shp"
outDriver = ogr.GetDriverByName("ESRI Shapefile")

# Remove output shapefile if it already exists
if os.path.exists(outShapefile):
    outDriver.DeleteDataSource(outShapefile)

# Create the output shapefile
outDataSource = outDriver.CreateDataSource(outShapefile)
outLayer = outDataSource.CreateLayer("states_extent", geom_type=ogr.wkbPolygon)

# Add an ID field
idField = ogr.FieldDefn("id", ogr.OFTInteger)
outLayer.CreateField(idField)

# Create the feature and set values
featureDefn = outLayer.GetLayerDefn()
feature = ogr.Feature(featureDefn)
feature.SetGeometry(poly)
feature.SetField("id", 1)
outLayer.CreateFeature(feature)
feature = None

# Save and close DataSource
inDataSource = None
outDataSource = None


texts: How to add grid labels on the border of the data frame?

texts:
Licence: Advanced
  Version: ArcMap 10.1

I'm trying to place the labels of my measured grid at the border of my data frame. I've tried to edit the labels in the "Data Frame Properties -> Grids -> Properties" but nothing works. Only the border of the polygon will be labeled.

How can i fix that?


texts: Calculate shortest paths for each output area in catchment through road network to a single point (i.e. store)

texts:
The title sums it up pretty much. I'm looking to derive trip distributions to a single store using population data and distance within a prescribed distance catchment. I have calculated the catchment area (polygon) using v.net.iso in GRASS and have clipped weighted population centroids that fall with the catchment accordingly. I therefore know the total population within the catchment, alongside the population of each OA centroid.

The ideal output would be:
OA Geocode (column 1); Distance to Store (column 2); OA population (column 3)

I could achieve the above using the 'distance to nearest hub'geo-algorythm in QGIS, however, it uses straight line distance which is not ideal.

I have tried to perform the following in GRASS using the v.distance function, however, i am unable to output results by OA (node), only by the road network (arc), which is not what i need.

As well as the above, my main goal is to be able to find a way of identifying whether certain paths utilise certain roads to access the store. This will give me my trip distribution for my study network (i.e. to model individual junctions).

The following would be ideal:

OA Geo Code (column 1); Distance to Store; Road Link ID (specific road link) 0 or 1 (1 if used); Road link ID ""; Road Link ID ""; Road Link ID ""...; OA Population.

The above would be a very simplistic gravity model. I can do those separately, the identification of link would be by hand a though and would be very laborious.


texts: Generate global extent hexagon grid tessellation in ArcGIS Pro

texts: I am trying to generate a global hexagon grid at 1 km2 resolution spanning the entire globe as an equal area grid using ArcGIS Pro's 'Generate Tessellation' function. I've tried it in different CRSs using a variety of global layers (polygons and points) to set the extend ('as specified below') but I keep getting an out-of-bounds error or the tessellation generated is empty. How can this best be achieved?


texts:
I work with arcview 10.3, and i try to extent to xy coordinate with arcpy?

I am using this code:

import arcpy, os, sys
from arcpy import env

env.workspace = r"C:\Project"
pnt = '178434 ,662446'
for mxdname in arcpy.ListFiles("*.mxd"):
    print mxdname
    mxd = arcpy.mapping.MapDocument(r"C:\Project\\" + mxdname)
    df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
    df.extent = pnt
    print 'extent'
    mxd.save()
del mxd

but i get en error:

  ValueError: Invalid string value for Extent argument


texts: How can I fetch related db records in a GeoExt webapp?

texts:
I am working on a web app which uses GeoExt, Geoserver, and PostGIS.

I would like to be able to identify features in this applications using the wms getFeatureInfo and then access a features related records (like photos, documents etc.).

Does anyone know how I could achieve this directly with GeoExt or would I have to venture into PHP or some kind of webservice?

edit - relationship between tables

Thanks

Ando


texts: Enabling Multiple Layers in same Infowindow of CartoDB.js?

texts:
I've a website where using Cartodb.js I'm adding two layers using

cartodb.createlayer(map, 'url').addTo(map);

Info window is on for both layers. While the features of one layer overlap with features of other layer. So, when I click the overlapping area, two info windows open. While I want to open one info window including data from both windows. Is there a way to customize info windows in my case.


texts: Issue with CRS in R's writeRaster and SAGA's raster-format

texts:
I have a normal RasterLayer in R and I would like to write it to a file and visualize it in QGIS. When I look at in in R it says:
class      : RasterLayer
dimensions : 3035, 1469, 4458415  (nrow, ncol, ncell)
resolution : 12.5, 12.5  (x, y)
extent     : 386042.8, 404405.3, 6293348, 6331286  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=19 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
source     : memory
names      : layer
values     : 2.220446e-16, 0.9845922  (min, max)

The I do a writeRaster(<rasterlayer>, <file.tif>, format = "GTiff", overwrite=TRUE). When I check the coordinate reference with gdalsrsinfo on the command line I get this information:
Driver: GTiff/GeoTIFF
Files: no_filter.tif
Size is 1469, 3035
Coordinate System is:
BOUNDCRS[
    SOURCECRS[
        PROJCRS["WGS 84 / UTM zone 19S",
            BASEGEOGCRS["WGS 84",
                DATUM["World Geodetic System 1984",
                    ELLIPSOID["WGS 84",6378137,298.257223563,
                        LENGTHUNIT["metre",1]]],
                PRIMEM["Greenwich",0,
                    ANGLEUNIT["degree",0.0174532925199433]],
                ID["EPSG",4326]],
            CONVERSION["UTM zone 19S",

However, when I drag and drop it to QGIS I get the Unknown CRS output.
Is there anything special that should be considererd when writing raster-files or something like a "best practice"?
UPDATE
SO when I just print the raster with
raster_objects[[1]] for example I get:
> raster_objects[[1]]
class      : RasterLayer
dimensions : 3035, 1469, 4458415  (nrow, ncol, ncell)
resolution : 12.5, 12.5  (x, y)
extent     : 386042.8, 404405.3, 6293348, 6331286  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=19 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
source     : memory
names      : layer
values     : 2.220446e-16, 0.9845922  (min, max)

When I do a rasterWrite(raster_objects[[1]], "test/test.sdat", format="SAGA"), I get:
class      : RasterLayer
dimensions : 3035, 1469, 4458415  (nrow, ncol, ncell)
resolution : 12.5, 12.5  (x, y)
extent     : 386042.8, 404405.3, 6293348, 6331286  (xmin, xmax, ymin, ymax)
crs        : NA
source     : /home/robin/geodata/hiwijason/workPkg_LslPrd_DEMsmoothing/maipo_basin/data/saga/raster_predictions/test.sgrd
names      : layer
values     : 2.220446e-16, 0.9845922  (min, max)

And I have no clue why the crs-line is saying NA in the writing-operation
UPDATE
Maybe this has to do with the raster beeing in SAGA-Format (.sdat)?! Because when I drag and drop the .sdat-file into qgis it gets georeferenced directly. However, when I import it into R and write it out again with:
r = raster("...sdat")
writeRaster(r, "...tif", format="GTiff")

I lose the crs somewhere on the way. Is there anything I'm doing completely wrong? Would readGdal maybe offer help?


texts: Problem with saving raster using writeRaster in R

texts:
I'm running some analysis in RStudio using raster and rgdal package.
Where I assign variable to some raster, or stack for example:

analysis = raster("M-34-91-B-a-1-4-4.img")

Then I run some analysis on this variable (classification, some function for every pixel), and I'm trying to save it using writeRaster function:

writeRaster(analysis, "analysis.tif")

My tif file is really small (53 KB) and when I'm trying to open it all I get is blank square with some black dots all around) <--- image

However, when I use plot function

plot(analysis, main = "Max height")

Everything in R works fine, and I'm able to export this plot using RStudio interface and resulting image is all fine (and weights more ~ 1 173 KB) (image below)

My question is whether is something wrong with my computer, or my RStudio. I feel really powerless at this point, tried reinstalling both R and RStudio, installing different file viewers like IrfanView etc. When I try to open written raster in ArcPro it works rather fine (image below)

I just noticed that black pixels from first image, are probably no data pixels after adding it to ArcPro, is it possible then that no data values (probably -9999) messed up the stretch in writeRaster function? I'll be thankfull for any help attemps, thanks


texts:
In old turfjs versions, it was possible to intersect a line and a polygon, like in this example.

But changing the turf.js version to the last one makes it fail, as shown in this other example.

The docs say that the intersect function works only with two polygons, so the error is documented.

Is it possible to use some other function to do that? Is there any example around there showing how to do it?