texts: Cloud mask in Surface Reflectance Landsat 8 test

texts:
I've noticed that my cloud mask wasn't working, so I've tried this simple test:
https://code.earthengine.google.com/50699c2eaa1a873ccd28f26c583c5a45

But my data uses the Surface Reflectance Landsat 8 imagery so, I just changed to that, and since there's no fmask band in this, I've changed to the 'pixel_qa' that also make those distinctions. I thought it was the same, but isn't working.

Code:

//Choose country using GEE Feature Collection

var region = ee.FeatureCollection('ft:1tdSwUL7MVpOauSgRzqVTOwdfy17KDbw-1d9omPw').filterMetadata('Country', 'equals', 'Portugal');

//Add region outline to layer ‐ for selected countries

Map.addLayer(region,{}, 'Portugal');

var landsat8 = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')

// Fmask classification values
var FMASK_CLEAR_GROUND = 0;
var FMASK_WATER = 2;
var FMASK_CLOUD_SHADOW = 3;
var FMASK_SNOW = 4;
var FMASK_CLOUD = 5;

var mosaic = landsat8
  .filterBounds(region)
  .filterDate('2017-08-01', '2017-10-11')
  .mosaic();

// Update the mask on our mosaic to mask cloud and cloud shadow pixels
var fmask = mosaic.select('pixel_qa');
var cloudMask = fmask.neq(FMASK_CLOUD).and(fmask.neq(FMASK_CLOUD_SHADOW));
var maskedMosaic = mosaic.updateMask(cloudMask);

Map.addLayer(fmask, {min:0, max:5, palette:'green, blue, black, cyan, pink, white'}, 'Fmask');
Map.addLayer(maskedMosaic.select('B4'), {min:0, max:0.5, palette:'yellow, green'}, 'Masked NIR');


texts: Converting real number fields to a simple integer field

texts: I have a shapefile with attribute table that has number fields that have five zeros after the decimal point, having been imported that way. I don't need decimals and want the figures as simple intigers. How can I convert them in QGIS?


texts:
I want to calculate the counts of each land use type based on the unique Landuse code from each raster dataset.

Raster dataset consists of six major land use and their unique codes are below:

Agriculture: Landuse code (2100, 2200, 2300, 2400)
Barren Land: Landuse code (7100, 7200, 7300, 7400, 7500, 7600)
Forest: Landuse code (4110, 4120, 4210,4220, 4230, 4311,
4312, 4321, 4322, 4410)
Urban: Landuse
code (1110, 1120, 1130, 1140, 1150, 1200, 1211, 1214, 1300, 1400, 1410, 1411, 1420)
Water: Landuse code (1419, 5100, 5190)
Wetlands: Landuse code (1461, 1711, 1750, 1850)

I have the code only for agricultural area to simply test whether the code is working or not for one category of land use.

When I ran the code neither did I get any error message nor the result for agricultural land. I know it's a pretty big description, sorry for that.

I am providing the code below:

 import arcpy, os
 from arcpy import env
 #To overwrite output
 arcpy.env.overwriteOutput = True
 #Set environment settings
 env.workspace = "C:/Subhasis/Project-06-02-14/Landuse/New folder/Merged/Landuse_raster"
 outws="C:/Subhasis/Project-06-02-14/Landuse/New folder/Merged/Landuse_raster/Table"
 inraster = arcpy.ListRasters("*", "GRID")

 for i in inraster:
     flds = ("VALUE", "COUNT")
     dct = {row[0]:row[1] for row in arcpy.da.SearchCursor(i, flds)}
     Agriculture=[]
     Barrenland=[]
     Forest=[]
     Urban=[]
     Water=[]
     Wetland=[]
     for j in dct:
         if (dct.keys()==2100) or (dct.keys()==2200) or (dct.keys()==2300) or (dct.keys()==2400):
             Agriculture.append(dct.values())
             print Agriculture
         else:
             break


texts: Parallel processing in QGIS3 on Windows system

texts:
I am working on a big vector data set with more than 2.4 million features containing 10x10m squares. I wrote a small script in Python to delete uninteresting fields and perform calculations on the remaining fields of interest, but this process comes to no end.
I am wondering how I can speed up processing.
I thought about the implementation of parallel processing using all four cores, but I didn't find any helpful information.
What I want to do is split the whole area in four parts, than for every tile I want remove to all fields where I'm not interested in and make some calculation on the remaining fields followed by a conversion to a raster. This is my code:
from qgis.core import QgsTrackedVectorLayerTools

layer = iface.addVectorLayer("[mypath]name.gpgk]","name","ogr")

##################################
###Filter Fields by selecting columns of interest

idx_keep = []
idx_all = []
idx_delete = []

fields = ["basalareap","basalareas","basalaread"]

for i in fields:
    idx_keep.append(layer.fields().indexFromName(i))

field_names = [field.name() for field in layer.fields()]

for i in field_names:
    idx_all.append(layer.fields().indexFromName(i))

idx_delete = [x for x in idx_all if x not in idx_keep]

layer.startEditing()

layer.deleteAttributes(idx_delete)

##################################
##Create new fields

pv = layer.dataProvider()
pv.addAttributes([QgsField('Share_SW',QVariant.Double), \
QgsField('Share_HW', QVariant.Double)])

layer.updateFields()

expression_SW = QgsExpression('("basalareap"+"basalareas")/("basalareap"+"basalareas"+"basalaread")')
expression_HW = QgsExpression('("basalaread")/("basalareap"+"basalareas"+"basalaread")')

context = QgsExpressionContext()
context.appendScopes(QgsExpressionContextUtils.globalProjectLayerScopes(layer))

#with edit(layer):
layer.startEditing()

for f in layer.getFeatures():
    context.setFeature(f)
    f['Share_SW'] = expression_SW.evaluate(context)
    f['Share_HW'] = expression_HW.evaluate(context)
    layer.updateFeature(f)

QgsVectorFileWriter.writeAsVectorFormat(layer,'[mypath]','utf-8',layer.crs(),'ESRI Shapefile')
iface.mainWindow().findChild(QAction, 'mActionToggleEditing').trigger()
#Here convert to raster

Can anybody give me a hint how to realize this?
Perhaps there is a tutorial.
I'm using Windows 10 with QGIS 3.10.
Update
I think the geopandas approach from olivier is faster and cleaner. I asked an seperate question for the geopandas package:
Parallel processing in Pandas and already get a good answer.


texts: Minimizing turns in ArcGIS Network Analyst routing?

texts:
I have a network of pedestrian routes and I want to provide them with clean directions, but Network Analyst chooses the first route that it can turn towards destination and gives me a zig-zag road.

Is there a way to minimize the number of turns?

Is there a way to count turns?


texts: How do I select rows in an attribute table of a vector layer grid which are outside of another vector layer?

texts:
So what I am working on is accessibility analysis, I want to know how many percent of the population of a country have access to certain health services within a given driving time. For example is it possible to visit a certain specialist within 45 mins driving time? So what I need to know is how much of the population falls within my drive time isochrones and how much outside of it.

I have my isochrone layers for driving distance and I have a Vector layer population grid (from Eurostat) where each population value is tied to a 1*1km cell. (I don´t know if this is important, but the data came in as a reference grid and a csv file I had to manually join to the grid by Grid ID). I tried using the "difference" vector processing tool in QGIS to substract the isochrone layer from my population grid layer. What I got out looks exactly how I wanted it to look, it shows me just the population which falls outside of the isochrones. However, the attribute table stayed the same, it has exactly as many entries as before so it is useless for calculations. I want to compare the underserved population to the overall population.

Do you know how to perform a similar difference operation on my grid that will also remove the connected attributes to the removed grids?


texts: Using Symbology Scaling when printing to PDF using QGIS

texts:
I am using QGIS 3.4.

When I go to print my map as a PDF the orange diagonals I added, which show up nicely in print composer, are scaled to the point of invisibility.

How can I fix this so that they are not scaled down when exported?

In Print Composer:

In PDF:


texts: Renaming Columns in Geotools FeatureType

texts:
I am attempting to build a simple java program that loads shapefiles into a Postgis 2.0 database. Everything is working out great except that I would like to be able to modify the FeatureType in the following ways:

1) change the geometry column name from "the_geom" to "geom"
2) set all columns to be lower case.
3) change the primary key from serial to bigserial.

The following is my current code which basically is just transferring the shpSchema to the pgSchema. I am fairly certain that what I am trying to do is fairly simple and common, but I can find no obvious way to do it or any documentation --or even any examples of someone trying to do this.

    Map shpConnectParams = new HashMap();
    shpConnectParams.put("url", shapefile.toURI().toURL());

    DataStore shpStore = DataStoreFinder.getDataStore(shpConnectParams);

    String shpTypeName = shpStore.getTypeNames()[0];
    FeatureSource shpFS = shpStore.getFeatureSource(shpTypeName);

    SimpleFeatureType shpSchema = (SimpleFeatureType) shpFS.getSchema();
    CoordinateReferenceSystem crs = shpSchema.getCoordinateReferenceSystem();
    Map pgParams = new HashMap();
    out.println(shpSchema.getGeometryDescriptor().getType());
    for (AttributeDescriptor attrDesc : shpSchema.getAttributeDescriptors()) {

        out.println(attrDesc.getName());
    }

    pgParams.put("dbtype", "postgis");                   //must be postgis
    pgParams.put("host", "localhost");      //the name or ip address of the machine running PostGIS
    pgParams.put("port", new Integer(5432));             //the port that PostGIS is running on (generally 5432)
    pgParams.put("database", "test");              //the name of the database to connect to.
    pgParams.put("user", "postgres");                    //the user to connect with
    pgParams.put("passwd", "test!");                 //the password of the user.
    pgParams.put("schema", "scratch");   //the schema of the database
    pgParams.put("create spatial index", Boolean.TRUE);

    SimpleFeatureTypeBuilder builder = new SimpleFeatureTypeBuilder();

    builder.addAll(shpSchema.getAttributeDescriptors());
    builder.setCRS(crs);
    builder.setName(shpTypeName);

    DataStore pgDatastore = DataStoreFinder.getDataStore(pgParams);
    pgDatastore.createSchema(builder.buildFeatureType());

    SimpleFeatureSource pgFS = pgDatastore.getFeatureSource(shpTypeName);

    SimpleFeatureType pgSchema = pgFS.getSchema();

    Transaction transaction = new DefaultTransaction("create");

    if (pgFS instanceof SimpleFeatureStore) {
        SimpleFeatureStore pgFeatureStore = (SimpleFeatureStore) pgFS;

        SimpleFeatureCollection features = (SimpleFeatureCollection) shpFS.getFeatures();
        pgFeatureStore.setTransaction(transaction);
        try {
            pgFeatureStore.addFeatures(features);
            transaction.commit();

        } catch (Exception problem) {
            problem.printStackTrace();
            transaction.rollback();

        } finally {
            transaction.close();
        }
        System.exit(0); // success!
    } else {
        System.out.println(shpTypeName + " does not support read/write access");
        System.exit(1);
    }

    pgDatastore.dispose();


texts: GeoServer Web Admin Interface: 404 error

texts:
I'm working through the OpenGeo workshop and I'm stuck with setting up a GeoServer. When I try to access the Web Admin Interface localhost:8080/geoserver/, Tomcat reports a "HTTP 404 Status" error page. The same holds for the GeoWebCache configuration page.

Some additional info:

My system: Ubuntu 13.04 with OpenGeo Suite 3.0.2 installed as described in the official reference. Things work fine on windows by the way.
I started Tomcat 6.0.35 via $/etc/init.d/tomcat6 start as the installation details point out.
The log-file of tomcat at /var/log/tomcat6 reports a "severe" FileNotFoundException for the file masterpw.xml:

SEVERE: Exception sending context initialized event to listener instance of class org.geoserver.platform.GeoServerContextLoaderListener
org.springframework.beans.factory.BeanCreationException:
Error occured reading security configuration;
nested exception is java.io.FileNotFoundException:
/usr/share/opengeo-suite-data/geoserver_data/security/masterpw.xml (No such file or directory)

At the end of geoserver.log in /var/log/geoserver the following warning is reported:

WARN [support.DisposableBeanAdapter] - Invocation of destroy method failed on bean with name 'geoServerLoader':
org.springframework.beans.factory.BeanCreationNotAllowedException:
Error creating bean with name 'GWCLifeCycleHandler':
Singleton bean creation not allowed while the singletons of this factory are in destruction (Do not request a bean from a BeanFactory in a destroy method implementation!)

I'm new to web mapping and the above steps together with an internet recherche didn't resolve my problem. Do you have any suggestions?

Update: Setting up a GeoServer outside of the OpenGeo framework works just fine. The key is that there is an OS-independent binary provided on GeoServer.org using a webapp bundled in Jetty instead of Apache Tomcat. The steps it took were:

Download and unzip geoserver-2.3.2-bin.zip.
Move the folder to your preferred installation path and cd into it.
Add the environmental variable JAVA_HOME to your PATH:

export JAVA_HOME=/usr/lib/jvm/java-1.6.0-openjdk-amd64/jre

Note: This path is distribution dependent (here: Ubuntu 13.04)
Shut down the Apache Tomcat servlet engine:

sudo service tomcat6 stop

Start up the Server by executing the provided java archive in the folder:

java -jar start.jar

The Web Admin Interface is now available via web browser under

localhost:8080/geoserver

However, I suspect that if one knows about how to assign port numbers etc it might still work with Apache as installed by the OpenGeo Suite. Suggestions are welcome.


texts: Is there a python-module to convert data to a Shapefile (w/o GIS)?

texts:
Input data can be excel, text, whatever. I would like to produce shapefiles from data on a computer without GIS installed. Is there anything out there?

Fiona seems to be such a module.  Are there any other modules? Is Fiona the best choice?

I have data in a list like this:

mylist = [["lat1","long1","lat2","long2","var1","var2","var3"],
         [55.000,7.000,54.000,6.000,30,40,50]
         [54.000,6.000,52.000,8.000,30,40,50]
         [52.000,8.000,53.000,7.000,30,40,50]
         [53.000,7.000,50.000,6.000,30,40,50]
         [50.000,6.000,55.000,7.000,30,40,50]]

And I would like to make a shapefile (polyline) out of that.


texts: Raster Distance from each cell to nearest cell with specific value

texts:
To analyse the "virginity" of an area I created a binary raster with

1 = human infrastructure (road, housing, ...)
0 = no human infratructure (forest, agriculture, ...)

I would no like to calculate the distance from each cell (especially the 0-cells) to the nearest cell with value 1.

How? I use QGIS 2.18.


texts: arcpy Select_analysis where_clause objectid

texts:
Trying to copy features from sde to fgdb but only when the objectid is below a certain number.

outpath = r'C:\Users\CHOK\Desktop\SDE Conversion Project\TEST.gdb'
query = '"objectid" < ' + "'121'"
for featureclass in dataset:
        if featureclass == 'Contour_10m':
            output = outpath + os.sep + featureclass
            arcpy.Select_analysis(featureclass, output, query)

The code executes but only the schema of the sde is copied in the new feature class created in the gdb and no data with objectid<121 that I expected would be copied over is present.

Any ideas on what the issue is?


texts: Overlaying a polyshape to faceted maps in R with ggplot?

texts:
I have managed to create a faceted map from my data following the approach proposed here http://spatial.ly/2013/12/introduction-spatial-data-ggplot2/
with the following code:

ggplot(data = plot.data, aes(x = long, y = lat, fill = value, group = group)) +
  geom_polygon() + coord_equal() +
  facet_wrap(~variable,labeller = as_labeller(names), ncol=2)

Where

plot.data <- merge(bound.f, mdata.melt, by.x = "id", by.y = "uid")

and bound.f is the fortified spatial polyshape (with sub-national boundaries)
and mdata.melt is the data used in map, in the the long format.

And I get the following output:

So far so good. Now I would like to add (overlay) the boundaries of the African countries (which I have a shapefile, same projection than bound) to both maps.

In other words, I want to add the shape, which boundaries should be black and the shape with no fill.

Does anyone have a hint on how I could code this or even has a smarter approach to get beautiful maps into my publication?


texts: Reading S3 Sentinel-2 image files with rasterio

texts:
I am currently experiencing an 'Access Denied' error after using rasterio to open the Sentinel-2 AWS .jp2 image files (currently hosted by Sinergise and can be found here).

My code looks like this:

import rasterio

access_key = '*****************'
secret_access_key = '***********************'
region_name = 'eu-central-1'

Session = rasterio.env.Env(aws_access_key_id=access_key,
             aws_secret_access_key=secret_access_key,
             region_name=region_name)

url = 's3://sentinel-pds/tiles/10/S/DG/2015/12/7/0/B01.jp2'

with Session:
    with rasterio.open(url) as dataset:
            print dataset

The response is:

---------------------------------------------------------------------------
CPLE_AWSAccessDenied                      Traceback (most recent call last)
<ipython-input-1-4f6da7274649> in <module>()
     12
     13 with Session:
---> 14     with rasterio.open(url) as dataset:
     15         print dataset

/Users/Barney/anaconda/lib/python2.7/site-packages/rasterio/__init__.pyc in open(path, mode, driver, width, height, count, crs, transform, dtype, nodata,  **kwargs)
    191         raise ValueError(
    192             "mode string must be one of 'r', 'r+', or 'w', not %s" % mode)
--> 193     s.start()
    194     return s
    195

rasterio/_base.pyx in rasterio._base.DatasetReader.start (rasterio/_base.c:2847)()

rasterio/_base.pyx in rasterio._base.DatasetReader.start (rasterio/_base.c:2799)()

rasterio/_err.pyx in rasterio._err.CPLErrors.check (rasterio/_err.c:1773)()

CPLE_AWSAccessDenied: Access Denied

I'm unsure if I am using rasterio aws environments correctly, is the aws session token key required even if the secret and access keys are provided? I should also add that this function works correctly if provided with a Landsat-8 url (i.e. http://landsat-pds.s3.amazonaws.com/L8/139/045/LC81390452014295LGN00/LC81390452014295LGN00_B1.TIF)

Versions (installed through Anaconda):

rasterio - 0.36.0

python   - 2.7

Sorry if the above is vague. My questions would be - why does this work for Landsat-8 but not Sentinel-2? How can rasterio be used to open Sentinel-2 AWS data?


texts: Enable 'on the fly' CRS transformation in QGIS 3?

texts:
In QGIS 2.8 there is a clear button to 'enable on the fly CRS transformation' in the "project properties" for a new project.  Seen here:

But in QGIS 3.0, I cannot find it to save my life.

Can anyone help?


texts: Can't connect Docker PostGIS container to QGIS

texts:
I recently started experimenting with Docker and PostGIS following this documentation. I have a locally installed version of QGIS which I am attempting to connect the default Postgres database to. Unfortunately I have not been able to achieve a successful connection?

Here is what I am attempting to connect with:

(https://imgur.com/a/lebIT)

I have attempted both localhost and my Docker container's IP address, as well as multiple iterations of Postgres, root, and the 'mysecretpassword' login provided in the Docker documentation but still no luck.

The only other bit of pertinent info I can think of is the fact that I am running all this on a Windows 7 / Ubuntu dual boot and I have an installation of postgresql/postgis/pgadmin on my Windows partition as well as a postgresql folder in my /var/lib/ directory on Ubuntu which (I think?) is from a previous installation that I removed.


texts: Error in clipping raster (.bil) using vector polygon in qgis

texts:
I am trying to clip the required portion of raster using polygon as a mask layer. The raster image is in .bil format. While loading the raster, I chose CRS as WGS84. Then I added the vector layer for which clipping was to be done. Using "Clipper" I tried to create clip layer. But it showed certain error as "Warning : the input vector layer has a SRS, but the source raster dataset does not.
Cutline results may be incorrect."
Both raster and mask layer have same CRS WGS84


texts: 'SELECT AddGeometryColumn' in C++ crashed

texts:
I am using MSVC 2018 developed a GIS program.  I use spatialite API to create a SQLite geometry database.  The program always crashed when executing "Select AddGeometryColumn'. I don't know if I need to include a preprocessor name or miss a library in my project file.  I included the following libraries,

spatialite_i.lib
sqlite3_i.lib
proj_5_2.lib
lwgeom_i.lib
freexl_i.lib
geos_c.lib
iconv.lib

My program as the following,

CString sSRID = sSRID.Format("%d", srid);
CString qrybuf = CString("") + "SELECT AddGeometryColumn('" + recs->code + "','geom', " + sSRID + ", 'POINT', 'XY');";
char *errmsg;
int result = sqlite3_exec(sqliteDS->getHandle(), qrybuf.GetString(), NULL, NULL, &errmsg);
if (result != SQLITE_OK) {
    acad_printf("failed to add geometry column. error: %s\n", errmsg);
    sqlite3_free(errmsg);
    return;
}

Can you advise?


texts: Qgis 2.14 Changing line and buffer at the same time

texts: I have created many buffers for a line shapefile, but now I have to change the line path and would like the buffers to change automatically. Is this possible?


texts: Advice on benchmarking my Tile Server (Apache2+mod_wsgi)

texts:
I'm wondering what advice people have for benchmarking a tile server?

Specifically, I'm running a TMS using TileLite, and Vector Tiles using TileStache. Both are running under Ubuntu Server 10.10 / Apache2 / mod_wsgi, using different virtual hosts. Both are pulling data stored in PostGIS, and are caching to disk. I'm displaying these tiles using Polymaps.

I'm interested in establishing a benchmark for TileLite and TileStache, then exploring options for increasing their performance/reliability.

I've tried out using ab on static (individual tiles), but I'm unsure how to simulate a user interacting with the map, causing numerous tiles to be rendered.

Anyways, any tips or sharing of experiences that people have had would be greatly appreciated!

thanks,
matt

N.B. Apparently I can't create new tags like 'TMS' or 'Apache2' because I'm too new to the forum :(


texts: Labels not displaying for layers inside group using ArcPy

texts:
I'm using a python script to create an empty group layer in the current map document, add some layers to this group and update their symbology from an existing .lyr file.

I also want to to turn on these layers' labels using layer.showLabels = True and it works fine for normal layers, but for some reason it doesn't work on layers inside a group layer. However, when I manually drag these layers outside of the group, the labels are suddenly displayed, and they are still displayed when I drag the layers back into the group.

Why does it happen and how can I programmatically display the labels even for layers inside a group?

Here is the relevant part (of a simplified version) of my ode:

mxd = arcpy.mapping.MapDocument("Current")
df = arcpy.mapping.ListDataFrames(mxd)[0]
groupLyr = arcpy.mapping.Layer(r'C:\Path\to\Group.lyr') # this is an empty group lyr file

arcpy.mapping.AddLayer(df,groupLyr)
currentGroup = arcpy.mapping.ListLayers(mxd,"Group")[0]

source_fc = arcpy.mapping.Layer(r'C:\Path\to\my.gdb\feature_class')
arcpy.mapping.AddLayerToGroup(df, currentGroup, source_fc, "BOTTOM")

lyr = arcpy.mapping.ListLayers(mxd,"[my layer's name]")[0]
lyr.showLabels = True # this line doesn't raise an error, but it does not turns the labels on


texts: How to define new classes from a change detection raster layer

texts: As part of my thesis, I'm looking at the atmospheric and landscape effects of the 2018 wildfire in the Algarve. As part of the process, I want to estimate total biomass loss so I have classified two Landsat 8 images (pre and post-fire, 5 and 6 classes respectively) and performed a land cover change analysis using SCP plugin. I'm still fairly new to QGIS/remote sensing so am struggling to understand the output from the change detection and how I can produce a layer showing classes which have changed/calculating total biomass loss from the vegetation classes.


texts: Resources on using R in QGIS for R users?

texts:
I'm finding that much of the QGIS documentation is oriented toward use of the GUI.

I want to be able to run R scripts in QGIS, but am having trouble finding a resource that will ground me in the relationship between QGIS and R data structures and how to pass structures between the two.  A related question (Is it possible to create and run custom R scripts in QGIS 2.0?) directed me to this page: http://www.qgis.org/en/docs/user_manual/processing/3rdParty.html, but the I'm finding the documentation there is coming at things from more of a QGIS-expert direction than what I'm looking for.

Any suggestions?


texts:
I am working on NDWI and MNDWI indices using Landsat 5,7 and 8 to identify permanent water bodies in two images by mosaicing them. When using Landsat collection 1 data I have performed atmospheric and radiometric correction on all of them. But recently I want to download Landsat Collection 2 Level 1 data because I only need to work with a few bands.
Do I need to do atmospheric and radiometric corrections for Landsat collection 2 level 1 data?


texts: Channel network from high resolution DEM for flat areas

texts:
I'm working on channel network extraction from a LiDAR DEM (square grid max resolution 0,25 cm), and my area of study is vary flat with some part below the sea level.

I applied the classical valley recognition method with D8 and ArcGIS software but the output, I think, is not correct (the fill tool modifies a lot the DEM).

I saw there are various method like, for example, Pelletier (2013) and Clubb et al. (2014) but they seem to be created for mountain areas and I can't find many methods applicable in flat areas.

What procedure can I use?


texts: Getting directions north/south/east/west with Google Maps API?

texts:
for example I have 5 markers and I want to get the direction (E,W,NE etc) for each one of them according to the main one(red in the center). Here is a pic.

Is there a way to do this with google.maps.DirectionsService or calculate it somehow?


texts: osmconvert cannot convert pbf file with -b option

texts:
I've downloaded the latest UA.osm.pbf file from

  http://data.gis-lab.info/osm_dump/dump/latest/UA.osm.pbf

I've installed osmconvert 0.8.5 to my Ubuntu 16.04
osmconvert successfully opens and creates statistic for UA.osm.pbf file

osmconvert UA.osm.pbf --out-statistics
-----------------------------------------------------------------
timestamp min: 2007-02-21T08:17:13Z
timestamp max: 2016-07-28T23:31:48Z
lon min: 19.1582119
lon max: 56.0369675
lat min: 41.6472555
lat max: 56.9060947
nodes: 44162664
ways: 5058105
relations: 110626
node id min: 487747
node id max: 4326298295
way id min: 4306002
way id max: 434671917
relation id min: 2865
relation id max: 6445175
keyval pairs max: 309
keyval pairs max object: relation 60189
noderefs max: 2000
noderefs max object: way 87717954
relrefs max: 3630
relrefs max object: relation 2094111
-----------------------------------------------------------------
Finished! Calculation time: 8s.

Also it can successfully convert pbf ==> 05m

But if I try to crop it by -b option I get file 167 bytes length

osmconvert UA.osm.pbf -b=47.45,32.94,49.19,36.94 --out-pbf -o=UA.osm_01.pbf

algot@algot-16:~/mapbuild/osmand-pbf$ ls -la
-rw------- 1 algot algot       167 UA.osm_01.pbf
-rw-rw-r-- 1 algot algot 380945644 UA.osm.pbf

What I'm doing wrong and how can I fix it?


texts: How can you send by Python console an expression to a layer property and activate its corresponding check box

texts:
Try this code but it has only worked for me with fixed and short values for example 'hello' or '777'
layer = iface.activeLayer()
idx = layer.fields().indexFromName('cp')
layer.setDefaultValueDefinition(idx, QgsDefaultValue(QString str = "aggregate('manzana','concatenate', "cve_mza",intersects(geometry(@parent), $geometry),',',"manzana" )", applyOnUpdate: bool = True))

There will be some code?
I´m using QGIS 3.12.


texts: Issue in loading external geojson file into leaflet map?

texts:
I want to load an external geojson file into leaflet map. I created a global variable (i.e. tracts) for json file and referenced this variable in the script file. But i receive an error message that variable is not defined. Please check the below html and JS code.

I also tried using L.GeoJSON.AJAX call, but still the geojson external file is not loading up on the leaflet map and i get an error message that 'L.GeoJSON.AJAX is not a constructor'.

## JavaScript

var map = L.map('mapid', {
center: [5, 28],
zoom: 3,
minZoom:2,
maxZoom:18
});

L.tileLayer('http://{s}.tiles.mapbox.com/v3/ianmule.bg2v5cdi/{z}/{x}/{y}.png', {
 attribution: "Mapbox"
}).addTo(map);

L.geoJson(tracts, {

style: style,
onEachFeature: onEachFeature
}).addTo(map);

# Using Ajax call.
 var mystyle = {
"color": "#ff7800",
"weight": 5,
"opacity": 0.65
 };

var geoJsonLayer = new L.GeoJSON.AJAX("kc-tracts.js", {
style: mystyle });
geoJsonLayer.addTo(map);


texts:
I am looking for a way to compare routes of two or more persons to be taken by car.

Person 1 has a starting point A1 and a destination point A2.
Person 2 has a starting point B1 and a destination point B2.
[...]

What I would like to accomplish:

Check, if A1 and B1 are close to each other.
Calculate walking and/or driving distance from A1 to B1.
Save the individual routes of Person 1 and Person 2 in a way, that allows for performant comparison to other routes / route points.
Identify, if the routes to be taken of Person 1 and Person 2 are similar enough, to have at least same route parts.

Are there already (open source) examples of code around that I haven't found yet?


texts: Determining errors in latitude and longitude?

texts: I have a large database that contains at grade railroad crossings (railroad crossings that intersect roadways). Many of these contain incorrect latitude and longitude values. I need to determine which coordinates are incorrect without entering them all in google maps (there are thousands of crossings). These coordinated need to be accurate, but an error up to 200ft is okay.


texts: Join is working in pgAdmin but not in DB Manager

texts:
I have two tables tbl1 and tbl2 without geometries which I'm trying to join (and later join result to a table with geometries when I get the first join to work):

Simply:

select * from tbl1
left join tbl2
on tbl1.status_id = tbl2.status_id

Is working in pgAdmin:

But not in (QGIS) DB manager:

It is driving me crazy. What could be the reason for this?


texts: Raster stretching on ImageMosaic with Geoserver

texts:
I am trying to demo serving rasters over GeoServer.  I have created the store and the layer but the layer preview comes up all black.  The layer is queryable however so I know the data is accurate because I see the correct values when I click on the raster.

The raster values are single-band grayscale from 0 to 1.  How do I stretch the colors on these from black to white the same way they would be if loaded into QGIS or ArcMap?  Is it possible to change color ramp altogether to pseudocolor?


texts: OGR set different features to column name in PDF

texts:
I'm trying to import shapefile data into a geospatial PDF using OGR and Python.

My shapefiles are structured as follow: shapefile1.shp: {stratt1: string1, intatt1: 1, ...} and shapefile2.shp: {stratt2: string2, intatt2: 2, ...}

From the documentation and other answers in SO I found the OGR_DISPLAY_FIELD=name command, but it only works for one shapefile. Thus, if I import several shapefiles at once, I get only feature1, feature2, feature3 and so on at the feature tree, as it can be seen in the following image:

This is odd, because with OGR_DISPLAY_LAYER_NAMES=names I'm able to change the layer name of several shapefiles at once like OGR_DISPLAY_LAYER_NAMES=  "points1, points2", but I had no luck finding a way to do the same for the features. Any ideas?

Because I'm doing it with Python, my code looks like this:

co = ["OGR_DATASOURCE=" + vrt, "OGR_DISPLAY_FIELD=" + "sat1, sat2"]
gdal.Translate(outpath, imagery, format="PDF", outputSRS="EPSG:4326", creationOptions=co)

(Again, the part failing is the "OGR_DISPLAY_FIELD=" + "sat1, sat2" one)

My VRT file (called vrt above):

<OGRVRTDataSource>
    <OGRVRTLayer name="points1">
        <SrcDataSource relativeToVRT="1">shapefile1.shp</SrcDataSource>
        <SrcSQL>SELECT stratt1 AS sat1, intatt1 as inat1 FROM shapefile1</SrcSQL>
        <LayerSRS>WGS84</LayerSRS>
    </OGRVRTLayer>
    <OGRVRTLayer name="points2">
        <SrcDataSource relativeToVRT="1">shapefile2.shp</SrcDataSource>
        <SrcSQL>SELECT stratt2 AS sat2, intatt2 as inat2 FROM shapefile2</SrcSQL>
        <LayerSRS>WGS84</LayerSRS>
    </OGRVRTLayer>
</OGRVRTDataSource>

Any ideas how could I achieve my goal with only one gdal.Translate() call?

TLDR:

In the screenshot of the PDF you can see that every feature in points1 is displayed as string1, the attribute from the column stratt1 in shapefile1. This is correct.

In the screenshot of the PDF you also can see that every feature in points2 is displayed as feature1, feature2, feature3... and so on. This is NOT correct. It should be called string2 following the same principle as above.

I can´t do "OGR_DISPLAY_FIELD=" + "sat1, sat2" because OGR_DISPLAY_FIELD=name doesn´t take lists, nor it is possible to call it two times in the creationOptions of gdal.Translate()

Edit to @cm1's answer


texts: Visualizing large datasets of polygons with Leaflet

texts: I am using leaflet to visualize several layers in the project. One layer(load from geojson format) contains thousands of small polygons, which has a very slow loading speed. I read this Visualizing large datasets with Leaflet, but it mainly deals with the markers. Is there a similar way to deal with polygon?


texts:
I am linking from QGIS to a FileMaker Pro record using links in HTML map tips. I can call the custom url scheme fmp:// with an IP address without issue, but as soon as I include the dollar sign ($) in the URL, the link no longer works.

I've tried encoding the symbol to %24 and concatenating strings but they don't work. I've also tried &#36; as suggested here, and variations of \$, \u0024 and \%24 with and without leading single and double backslashes.

Examples given by FileMaker:

fmp://$/Clients.fmp12?script=ListClients
fmp://192.168.10.0/Clients.fmp12?script=ListClients

The $ is significant, as it calls the currently open database, rather than trying to reopen the file, and removes the need to hardcode an IP address.

Update

I've double checked that $ is a valid symbol to include in a URL as mentioned here. I also managed to get a log output when calling the url: invalid escape sequence \$

qt seems to offer a range of options for encoding urls and strings, with reserved, unicode and byte codes for characters. I'm not sure which QGIS is using internally.

Even typing backslash escapes as a string in the expressions editor turns $ into the ? symbol.


texts: Is there an comparison of the various mapping libraries for Android?

texts:
I am currently looking at the various libraries that I can use to show a map in my custom Android App. I have found several like Mapsforge, OSMDroid, ArcGIS runtime for Android, Google Maps for Android and so on. I'm sure that I am missing a few.

Is there a comparison chart somewhere about the differences in functionalities and pros and cons of each one?

If there isn't one, we could make this question community wiki, to serve alongside Comparison of Open Source Desktop GIS Packages & Comparison of JavaScript mapping libraries.


texts: Create a ee.FeatureCollection out of a .tif file

texts:
I am using the Earth Engine API in a Python Jupyter notebook.
I have several .tif files in my local folder and I would like to create ee objects (FeatureCollection a priori) to manipulate them. Is possible without adding them in my assets?


texts: Calculating monthly averages from daily precipitation data

texts:
I have a raster brick for daily precipitation data and I need to convert it to monthly data, how would I go about doing this? Each layer is a day and there are 34698 days (i.e. too many to go through manually and create monthly averages by listing mean(raster1, raster 2... etc)
class      : RasterBrick

dimensions : 64, 58, 3712, 34698  (nrow, ncol, ncell, nlayers)

resolution : 0.44, 0.44  (x, y)

extent     : 12.1, 37.62, -45.54, -17.38  (xmin, xmax, ymin, ymax)

crs        : NA

source     : C:/Users/teyah/AppData/Local/Temp/RtmpQBzpnO/raster/r_tmp_2020-08-13_121709_23408_66253.grd

names      : Precipitation.1, Precipitation.2, Precipitation.3, Precipitation.4, Precipitation.5, Precipitation.6, Precipitation.7, Precipitation.8, Precipitation.9, Precipitation.10, Precipitation.11, Precipitation.12, Precipitation.13, Precipitation.14, Precipitation.15, ...

min values :    0.000000e+00,    0.000000e+00,    0.000000e+00,    0.000000e+00,    0.000000e+00,    0.000000e+00,    0.000000e+00,    0.000000e+00,    0.000000e+00,     0.000000e+00,     0.000000e+00,     0.000000e+00,     0.000000e+00,     0.000000e+00,     0.000000e+00, ...

max values :    6.533941e-04,    1.833772e-04,    8.561218e-04,    1.070062e-03,    7.177751e-04,    6.058318e-04,    5.470571e-04,    1.588907e-03,    1.674356e-03,     1.091221e-03,     5.895557e-04,     3.984022e-04,     4.418955e-04,     3.041818e-04,     2.085147e-04, ...


texts: how to access Canvas mouse point on action of button from another window Pyqt5

texts:
I am PyQt5 and PyQgis novice. Below code snippet is my main window application. Below code throw a main Window with QGsMapCanvas as central widget as canvas, which displays maps.

On menu action and i can display Pyqt5 window , there are 4 such GUI. all are runnning successfully. But..  Mainwindow (canvas) and GUI (from menu bar action) are two different window.

I want to access polyline or point or polygon drawn in canvas on Action of QPushButton in GUI.

that is how to access mouse points or QgsPointXY of Canvas to be shared with other Pyqt5 GUI (on action of button).

Note: i am able to draw line in canvas using custom LineTool (QgsMapToolEmitpoint) from toolbar action in mainwindow, but can not shared with other GUI on button event.

import sys,os
from qgis.core import *
from qgis.gui import *
from PyQt5.QtCore import *
from PyQt5 import QtGui, QtWidgets, uic
from PyQt5.QtGui import *

# Environment variable QGISHOME must be set to the install directory
# before running the application
qgis_prefix = os.getenv("QGISHOME")

class Ui(QtWidgets.QMainWindow):
    def __init__(self):
        super(Ui, self).__init__()
        uic.loadUi('MainMenu.ui', self)
        self.setWindowTitle("Standalone Desktop Application")
        self.show()

    self._mks=None
    self._mkb=None
    self._mkw=None
    self._mka=None

    # Find & connect Menu1
    self.item1 = self.findChild(QtWidgets.QAction, 'menuItem1')
    self.item1.triggered.connect(self.showGUI1) # Remember to pass the definition/method, not the return value!

    # Find & connect  Menu2
    self.item2 = self.findChild(QtWidgets.QAction, 'menuItem2')
    self.item2.triggered.connect(self.showGUI2) # Remember to pass the definition/method, not the return value!

    # Find & connect Menu3
    self.item3 = self.findChild(QtWidgets.QAction, 'menuItem3')
    self.item3.triggered.connect(self.showGUI3) # Remember to pass the definition/method, not the return value!

    # Find & connect Menu4
    self.item4 = self.findChild(QtWidgets.QAction, 'menuItem4')
    self.item4.triggered.connect(self.showGUI4) # Remember to pass the definition/method, not the return value!

    # Find canvas object
    self.canvas = self.findChild(QgsMapCanvas, 'canvas')
    self.canvas.enableAntiAliasing(True)

    # Load A Test Vector layer and display in canvas
    self.loadGISLayer()

def showGUI1(self):
    # This is executed when the menu item is pressed
    print(' GUI 1 To be Launched')
    self._mks=showUI1()
    self._mks.show()

def showGUI2(self):
    # This is executed when the menu item is pressed
    print(' GUI 2 To be Launched')
    self._mkb=showUI2()
    self._mkb.show()

def showGUI3(self):
    # This is executed when the menu item is pressed
    print(' GUI 3 To be Launched')
    self._mkw=showUI3()
    self._mkw.show()

def showGUI4(self):
    # This is executed when the menu item is pressed
    print('  GUI 4 To be Launched')
    self._mka=showUI4()
    self._mka.show()

# Class for UI1
class showUI1(QtWidgets.QMainWindow):

    def __init__(self):
        super(showUI1, self).__init__()
        uic.loadUi('UI1.ui', self)

# Class for UI2
class showUI2(QtWidgets.QMainWindow):

    def __init__(self):
        super(showUI2, self).__init__()
        uic.loadUi('UI2.ui', self)

# Class for UI3
class showUI3(QtWidgets.QMainWindow):

    def __init__(self):
        super(showUI3, self).__init__()
        uic.loadUi('UI3.ui', self)

# Class for UI4
class showUI4(QtWidgets.QMainWindow):

    def __init__(self):
        super(showUI4, self).__init__()
        uic.loadUi('UI4.ui', self)

app = QtWidgets.QApplication(sys.argv)
QgsApplication.setPrefixPath(qgis_prefix)
QgsApplication.initQgis()

print(qgis_prefix)
print("------------------------------------------------")
print(QgsApplication.showSettings())

window = Ui()
app.exec_()


texts: Export to images subsets of Spatial Points Data Frame

texts:
I have a csv file, which I converted it to Spatial Points Data Frame. And I want to split the spatial data and export each map to a jpg file.

I did the following:

SB <- read.csv("SB_040914.csv", header=TRUE, sep = ";")
SB <- SB %>% group_by(Fecha, CodPar, x, y) %>% summarise(mean(Sube),mean(Evade), mean(Baja))
SB <- data.frame(SB)
SB <- SB[complete.cases(SB),]

coordinates(SB) <- ~ x + y
proj4string(obj = SB) <- CRS("+init=epsg:32719")

SB <- split(SB, f = SB$Fecha)

The Split function, shows me the next warnings:

Warning messages:
1: In min(x) : no non-missing arguments to min; returning Inf
2: In max(x) : no non-missing arguments to max; returning -Inf
3: In min(x) : no non-missing arguments to min; returning Inf
4: In max(x) : no non-missing arguments to max; returning -Inf

And for trying to export I did this:

for (i in SB) {
  jpeg(paste(SB[[1]],".jpg",sep=""))
  bubble(SB[[i]], "mean.Evade.", col=c("#00ff0088", "#00ff0088"),
         main = "Evasión" )
  dev.off()
}

Which shows me the next error:
Error in as.character.default() :
  no method for coercing this S4 class to a vector
6 as.character.default()
5 as.character()
4 paste(SB[[1]], ".jpg", sep = "")
3 gsub("%%", "", s)
2 checkIntFormat(filename)
1 jpeg(paste(SB[[1]], ".jpg", sep = ""))

As you can see, this is made in R.

The link to the data: https://www.sendspace.com/file/nd6j8h


texts: Is there something like ArcGIS's Attribute Assistant available in QGIS?

texts:
I'd like to update attributes from an excel table based on primary/foreign key relationships.

In ArcGIS´s attribute assistant, this method exists:
Copy Linked Records Updates an attribute of a feature with a value from a related table.


texts: Find polygons that has no points in PostGIS

texts:
I have a layer with points and a layer of polygons. I can find all the polygons that contains at least one point inside itself by:
SELECT polygons.id
FROM polygons, points
WHERE ST_Intersects(polygons.geom, points.geom)

But what if I just want the opposite of this: select all polygons that contains no points?
SELECT polygons.id
FROM polygons, points
WHERE NOT ST_Intersects(polygons.geom, points.geom)

This last query actually returns all polygons since the condition NOT ST_Intersects(polygons.geom, points.geom) can be true for some points/polygons pairs (there is always a point somewhere which is not inside a polygon).


texts: Trying to set the OpenLayers MouseOver position to be outside the Map

texts:
This is my code and I am trying to obtain the MouseOver to be displayed outside the Map container.

<!DOCTYPE html>

        <title>My Map</title>
        <link rel="stylesheet" href="openlayers/theme/default/style.css" type="text/css">

My Map

Position

I am using the , the mouse over co-ordinates are only appearing at the bottom right of the map (within the map container).


texts: How to best execute the same query on many different vector layers?

texts: I have some vector layers in my project and I need to separate the polygones by area (I need a copy of all shapes > 2ha). I know how to do this for each separate layer (query) but is the a tool to do this for all my layers in one step?


texts: Geoserver WMS limiting map visibility by scale

texts:
In ArcGIS I can manually limit visibility of a map layer by specifying a min and max scale (see below).

Is it possible to instead set this server-side on Geoserver? I've seen something similar in Tile Caching but didn't see it enforced when i tested it in ArcGIS.


texts: How to fix Excel/Arcmap "data type mismatch" error?

texts:
I've had a polygon shapefile joined to an XLS file with no issue until today. Although the category of the column category is designated as "Number" in Excel, ArcMap recognizes the column as "String" when the table is added into Arc or when a Join is performed. This is good since it allows for the Join to be performed.

However, something must have happened today as Arcmap now reads the column as "Double". This of course prevents the Join between the shapefile and XLS. I suspect I entered some value incorrectly, but I can't seem to find the problem cell.

I firstly tried changing the column's category in Excel from "number" to "text", but that didn't solve it, nor did changing the category to "general". So I then created a new column ("B", whereas the column I'm joining with the shapefile is column "A"), and used the formula " =TEXT(A2,"0") " . ArcMap will now allow me to join the shapefile to this newly created column.

Although the problem is solved, I want to pinpoint which cell may have causing the original problem. I have thousands of 5 digit numeric values in this column which were joining with no issue before today.

My process is that I work from an XLSX file and then, when I want to see the updates reflected in the joined shapefile, I close Arc (to remove the schema lock from the XLS), and then save the XLSX to XLS. My XLSX file does have a "highlight duplicates" conditional formatting rule, but I don't think that was the problem since it has otherwise worked up to this point. Additionally, the Excel file has a "custom sort" option with two levels of sorting rules, but again that hasn't caused any issues before.

So I just want to know whether there is some way to check for which cell is causing Arc to read the column as "double" as opposed to "string" as it was before.


texts: rgdal / ggmap : rotate map data and clip map

texts:
I'm trying to plot a map of Toronto with ward boundaries and geocoded data on events over years.  But the layout is such that I end up with mostly wasted space in the map.  To increase resolution for the data I would prefer to rotate all the data about 12-15 degrees clockwise to make the top boundary approximately horizontal and clip the map to a bit more than the city boundary shown below.

The ward boundaries come from http://www1.toronto.ca/City_Of_Toronto/Information_Technology/Open_Data/Data_Sets/Assets/Files/wards_may2010_wgs84.zip and are read in using rgdal:

dir <- "C:/Dropbox/Documents/consulting/Abigail"
folder <- "wards_may2010_wgs84"
library(rgdal)
wards <- readOGR(paste0(dir,folder), layer="icitw_wgs84")
plot(wards, col="gray95")

I get the map for my main plot with get_map:

library(ggmap)
map <- get_map("Toronto", zoom=10)
p <- ggmap(map)
# convert wards to a data.frame for use with ggplot2/ggmap and plot
wards.df <- fortify(wards)
p +
    geom_point(data = section37, aes(x = long, y = lat, group = Ward), size=0.7, color="red") +
    geom_polygon(data=wards.df, aes(x = long, y = lat, group=group),
                 color="blue", alpha=.1)

I'm fairly new to mapping applications in R, so could use some help with this task.


texts: What is the order of raster after stacking them?

texts:
When I merge a lot of raster files I don't know what is the order that they forward in the new file (stacking).

How do I know the name of band 1 (and others ones) in the other file that I have had created.

Is there a default rule to it?


texts:
I built my own tile server and when I try to prerender zoom level 10 or 13 modtile/renderd just stops working (like when I press Ctrl + C in the terminal).
The other levels work fine.

OS: Lubuntu 14.04

Style: OpenTopoMap

Command to start modtile/renderd:

sudo -u renderuser renderd -f -c /usr/local/etc/renderd.conf

Command to start prerendering (for example lvl 10):
sudo -u renderuser render_list -a -z 10 -Z 10 -s /var/run/renderd/renderd.sock

So I start modtile and after that render_list. It starts doing stuff and rendering some tiles. But after some time both stop working, as you seen in the following screenshot:

Is there a way to see why it stopped?

Has somebody faced the same problems?

Has somebody a idea how to get it working?


texts: Arcmap Update field from a different table

texts:
I want to update a field in one table using two matching values in another table. Let's call the fields A, B and C in one table, and X, Y and Z in the other table.
So, where A=X and B=Y then update C with Z


texts: Merging small polygons into larger adjacent polygons using QGIS

texts:
I wish to select just the small green polygons (see image) and merge these with the pink polygon. However, I wish to leave the large green polygons that are over 1km2. However, all the green is one polygon, when I select one small area every green polygon is selected.
Any ideas?


texts: How to iterate through maps & add new layer to data frame and legend of each using ArcPy?

texts:
Newbie for Python and wondering if the following is doable.

Problem is like this:

I have two sets of data:

one mxd file for each of multiple cities (all stored in the same folder) containing links to multiple layers of data.
additional dataset (multiple fields; point shapefile) with data for all cities stored in a file geodatabase (*.gdb)

My goal is to have a python script to:

open all mxd file iteratively
open the additional dataset and select the data based on a query (City name) and add to the mxd file.
2) Label the added data in the map by selecting desired field
2) add the legend of the added data at the bottom of the original legend.

Help with scripting would be highly appreciated.


texts: Slow retrieve PostGIS raster value through ST_Value

texts:
I have a PostGIS DB containing a set of daily raster data from 1981 up to 2021. The inserted raster tile size is 83x90.
Problem: I wanna query the data through ST_Value for a specified Lat Lon, however, I notice it's so slow according to the coordinate I use. The first query execution takes so long (around 1.5 minutes) to be executed but if I query again passing a Lat Lon nearby to the latest query, it is faster (around 10 sec).
Basically, this is the way I'm querying.
SELECT
    date,
    ST_Value(rast, 1, (ST_Transform(ST_SetSRID(ST_MakePoint(-119.167504,44.962215 ), 4326),4326)))
FROM
    myraster
WHERE
    AND date >= '19810101' AND date <= '20210101'
    AND ST_Intersects(rast, ST_Transform(ST_SetSRID(ST_MakePoint(-119.167504,44.962215), 4326),4326))

Table:

column
type
Nullable

id
bigserial
not null

variable
character varying(10)
not null

date
date
not null

rast
raster
not null

filename
text
not null

Explain analyze result:
1. Index Scan using myraster_rast_gist on myraster  (cost=0.41..7.18 rows=1 width=12) (actual time=4.881..60347.936 rows=14618 loops=1)'
2. Index Cond: ((rast)::geometry && '0101000020E61000007D04FEF0F39354C0A25D85949FA03D40'::geometry)'
3. Filter: _st_intersects('0101000020E61000007D04FEF0F39354C0A25D85949FA03D40'::geometry, rast, NULL::integer)
4. Planning time: 0.272 ms
5. Execution time: 60353.772 ms

Looks like the raster structure maintains a kind of cache for each query, where,  as farther away is the coordinate from the latest fetched point, the query takes so long to be executed. I'm trying to figure out a way to reduce the 1.5 minutes.


texts: from geometrycollection to polygon geometry with python

texts:
I have a geojson file containing a feature with geometrycollection type (polygon and linestring).
I need to produce two ESRI Shapefile from the geojson geometrycollection: one with Polygon only and a second with Linestring only.

The code must be written in python.

There is any solution?

I met another problem in the "from geojson to shapefile" process.
My geojson file include some features with geometrytype = MultiPolygon.
These MultiPolygon are constructed in such a way that the first element includes the following elements and these never intersect.
I enclose a picture , showing the geojson file opened into QGIS, that clarifies the situation:

When I use the ogr2ogr to convert geojson to shapefile I obtain the following result, shown in the following figure:

There is an ogr2ogr function allowing to explode the multipolygon into polygon?

I can't use the QGIS ftools "from multi parts to single parts" because I must have the instance of python I use to develop separeted from the QGIS folder.
I try with fiona, shapely and json:

    with fiona.open(settings.L_json, 'r') as data_file:

    data = json.load(data_file)
    for f in data_file:
        next = f.next()
        multipol = shape(next['type']['geometry']['properties'])
        for poly in multipol:
            feature = {}
            feature['type'] = poly['type']
            feature['geometry']['type'] = 'Polygon'
            feature['geometry']['coordinates'] = poly['geometry']['coordinates']
            feature['properties'] = poly['properties']
            data['features'].append(feature)
            if settings.L_json['*L_*']:
                with open('G:\\geojson\\geojson_Rep\\LOW_'+settings.now+'.geo.json', 'w') as f:
                    f.write(json.dumps(data))

but it does not work.
Any simplest solution? ogr2ogr can do it? It is ok also obtaining from the first geojson other geojson splitting the multipolygon geometry.

Simplifying, I need to get n polygon from a multipolygon with n elements, in one of the following formats: .geo.json, .shp, .wkt.

Any hint?

Here my solution:

with fiona.open(input, 'r') as source:

    count = 0
    for f in source:
        c = f['geometry']['coordinates']

        if count == 0:
            outputfile = 'G:\\SVN EO\\Python\\Utility\\geojson\\tmp\\mondometeo_20160930_145814_MEDIUM.geo.json'
            with fiona.open(outputfile, 'w', driver=source.driver, schema=source.schema) as output:

                output.write({
                    "type": f['type'],
                    "geometry": {"type": "Polygon", "coordinates": c[count]},
                    "properties": f['properties']})
            return outputfile
        elif count >0:
            outputfile = 'G:\\SVN EO\\Python\\Utility\\geojson\\tmp\\mondometeo_20160930_145814_TOP'+str(count)+'.geo.json'
            with fiona.open(outputfile, 'w', driver=source.driver, schema=source.schema) as output:
                output.write({
                    "type": f['type'],
                    "geometry": {"type": "Polygon", "coordinates": c[count]},
                    "properties": f['properties']})

            return outputfile
        count +=1

Anyone has some other solution or hint (also performing the if-elif cycle)?


texts: Swapping start and end points of line (i.e. direction) of LineString in PostGIS?

texts:
I know that postgis has the function to find the start and end point of a line.
Is there a function that can swap the start point and end point?

(Note: I do not mean literally changing the co-ordinates).

You could also describe it as switching direction of a linestring.

Or in technical terms I think we just have to reverse the sequence of coordinates.

I know a gis tool which does this, but is there postgis counterpart?


texts: GeoTIFF: How to remove Internal Overview without recompression

texts:
"gdaladdo -clean" does not remove the internal overview from the file see http://trac.osgeo.org/gdal/ticket/4940

How do I remove GeoTIFF overview(s) without having to re-compress the main image.


texts: Exporting polygon features to individual image pngs

texts: I have a shapefile with a few thousand polygons, ArcGIS [9.3] (and even PostGIS). I would like to export each of those polygons into its own thumbnail image file without any real context. What is the best approach?


texts: Checking for equations used in ArcGIS field calculations?

texts:
I was using the field calculation to seek some information. I closed my table of attributes and reopened it and wanted to verify that I had used the right equation for my field statistics. Although the correct information was still in the table, I wasn't able to see the equation I had used. Is there a way for me to retrieve or find the equation?

I'm using ArcMap Desktop, Version 10.2.2.


texts: Performing IDW interpolation and mean estimation on county level using ArcGIS Desktop?

texts:
I am absolutely new to GIS.

I have imported a map with information in the counties of Germany. Additionally I have added a layer with information on SO2 concentration at various measuring stations in Germany:

Additionally, I have already interpolated the SO2 data using IDW (based on the suggestions from the Federal Environmental Agency, who distributed the data):

I am now struggeling with calculating area means. That is, for my analysis I need one value of SO2 concentration for each county in Germany, which I can export to R.


texts:
I have a question connected with ArcPy. I was trying to make a nested loop with 2 cursors  - 1st with 2 elements and 2nd with 5 elements. The statement should return 10 connections, e.g. 1-1, 1-2, ... 2-5.
Here is part of my code:

import arcpy
import os.path

folder = "Z:\WOLGA"
na="na"
baza = "testy.mdb"
dataset = "network_analyst"
arcpy.env.workspace = os.path.join(folder, na, baza, dataset)
shapePocz = os.path.join(folder, na, baza, dataset, "dom")

listaPktPocz = arcpy.SearchCursor(shapePocz, "", "", "", "")
shapeKonc = os.path.join(folder, na, baza, dataset, "punkty_docelowe")
listaPktKonc = arcpy.SearchCursor(shapeKonc, "", "", "", "")
shapeKoncName = arcpy.Describe(shapeKonc).shapeFieldName
shapePoczName = arcpy.Describe(shapePocz).shapeFieldName

for pktPocz in listaPktPocz:
    for pktKonc in listaPktKonc:
        print str(pktPocz.numer) + " - " + str(pktKonc.numer)

For now I got such connections as: 1-1, 1-2, 1-3, 1-4, 1-5 only. The outer loop is not being executed properly.


texts: How do I convert an Est file to a shapefile?

texts: I need help converting an Est(Streets and Trips) file to shapefile or some other format that ArcGIS can read. I haven't been able to find any luck using Google and my go to program Alteryx can't read Est files. Does anyone have any insight on how I can convert this data format?


texts: How do you find the appropriate spatial weight to use for a spatial analysis?

texts:
I am currently about to run a spatial regression (OLS, Lag, Error etc) for my dataset. My dataset contains 1000 randomly generated points throughout the utilisation distribution for my animals (so, there are more or less high and low values attached to whether an animal goes to this area frequently). Since the species that I study has a very specific core range (where most of its time is spent), there will be hotspots (used often) and cold spots (used infrequently). I hope to determine whether their utilisation distribution is primarily in response to predation risk, the location of "competitive" individuals (i.e. other animals) as well as food availability.

In so far, I ran an OLS, Lag Model (there is spatial auto-correlation in my response) and an Error model in the package spdep in R. The Error model usually has the lowest AIC. In addition, I made spatial weights based off of euclidean distance in Geoda.

One thing that has really been bugging me since I have yet to find a guide on this, and I am rather new with spatial analytics involves well...choosing the right weights (i.e. K neighbour selection, queen, rook, distance). It should be noted that since I am dealing with 1000 random points for this analysis...and, they are throughout the landscape (for the animals home range) there are obvious "hot spots" of clusters for high and low values. Attached is a photo from geoda that show where these values are.

Anyways, I apologise for the daft question. However, given the fact that I have the following:

animal home range data (0.1-99) that were then made into 1000 points (my response)
these 1000 points show clear clusters (primarily at feeding and sleeping sites)
predictors (attached to these points through arcgis) including food availability and predation.

I am truly concerned about the appropriate spatial weights to use given this. Digging into the literature hasn't really helped and to be honest, I am beginning to become more confused than informed!


texts: How to visualize point accuracy with buffers in QGIS?

texts:
I have a CSV file that contains 3 columns - Lat, Long, and Accuracy.  The Lat and Long fields contain point data in decimal format.  The Accuracy fields contain a numeric values that indicate the accuracy of the corresponding point.  i.e. the value '50' would indicate that the point values were accurate to within 50m.

I can put the points onto my map, but would like to also show how accurate each point is by varying the size of the point marker according to the value in the accuracy field.  In other words if a point has an accuracy value of 50, the point marker should have a 50m radius.

I am very new to GIS and have been trying to nut this out for a few days now with no success.  I have found several references to using the style options to change points according to some value, but none that will allow me to control for the exact size of the point marker.


texts: Calculating length of polyline defined by list of coordinates using Python

texts:
I have a list of coordinates, knowing that I am working using EPSG:4326 and using Python I would like to calculate the length in kilometres of the polyline defined by these coordinates.
Coords = [(0.0, 50.787944),
 (0.0, 50.787944),
 (-0.20159865271498856, 50.824569950535725),
 (-0.40044683717220364, 50.803694),
 (-0.599761967889834, 50.78975316549538)]

I've tried with GeoPandas and shapely but reading the documentation I can't find an easy way to do it. Is there any function that allows me to calculate the length in kilometre knowing Coords and the reference system used?


texts:
When we convert data from bessel1841 to WGS84, parameter setting must be used because it is impossible automatic convert on QGIS.

Reference: How to convert a shapefile from bessel1841 to WGS84?

Is it possible to convert data from WGS84 to Bessel1841?
Your advice will be big support for me.


texts: What's the parameters of heatmap formula on QGIS?

texts:
I want to calculate the values by hand to know which parameters effect the result in the image, but i couldn't able to do that because i didn't understand it. I tried with 500 mt and 1000 mt radius and here are the results and i want to know the correlation between radius, cell size and the result.

This is 500 mt radius.

and this is the 1000 mt radius results.


texts: While dragging a map in Leaflet it marks a point for the feature

texts:
In my Leaflet application, When I draw a feature (Line/Polygon) and if I need to drag the base map to mark a point outside the visibility of the map then Leaflet marks additional point where I tap on the map to drag, which is not the expected behavior.

Code for marking line in typescript:

this.l_PolylineDrawer = new L.Draw.Polyline((this.map as L.DrawMap),{showLength : true,icon : this.icon,
      shapeOptions: {color:'#ffcc33',stroke:true}});


texts: I have a GP service that generates a PDF report in reportlab based on a country clicked on a web map. I've placed the HTML/Javascript document on our web server for testing, and it's working (although we can only generate one report at a time, but that is a separate issue). I've tried to place this code in an HTML widget on our website and test it, and I keep getting the "XMLHttpRequest cannot load  No 'Access-Control-Allow-Origin' header is present on the requested resource." I've tried to follow every thread I've found on this (there are dozens) and added "add name="Access-Control-Allow-Origin" value="*" to the web config files on our proxy server and the arcgis web config in the wwwroot. I'm still getting this error when I run the tool. We have two other pages on our website that use a GP tool to generate charts and these pages do not throw up this error. I'm not sure how to get this working.


texts: QGIS3.10 - Deleting duplicated and overlapping geometries based on attribute comparison

texts:
I'm having a bit of a headache on this one, I hope you'll be able to help me out. (I'm on QGIS 3.10)
Problem
Here is the case : the national Dutch zoning plan regulation is delivered with shapes in which the maximum building height is stated. Unfortunately, they keep in the dataset current and previous zones. I would like to keep only the ones that are applicable on a national scale.
I have multiple features (in my case buildings) that are overlapping and/or being duplicated in one single layer. When there is an overlap, I would like to keep the geometries that are with the attributes the closest to current date.
My understanding would be that overlapping geometries should be compared and the one with the date attribute closest to today's date should be kept.
Done so far
But so far no success, with multiple methods I could try :

Filtering out by dates doesn't make sense as depending on where you
are in the country, the latest plan might be more or less old.
There is no exhaustive list of attributes which could make the trick by
filtering out the ones that are considered outdated by the
municipalities.
Join Attributes By Location doesn't work in this case.
Delete Duplicates works but there is no control over which one is deleted, meaning that some times it is an outdated rule that is being kept.
Vector Overlay hasn't worked neither.

As the data I'm dealing with is country scale, I cannot rely on manual sorting... :/
Any help would be AMAZING!


texts: Meaning of "feature has null Z values" when using Move tool of Editor in ArcMap?

texts: i am trying to move a polygon from a location to another location. when using the move tool (editor) in ArcGIS, I place the x and y values to where i will move the polygon. however, when pressing enter, a message stating that the feature has null z values appear. how do I solve this problem?


texts: Disabling map in ArcGIS API for Javascript

texts:
Here is my code for my map:

 omap = new Map("overviewMap", { extent: initExtent, logo: false, slider: false });
tmoverview = new ArcGISDynamicMapServiceLayer(mylayer, { "opacity": 1.0, "visible": true });
omap.addLayer(tmoverview);
omap.disableMapNavigation();
omap.disableKeyboardNavigation();
omap.disablePan();
omap.disableRubberBandZoom();
omap.disableScrollWheelZoom();

Yet, when I run it, the map navigation is NOT disabled, the keyboard navigation is NOT disabled, the pan is NOT disabled, the rubber band zoom is NOT disabled, and the scroll wheel also continues to operate.  Please tell me how to disable all of this functionality in my web map.


texts:
I just made some experiments with Matplotlib and what i noticed is that all that is generated using matplotlib appear in a new window that doesn't seem to be tightly linked to qgis. The only "link" i see is a save button enabling plot export as raster images.

My aim would be to stack(/merge) graphics/plots or (georeferenced) shapes generated with matplotlib with layers in a qgis composer.

Is there a way to do that ? (i mean a tighter/quicker way than adding image objects in the composer).

Would it be a way to share a common "coordinate system" between map objects and rasters generated by Matplotlib ?

If possible could someone point me in the right direction ?


texts: Creating batch file from ERDAS IMAGINE Spatial Modeler?

texts:
How to create batch file from the ERDAS IMAGINE Spatial Modeler?

I want to do the same process to different input images, something similar which we can do by batch processing of the existing programs in Erdas.

If there is no inbuilt function then one can use model maker to make *.gmd (graphical model) file but how to batch it?


texts: Enterprise Geodatabase Replication for Offline Use

texts:
I need to have a consistently updated geodatabase that will go offline in Collector.  I am open to new ways to do this.  What I have been working on is below.  Our environment: SQL Server 2012, ArcGIS 10.2.1, ArcGIS Server 10.2.2., Collector on iOS.

Our production sde database does not meet the requirements for offline use so a new database has been created solely for the offline purpose, "OfflineDBProd".  Since the script we use to update the offline database takes a few hours to run and in an effort to have OfflineDBProd down for as little time as possible, a staging database, "OfflineDBTemp," was created on a different server instance for this purpose. In reality, both OfflineDBProd and OfflineDBTemp have the same name but to be less confusing here I have them named differently.

A script is used to, daily, copy and transform data from the production sde database to OfflineDBTemp.  Contents of OfflineDBTemp are deleted and is refilled with the latest data, archiving enabled, and global IDs added, among other things.

A different script backs up OfflineDBTemp, copies the .bak file to the production server, and OfflineDBProd is restored with OfflineDBTemp.

To me it seems like the only changes in the OfflineDBProd (as viewed from the Collector app) would be the edits made in the production sde database the day before.  Collector will download a fresh copy of the data but will not sync the following day, after the script and restore was carried out.  The error message is:

  Invalid or missing input parameters.  Replica with this GUID does not exist on the server.


texts: When generating WMS service in ArcGIS random numbers are added to the layer name

texts:
I have an mxd file. Each time I share/publish a WMS service to the ArcGIS server as WMS, random numbers are added to the layer name. For example here is added number 8950 at the end of the <Name> tag. Even though the original name is down in the title and there is no number at the end. I pair these services with other services by layer name and need to prevent ArcGIS to add random numbers to the names and let them stay as they are.
<Layer queryable="1">
<Name>Povodeň_z_topiaceho_sa_snehu_gc8950</Name>
<Title>
<![CDATA[ Povodeň z topiaceho sa snehu gc ]]>
</Title>

I also checked the checkbox Use layer names from the map document.
I don't know what else to do to prevent this from happening...


texts: Creating a Grid For Raster When Measuring Land Cover In Accuracy Assessments

texts:
I am assessing the accuracy for the 2011 National Land Cover Map (NLCD).  Each pixel for the NLCD measures 30m x 30m.  For my reference data, I'm using two sources:

1)  Google Earth (~20 m resolution)
2)  NAIP hi res imagery (1-2 m resolution)

For each randomly sampled point in my hi res imagery, I want to create a grid that lays over it and simulates the area of the pixel it represents on the NLCD (30m x 30 m).  This will allow me to label the land cover based on the rules for the NLCD classification scheme (if more than 30% of the 30 m x 30 m pixel is covered in leaf foliage, it is forest).

Does anyone know of a tool for creating a grid like this in Arc that works on raster?  The create fishnet tool is for vector.


texts: Representing clustered points above 80,000 in map using ArcGIS 10.2.1 in such a way that their labels and points themselves are clearly visible?

texts:
I have places of Pakistan which are 85949 in number, I want to display them on a static map in such a way that points and their labels are readable and they don't overlap or map doesn't look messy. I have tried Disperse tool, Aggregate tool and have also tried the following method:

Firstly created buffer on points by the radius of 5 km. Dissolve Type was ALL
Exploded the buffers with Multipart to Singlepart
Used a Spatial Join through Join Data to join the original points to the exploded buffers, giving you buffers with a point count.

I thought this would create polygons with count of places which it did but, around 8000 points came in one polygon which is again difficult to represent with labels in a map.
I am attaching picture as reference.


texts: Django: how to set all of the model's fields using LayerMapping?

texts:
I have a GeoDjango model like below (in "data" app):
from django.contrib.gis.db import models

class FirstModel(models.Model):
    data_id = models.IntegerField(blank=True, null=True)
    model_base_date = models.DateTimeField(blank=True, null=True)
    mpoly = models.MultiPolygonField(blank=True, null=True)

    class Meta:
        db_table = "first_model"

Here is how I used this model in LayerMapping class to import a test shapefile into my PostGIS database:
from django.contrib.gis.utils import LayerMapping
from data.models import FirstModel
from datetime import datetime

mapping = {
    'mpoly' : 'MULTIPOLYGON',

}

world_shp = 'po/my_shapefile.shp'

lm = LayerMapping(FirstModel, world_shp, mapping, transform=False, )
lm.save(strict=True, )

After running the code above, I look inside the table in the database, and all of the model_base_date and data_id fields are empty. Here is a picture from the table:
How can I add the model_base_date and data_id when importing the shapefile into the table?


texts: Finding spatial coordinates in geo-tagged image

texts:
I have a geo-tagged jpg. I know image resolution, GPS location and height, even some angles.

How can I calculate lat/lon for pixel, bounding box, etc from this data?

Can someone give me a hint, maybe formula for this?


texts: Faster alternative to `Intersection` in QGIS to just retain overlapping features without cutting

texts:
I want to keep only the features from my data layer (polygon A) within my area of interest (AOI). If I compare two methods to perform that, the processing time is very different.

Intersection tool

I find the function Vector--> Geoprocessing Tools-->Intersection can be very slow (sys. time. 3981.78 secondes or 660 times longer). It trims each polygon from the data layer (polygon A) by the limits/edges of the second polygon (AOI) and adds the information from the attribute table of the other layer. For my purposes, this is not necessary.

Select by location... method

In contrast, if I used Vector--> Research Tools-->Select by Location..., it can run very quickly (sys. time. 6 sec.) + selecting the layer and exporting only selected features. This method selects the polygons that overlap (or touch, depending on the settings) without cutting out the portions of the polygons that do not intersect.

Is there another tool that will perform the algorithm efficiently with just one step? I don't need to trim the edges of each polygon in the layer I want to cut, just keep those that over within the area of interest (I think that is the difference in processing time, given that that is the difference in the result between the two methods.


texts:
I wanted the following Code to show green icons when the value is < 6.0 and nothing shown when it is >= 6.0 But it is returning Red icons for values >= 6.0 as in screenshot attached. What should be the Else statement for this?

import folium
import pandas

data=pandas.read_csv("EQ_CSV.TXT")
lat = list(data["LAT"])
lon = list(data["LON"])
elev = list(data["MAG"])

fga = folium.FeatureGroup(name="Magnitude below 6.0")
def color_producer(xyz):
    if xyz < 6.0:
        return 'green'
    else:
        return False

map = folium.Map(location=[24.01, 92.01], zoom_start=6, tiles="Mapbox Bright")

for lt, ln, el in zip(lat, lon, elev):
    fga.add_child(folium.Marker(location=[lt, ln], popup= "Magnitude of "+str(el),
    icon=folium.Icon(color=color_producer(el))))


texts:
I'm having trouble exporting an SVG(or a PDF for that matter) from the QGIS print composer. I have a relatively simple map that should, when exported as an SVG, be maybe around 1MB in size, tops. But QGIS is exporting a 450MB file that's way too large for me to do anything with. I'm assuming it's somehow including a whole lot of data that isn't on the page but I can't verify because the file is too large for my computer to open in Inkscape.

I've exported good, small SVGs in the past where everything that didn't show on the print composer page was simply not included. Any ideas what's going wrong this time? Thanks!


texts:
I made this webpage with JavaScript and CartoDB.

The layer selector works like I expected it to but I noticed that the legend are stacking on top of each other and not being removed when the layer is changed. I tried a few different things but to no avail.

A few posts on the internet also relate to a possible bug in CartoDB that may be causing this. Anyone have experience with this?

If you open the map switch layers several times then examine the DOM you will see that the map.remove(); option is not totally clearing all elements of the layers for the previous selection.

Example:

Any suggestions?


texts: Published ESRI layer scale dependency problems

texts:
I have published a layer with 3 scale dependency settings so regardless of the scale one of the three layers should be visible at any one time.
When I view the layer using the "View In: ArcGIS JavaScript" link (from the ArcGIS REST Services Directory page) the layer works correctly, I zoom in and out and the features are rendered correctly.
When I then switch to my development map and include the layer in the map layers list and load the map up, the layer appears on the zoomed out scale fine (global view).  When I then zoom in past the threshold my layer completely disappears.  If I inspect the JavaScript I can navigate to map.layers.item[x] and the visible property is set to true but the layer simply does not render on the map.
Does anyone have any idea what is going on?
Layer caching is disabled.
The only difference I can see between the ArcGIS JavaScript view and my development map is the layer is a FeatureLayer in my dev configuration.


texts: Leaflet routing machine Origine to destination

texts:
I want to do routing from a point to another the points are listed in legend control .the data retrieved from listbox does not work with routine machine may be because of type compatibility.

  var mymap=L.map('map',{attributionControl: false,center : [33.37,6.84],zoom: 12,minZoom:7,maxZoom: 17,zoomControl:false});
    var GRaster = L.tileLayer('http://{s}.google.com/vt/lyrs=s&x={x}&y={y}&z={z}',{subdomains:['mt0','mt1','mt2','mt3']}).addTo(mymap);
    function Itineraire(){
        var x = document.getElementById("Origin").selectedIndex;
        var y = document.getElementById("Origin").options;
        var x1 = document.getElementById("Destination").selectedIndex;
        var y1 = document.getElementById("Destination").options;

    //clear route
    Routing.getPlan().setWaypoints([]);
    //retrieve data
    var point1=y[x].value;
    var point2=y1[x1].value;

     alert('Routing :'+"pt1: "+y[x].value+" pt2: "+y1[x1].value);

    Routing.getPlan().setWaypoints([
    L.latLng(point1),
    L.latLng(point2)
    ]);
    }
    var Routing =L.Routing.control({
    waypoints: [
    L.latLng(33.357464,6.864267),
    L.latLng(33.374304,6.871693)
    ],
    routeWhileDragging: true
    });
    Routing.addTo(mymap);
    Routing.getPlan().setWaypoints([]);

    var legendbox = L.control({
    position: 'topleft'
    });

    legendbox.onAdd = function (e) {
    var div = L.DomUtil.create('div', 'info legend');
    div.innerHTML = '<select id=Origin><option value=33.374304,6.871693>Zone1</option><option value=33.367979,6.837744>Zone2</option><option value=33.372684,6.839808>Zone3</option></select><select id=Destination><option value=33.374304,6.871693>Zone1</option><option value= 33.367979,6.837744>Zone2</option></select>';
    div.innerHTML += '</br><button type = button onclick = Itineraire()> Find route </button>';
    div.firstChild.onmousedown = div.firstChild.ondblclick = L.DomEvent.stopPropagation;
    return div;
    };
    legendbox.addTo(mymap);

My running code on plunker


texts: Assigning points to their respective Reaches using Stream Gradient?

texts:
This is a continuation of the question I initially posed here.

Also, I've attached a shapefile on my Drop Box account here.

I'm trying to break a stream network up into 100m reaches, assigning elevation values pulled from a DEM, to the start and end of each segment, the goal being to determine gradient for each segment.

Process-wise, I'm not approaching this correctly. So far, I haven't been able to those z-values to the start and end nodes. What I've been doing is to segment the stream network, and using XTools, convert Features to Points, again, at 100m intervals. They do line up with the segments, as you can see below. From there, I can easily use XTools to pull raster values, and assign them to the points.

What I'm not able to do, is to assign those points back into the reaches. I can get one of them, sure, but not both. Thus, I can't calculate gradient.

Again, any help would be MUCH appreciated, particularly any sort of working example.


texts: Are there a list of locations or features with known dimensions?

texts:
I'm plotting lines and objects on a map and want to verify scales with real, physical objects on the map. Is there a list of physical features that are known distances?
Perhaps like the runway at Laguardia is 7000 feet long, exactly from this point to this point.
Or, this particular road is exactly due North, and this other is East-West.
Where is the calibrated list of map features bridges, roads, runways or road intersections that have known dimensions?


texts: Running Near 3D + Alter Field sequentially

texts:
I am asking for help with generating an automated process to do what I am tediously doing manually at the moment. I want to know the 3D distance between buildings and a tunnel. I have several building points and a set of tunnel points (one point per meter inside the tunnel).

This is how I have been doing it so far: I am running Near 3D tool to find distances between each building point and each tunnel point (I only have the basic license so I cannot use Generate Near Table tool). In Near 3D the tunnel point set is the input feature and one building point is the near feature. Afterwards I run Alter Field tool to change the name of the 3D distance field in the attribute table (NEAR_DIST3) to the address (and filename) of the building point:

I then do this with every single building point, in model builder, so the attribute table of the tunnel points gets several fields named after the given building point.

I have attempted writing a script to solve this myself, but I am not very experienced with python so it is just from the model builder:

# Near 3D
arcpy.Near3D_3d(inpoints_N2, "'C:\\GIS\\distance.gdb\\Streetname123Apt456'", "2000 Meters", "NO_LOCATION", "NO_ANGLE", "NO_DELTA")
# Alter Field
arcpy.AlterField_management(inpoints_N2__9_, "NEAR_DIST3", "Streetname123Apt456", "", "DOUBLE", "8", "NULLABLE", "false")
# Near 3D (2)
arcpy.Near3D_3d(inpoints_N2__2_, "'C:\\GIS\\distance.gdb\\Streetname123Apt457'", "2000 Meters", "NO_LOCATION", "NO_ANGLE", "NO_DELTA")
# Alter Field (2)
arcpy.AlterField_management(inpoints_N2__6_, "NEAR_DIST3", "Streetname123Apt457", "", "DOUBLE", "8", "NULLABLE", "false")

Writing this manually also seems unnecessarily time consuming. Does anyone have a suggestion how I can setup a model in model builder or rewrite the script to avoid doing this point by point?

I am using ArcGIS 10.3 with basic license and 3D Analyst and Spatial Analyst extensions.


texts: Arcpy Error 010423 when using Zonal Statistics to Table in a tool

texts:
I want to include Zonal Statistics to Table to my tool. But if I do I get the following error:

Traceback (most recent call last):
  File "C:\Users\xyz\Desktop\test.py", line 16, in <module>
    outZSaT = ZonalStatisticsAsTable(inshape, "field_id", ndvi, outTable, "DATA", "MEAN")
  File "c:\program files (x86)\arcgis\desktop10.1\arcpy\arcpy\sa\Functions.py", line 5836, in ZonalStatisticsAsTable
    statistics_type)
  File "c:\program files (x86)\arcgis\desktop10.1\arcpy\arcpy\sa\Utils.py", line 47, in swapper
    result = wrapper(*args, **kwargs)
  File "c:\program files (x86)\arcgis\desktop10.1\arcpy\arcpy\sa\Functions.py", line 5828, in wrapper
    statistics_type)
  File "c:\program files (x86)\arcgis\desktop10.1\arcpy\arcpy\geoprocessing\_base.py", line 498, in <lambda>
    return lambda *args: val(*gp_fixargs(args, True))
ExecuteError: ERROR 010423: C:\Users\xyz\AppData\Local\Temp\arc5F5D\t_t646 does not have valid statistics as required by the operation.
Failed to execute (ZonalStatisticsAsTable).

Failed to execute (test).
Failed at Tue Nov 07 08:31:06 2017 (Elapsed Time: 0.00 seconds)

Here is the code:

import os
import os.path
import sys
import datetime
import arcpy
from arcpy import env
from arcpy.sa import *

inshape = arcpy.GetParameterAsText(0)
ndvi = arcpy.GetParameterAsText(1)
out_location = arcpy.GetParameterAsText(2)

outTable = out_location + "\\outtable"
outZSaT = ZonalStatisticsAsTable(inshape, "field_id", ndvi, outTable, "DATA", "MEAN")

inshape: Feature layer, ndvi: Raster layer and out_location: Folder

If I replace the arcpy.GetParameterAsText with the actual path to the files and location and run in the python window it works fine. So where that t_t646 comes from?


texts: Converting geographical projection of extremely large file using GDAL?

texts:
My Task
I need to convert 9 x 14Gb GeoTiff files to WGS84 geo projection.  The rough size in pixels is about 200k x 100k.

The data-type is byte and it is not a subject to interpolation, it is some logical attribute.  For example, if A = 20 and B = 30 does not mean that overall region is 25 - it has completely different meaning.

What I did

Started gdalwarp on Windows (takes lots of time, fails on some machines, never finishes on other)
Started gdalwarp on Ubuntu (Intel(R) Xeon(R) CPU E5-2660 0 @ 2.20GHz, 56Gb RAM, takes ~3 days for 50% and stops). Takes 100% of one core, low RAM and low HDD io, and also puts Potential thrashing on Band1 message in debug-log.
Modified arguments, e.g. CACHE_MAX, wa were set to some really big values.
The command I use now:
gdalwarp -t_srs EPSG:4326 --config GDAL_CACHEMAX 1000 -wm 1000 --debug on $1 $2

Question
What are my other options?

I have successfully downloaded and compiled GDAL Beta 2. I have also placed entire source tif file in tmpfs (RAM). My current cmd is

gdalwarp -t_srs EPSG:4326 --config GDAL_CACHEMAX 99000 -wm 2000 -multi -co TILED=yes -wo NUM_THREADS=ALL_CPUS -wo STREAMABLE_OUTPUT=TRUE -oo GTIFF_DIRECT_IO=YES -oo GTIFF_VIRTUAL_MEM_IO=YES ram/2014/2014/2014_30m_cdls2.tif /mnt/test2.tif > test2.log 2>&1 &

It does ~40% pretty quick (5 mins). Other 5 minutes, that I waited, it did not provide any output. From that point, gdal stops using more than one core.


texts: Calculating distance in polygons using ArcGIS Desktop?

texts:
I want to show on map what's the distance between polygon A and B, but the result should be specified in polygons.

Here's the result similar to my idea:

I am working on ArcGIS Desktop.

Which tools should I use to get this result?


texts: Creating custom Operation for encapsulating some map algebra using GeoTools?

texts:
In GeoTools, I need to perform a one-to-one pixel computation between a bunch of GridCoverage2D. In short: some map algebra.
As I understood I can do things like summing values of coverages using:

Coverage summCoverage = Operations.DEFAULT.add(coverage1, coverage2);

I would like to create an Operation involving multiple coverages.. ..but I didn't find any hint about how to proceed.


texts: Are there any OpenSource ArcObjects Projects?

texts:
One of the major benefits of developing GIS applications based on opensource libraries, is that many projects spring up on top of them which are also opensource.

Well designed opensource projects are great for seeing:

which coding patterns should be used
what coding styles to follow
how to structure larger projects
what unit tests to include
autodocumentation tools and procedures
sourcecontrol workflows

The ArcObjects Code Gallery seems the best place to look, but most examples I've looked at are very small pieces of functionality or single tools, rather than full applications.

The SharpGIS project is a nice .NET example to follow, but are there any ArcObjects specific projects?

Updated with Examples:

http://eagl.codeplex.com/

  Easy ArcGIS Library is a set of C# .net classes that wrap the common
  functionality of ArcObjects, that help ArcGIS developers do a lot of
  common functions in less time and direct way.

And from the same developer:

http://arcgisslcontrols.codeplex.com/

  ArcGISSLControls is a set of ArcGIS Silverlight controls that is built
  on top of the ArcGIS Silverlight SDK, it is developed in C#
  programming language.


texts: Working with shapefiles in MapServer: adding color, smooth edges

texts:
I am successfully serving Koeppen-Geiger climate shift data from a MapServer instance. Here is an example WMS request to show the current state.

Unfortunately this looks nothing like the ArcGIS Koeppen-Geiger map. For comparison here is a full screenshot of the layer I am serving from the raw data:

I downloaded the map from here. There are a lot of files in the shapefile zip, and I'm not sure if any of them can help me get the desired output. Here is a screenshot of the files:

Additionally, when I load the shapefile into QGIS by itself it looks the same as delivered from MapServer.

Can anyone explain how I improve the map's aesthetics?


texts: QGIS - using fill ring and snapping to segments not working

texts:
When using the fill ring tool and having to run along the edge of the polygon the fill ring tool works perfectly when snapping to a vertex, but when I try and snap to a segment, which appears to work fine until you finish the fill ring and you get an error saying the ring was not contained in the feature.
What other option do I have to cut out an area of the polygon along the boundary of the polygon while still being offered the attributes table to name that new feature?


texts: Calculate weight using distance from point for Weighted Reductions in Earth engine

texts:
I want to calculate the total yearly NO2 emission from the power plants by using Sentinel-5p NO2 data. First I defined my aoi and sum all 365 images into one image using ee.Reducer.sum(). After that reduceRegion (with sum reducer) to find the total yearly sum of no2 of my aoi. But this picks a lot of noise from other emission sources like transport.

For avoiding those, I have planned to use Weighted Reductions. I want to use weight in reduceRegion in a way where pixels away from power plants will receive less weight and higher weight to nearest pixels to power plants.

From the documentation, we can see how to calculate weights for clouds ( the more cloud, the less weight each pixel will get ). But how do I calculate weights in my case where each pixel will receive weights according to distances to nearest power plants.


texts: Double labeling line problem using CSS for OSM data from Geoserver

texts:
I had created a CSS for postgres planet_osm_line source on geoserver, for roads displaing
in client, using WMS.
CSS seted like:

/*.....*/

 /* @title Residential */
 [@scale < 20000] [highway = 'residential'] {
    stroke: #e2e1e1, #fefefe;
    label: [name];
    z-index: 0, 1;
    stroke-linecap: round;
    font-family: DejaVu Sans;
    font-weight: bold;
    font-size: 10;
   label-follow-line: true;
     label-max-angle-delta: 90;
     label-max-displacement: 250;
     label-repeat: 300;
  }

 /* @title Tertiary */
 [@scale < 20000] [highway = 'tertiary'] {
    stroke: #e5e5a8, #f8f8ba;
    stroke-width: 10px, 8px;
    label: [name];
    color: #dadada;
    z-index: 2, 3;
    stroke-linecap: round;
    font-family: DejaVu Sans;
    font-weight: bold;
    font-size: 10;
   label-follow-line: true;
     label-max-angle-delta: 90;
     label-max-displacement: 100;
     label-repeat: 100;
  }
 /*.....*/

But a unnormal artefact was detected: some line labels displays twice after n Zoom level:

red line shows double label, green normal label.
How to fix it using css or geoserver or OpenLayers configs?


texts: Is it possible to have a 1:n join which takes the min, max, or average value?

texts:
I need to join the attributes of a shapefile with the one contained in a .dbf (or .txt) file. I currently use the "join" in the properties menu, but I have a problem. In the shapefile I have one identity code for each building; in the .dbf (or .txt) I have several noise levels related to the several façades of the same building. So when I use "Join" I miss some data because the results is only one value that seems to be taken from the last point considered in the merge of the properties. Is it possible to do this operation in a different way? For instance taking the minimum, the maximum and/or the average value of the noise level...
Thanks a lot in advance for your useful suggestions!


texts: Postgres returning incorrect results with distance join queries

texts:
I created a table containing locations along with their spatial coordinates expressed in (latitude, longitude (in degrees)) in postgres. Commands which I used for the same are:

create table spatialTest(name character varying(100), the_geo geography);

\copy spatialTest(name,the_geo) FROM 'testSpatial.csv' DELIMITERS E'\t' CSV HEADER;

I created a table containing locations along with their spatial coordinates expressed in (latitude, longitude (in degrees)) in postgres. Commands which I used for the same are:

create table spatialTest(name character varying(100), the_geo geography);

\copy spatialTest(name,the_geo) FROM 'testSpatial.csv' DELIMITERS E'\t' CSV HEADER;

testSpatial.csv contains the following values:

A   SRID=4326;POLYGON((0.178773,-127.681841|0.178711,-127.681962|0.179125,-127.682083|0.179176,-127.682006|0.179153,-127.681986|0.179143,-127.681962|0.179147,-127.681935|0.179166,-127.681913|0.179195,-127.681897|0.179244,-127.681886|0.179284,-127.681887|0.179336,-127.681904|0.179464,-127.681757|0.179489,-127.681736|0.179429,-127.681680|0.179370,-127.681516|0.179221,-127.681331|0.179184,-127.681185|0.179051,-127.681264|0.178822,-127.681499|0.178761,-127.681698|0.178796,-127.681703|0.178839,-127.681721|0.178857,-127.681736|0.178861,-127.681740|0.178871,-127.681756|0.178873,-127.681782|0.178859,-127.681809|0.178843,-127.681825|0.178812,-127.681839|0.178773,-127.681841))
B   SRID=4326;POINT(0.628912,-127.700922)

Now I want to find all spatial locations which are within a distance of 50 km of each other. For doing so I used the following command:

select s1.name, s2.name  from spatialTest s1,
     spatialTest s2  where ST_DWithin(s1.the_geo, s2.the_geo, 50000);

However, to my surprise I found that although A and B are separated from each other by a distance greater than 50 km (50.0995 km to be precise. Found using Chris Veness's geodesy formulae (Calculate distance between a point and a line segment in latitude and longitude)), yet they are returned by postgres as results.

Where am I going wrong?

I am using PostgreSQL 9.6devel and Postgis version which I am using is: POSTGIS="2.2.1 r14555"


texts: How to display GeoJSON file over OSM, with Google satellite as a base layer (at correct location) in open layer?

texts:
How can you display a GeoJSON file over OSM, with Google satellite as a base layer (at the correct location) in OpenLayers?

var map = new OpenLayers.Map('map',{
            projection: "EPSG:900913",
            displayProjection: "EPSG:4326",
    });
    layer = new OpenLayers.Layer.WMS("districts_2011census",
                    "http://localhost:8081/geoserver/satya/wms?",
                { layers: 'districts_2011census',
                  format: 'image/png'
                },
                { layers: 'basic' }
            );
var gmap = new OpenLayers.Layer.Google("Google Streets",
      { sphericalMercator: true,numZoomLevels: 20} );
      geojsonayer1 = new OpenLayers.Layer.Vector("GeoJSON23", {
          Projection: new OpenLayers.Projection("EPSG:4326"),
                strategies: [new OpenLayers.Strategy.Fixed()],
                protocol: new OpenLayers.Protocol.HTTP({
                url: "../geojson/result_routing.geojson",
                format: new OpenLayers.Format.GeoJSON()
                }),
                strategies: [new OpenLayers.Strategy.Refresh({force:true}),
 new OpenLayers.Strategy.BBOX({ratio:2, resFactor: 3})],
                style: {
                    graphicName:'cross',
                    pointRadius: 5,
                  'strokeWidth':4,
                  'fillColor':'#ff0000',
                  'strokeColor':'#ff0000'
                }
                });

        var bounds = new OpenLayers.Bounds(68,8, 95, 34);
        map.addLayers([layer,gmap,geojsonayer1]);

        map.addControl(new OpenLayers.Control.LayerSwitcher());
        map.addControl(new OpenLayers.Control.MousePosition());
        map.addControl(new OpenLayers.Control.ScaleLine());
        map.addControl(new OpenLayers.Control.KeyboardDefaults());
        map.zoomToExtent(bounds);