texts: Java GeoTools API: getInteriorPoint does not work on an intersection between Polygon and MultiLineString

texts:
I'm working on a Java project using GeoTools API that want to get the interior point getInteriorPoint(), from the geometry coming from an intersection between a Polygon and a MultiLineString. Please note that the same approach works perfectly between a Polygon and MultiPolygon, so between objects with the same number of dimensions, I guess.

Let me explain better with a snippet of code (java GeoTools API version 15.1):

if(feature.getDefaultGeometry() instanceof MultiLineString){

   MultiLineString multiLineString = ((MultiLineString) feature.getDefaultGeometry());

   //Interior point on intersection
   Coordinate intersectionInteriorPoint = this.getInteriorPointForIntersectedArea(polygon, multiLineString, sourceCRS);

}

where the getInteriorPointForIntersectedArea(..) is as this:

   /**
     * Compute the interior point for the overlapping area
     * between the search geometry and an intersected geometry
     * @param one
     * @param two
     * @return the coordinate in the target output reference CRS
     */
    private Coordinate getInteriorPointForIntersectedArea(Polygon one,
                                                          MultiLineString two,
                                                          CoordinateReferenceSystem sourceCRS,
                                                          Coordinate waypointCoordinate){

        com.vividsolutions.jts.geom.Geometry lineString = (com.vividsolutions.jts.geom.Geometry) one.intersection(two);

        System.out.println("Geometry type: " + lineString.getGeometryType());

        Coordinate interiorPointCoordinate = lineString.getInteriorPoint().getCoordinate();

        Coordinate interiorPointCoordinateChanged = null;
        try {
            interiorPointCoordinateChanged = CoordinateReferenceSystemTools.FromToCRS(sourceCRS,
                                                                                org.geotools.referencing.crs.DefaultGeographicCRS.WGS84,
                                                                                interiorPointCoordinate,
                                                                                true);
        } catch (FactoryException |
                 TransformException |
                 NullPointerException |
                 IllegalArgumentException e) {

                 interiorPointCoordinateChanged = new Coordinate(0,0);
                 System.err.println("Unable to extract feature interior point: ");
                 e.printStackTrace();
        }

        System.out.println("MultilineString interior point lat: " + interiorPointCoordinateChanged.y + " lon: " + interiorPointCoordinateChanged.x);
        return interiorPointCoordinateChanged;
    }

I apply this on multiple features and sometimes the same code return the correct interior point. Some other times, the method below fails with exception:

*java.lang.IllegalArgumentException: Argument "source" should not be null.
    at org.geotools.geometry.jts.JTS.ensureNonNull(JTS.java:135)
    at org.geotools.geometry.jts.JTS.transform(JTS.java:460)*

Looking into JTS at line 135, I find that name = "source" coming into this method, as null value in object when the exception happens:

private static void ensureNonNull(final String name, final Object object)
            throws IllegalArgumentException {
        if (object == null) {
            throw new IllegalArgumentException(Errors.format(ErrorKeys.NULL_ARGUMENT_$1, name));
        }
    }

I have discovered so, that the exception occurs when the Point coming out from the call:

lineString.getInteriorPoint()

is a POINT EMPTY

Does anyone of you knows why sometimes, getInteriorPoint() returns a POINT EMPTY ? This happens to me only when I try to call getInteriorPoint() on a geometry coming from the interesection between w Polygon and a MultiLineString.


texts: Issue using pgr_dijkstra from Postgres pgRouting

texts:
I am using Postgres with Postgis extension to retrieve streets (only for pedestrian use) and several "natural" points. I need to make updates every minute and i am working on a small bounding box.
For rendering I am using leaflet that makes request to a java servlet and all is working well...Postgres is populated with "osm2pgsql".

I would implement some kind of routing service, unfortunately i am not expert with PLPGSQL dialect  but i found this excellent post in LOUDHUSH BLOG.
To make it work I had to change the name of the last function to "pgr_createTopology":

drop table if exists network;
create table network(gid serial, osm_id bigint, name varchar, the_geom
geometry, source bigint, target bigint, length float);

CREATE OR REPLACE FUNCTION compute_network() RETURNS text as $$
DECLARE
streetRecord record;
wayRecord record;
pointCount integer;
pointIndex integer;
geomFragment record;
BEGIN

-- for each street
FOR streetRecord in select way, osm_id, name from planet_osm_line where
highway is not null LOOP

-- for each street in the region of interest
 SELECT * from planet_osm_ways where id = streetRecord.osm_id into wayRecord;
 FOR pointIndex in array_lower(wayRecord.nodes, 1)..array_upper(wayRecord.nodes,1)-1 LOOP

RAISE NOTICE 'Inserting name % source %, target %', streetRecord.name, wayRecord.nodes[pointIndex], wayRecord.nodes[pointIndex+1];
  select st_makeline(st_pointn(streetRecord.way, pointIndex), st_pointn(streetRecord.way, pointIndex+1)) as way into geomFragment;
  insert into network(osm_id, name, the_geom, source, target, length) values(streetRecord.osm_id, streetRecord.name, geomFragment.way,
  wayRecord.nodes[pointIndex], wayRecord.nodes[pointIndex+1], st_length(geomFragment.way));
 END LOOP;
END LOOP;
return 'Done';
END;
$$ LANGUAGE 'plpgsql';

select * from compute_network();
select pgr_createTopology('network', 1, 'the_geom', 'gid');

"network" table is populated in the right manner..i think,
before to run Shortest Path Dijkstra i changed manually "source" and "target" column to integer type but the query below return an empty table

CREATE TABLE route AS SELECT seq, id1 AS node, id2 AS edge, cost, b.the_geom FROM pgr_dijkstra('
            SELECT gid AS id,
                     source::integer,
                     target::integer,
                     length::double precision AS cost
                    FROM network',
            41, 50, false, false) a LEFT JOIN ways b ON (a.id2 = b.gid);

can you tell me what am I doing wrong?


texts: Subsetting a curvilinear netCDF file (ROMS model output) using a lon/lat bounding box.

texts:
I am attempting to get salinity data from one of the ocean models in our region and have run into a roadblock. This is ROMS model output in NetCDF form, but the grid is curvilinear (lon and lat arrays are 2D).

I want to get a subset of the data for the bounding box for the area I am interested in, however I am unsure how to go about doing this. Given a lat/lon/depth for a bounding box of interest, what are the grid cells and vice versa.

I didn't know if using PyDap would make it any easier to get the data of interest, or if there was another toolset I should look at. I saw Rob Hetland's Octant, but it looks like it works on NetCDF files and I want to avoid pulling the files down if possible since they are fairly substantial.

Thanks

Dan


texts: Recommended specs for a single-country OSM tile server

texts: What specs would you recommend for an OSM tile server dedicated to serving tiles for a single, small country (e.g., Portugal only)? Total server load is estimated to be no more than two/three users accessing tiles at the same time.


texts: Iterate a clip over multiple feature classes using Model Builder or ArcPy

texts:
I have multiple feature classes (buildings, airport, land, roads etc) that I would like to clip using a Local Authority boundary. I am trying to create a process in model builder or Arcpy that will iterate the clip operation through the feature classes and create separate output (clipped building, airport, land, roads etc) for each feature class.

Please can anyone advice me on how to automate this than having to set separate clip operations for each feature class?

I've tried using the Iterate feature class in model builder but it only selects and clips the first feature class (Airports) not all feature classes.


texts: Why GRASS 7.0.0 can't create a table while importing shapefiles from QGIS?

texts: I use GRASS7.0.0Beta1 and QGIS 2.4.0 on Windows8. With v.in.ogr, when I import a shapefile, I get the error message "Unable to create table". When I put the option of not importing the table, the vector file is imported. It used to work with GRASS6...


texts: Debugging ArcPy ERROR 000358: Invalid expression?

texts:
I've tried several iterations of this SQL statement in a stand alone script:

    QF5 = "Week"
    QF6 = "Wk"
    RV3 = ""
    WC4 = """" {} <> {} AND {} <> '{}' """.format(QF5,QF6,QF6, RV3)

I'm aiming for an expression that would read int he field calculator:

Select From(*) Table Where:
Week <> Wk AND Wk <> ''

This is the error I'm getting:

ERROR 000358: Invalid expression
Failed to execute (SelectLayerByAttribute)


texts:
I have a db with two schemes "public" and "versioned" (the latter generated with the fastversion plugin).
In the "public" scheme I have shp with related styles_table, in "versioned" I have other shp that always look to the styles_table in "public" schema. Everything works.
The problem is that this plugin generates versions of my shp through the views (with automatic names) and these views do not recognize the basic shp style (table_style in a "pubblic" schema).

There is a way to give at the views the styles present in the "public" schema automatically?


texts:
I'm currently thinking about ways to develop scripts with openCV library and use with QGIS.
I've installed python directly from Anaconda distribution and installed new packages from it (never had a problem with anaconda, and in fact I love it).
I wonder that the python from QGIS is installed directly from QGIS, and not from other sources in my computer.
Is there a way that I can link QGIS Python to Anacondas packages or Do I need to install 'new' packages for QGIS?


texts: Reclass in ModelBuilder - Stopping at Random Raster

texts:
I've created 5,000 viewshed rasters for ~5,000 sites in my study area. Using ModelBuilder, I am trying to reclass these rasters so that all "0" values are now "NoData".

When I iterate across all rasters in my folder, it seems to work well until I get to the 4,001th raster. Then the operation stops and I receive an error message. This has happened three times.

It doesn't seem to be a matter associated with naming, as the operation works for the first 4,000 rasters. Does anyone have any idea why this may be happening? I don't have the error message recorded, unfortunately, but I'm wondering if anyone has come up on similar issues.


texts:
We are using ST_Geometry for Oracle 11g at the moment.

And we need to query features (points,lines,polygons in a view) inside a given extent(either a circle or a polygon or a rectangle), and we use the st_envintersects to do the spatial calculation.

However since some of the features are poygons which make the query rather slowly even we create indexes for the geometry column of each feature.

For example there is a view feature LView with columns:

objectid
lon(
lat
geom(the geometry, point or polygon or polyline)
...

And we use the following sql to do the query:

select * from Landuse l where sde.st_envintersects(l.geom,sde.st_geometry('POLYGON((....))',3)) = 1;

It will take almost 2-3 minutes to do the querying and there is only 2000 records in the view.

Then when we re-create the LView and remove the features of polyon and polylines whose geometry is too complex.

The same query will take only a few seconds.

It seems that the spatial calculation for complex the polygon will slow down the query. So I want to treat the polygons as points(by the lon lat column) during the query.

Something like this:

select * from Landuse l where sde.st_envintersects(sde.st_geometry('Point(l.lon,l.lat)',3),sde.st_geometry('POLYGON((....))',3)) = 1;

Is this possible

Update:

I am using the Arcgis 10.0.

A complex polygon(multiple polygon) may contain 2000+ points.

Index: CREATE INDEX spatialidx_landuse_polygon ON landuse_polygon(geom) INDEXTYPE IS sde.st_spatial_index PARAMETERS('st_grids=1,0,0 st_srid=3')

View:

create or replace view LView as
  select objectid,geom,'landuse_polygon' as table_name, lon, lat from landuse_polygon
union all
    select objectid,geom,'school_point' as table_name, lon, lat from school_point
    ....................

The view contains almost 20 + tables.

I have not register the view with sde. ( I did not know this at all).

select count(*) from LView will give result of 2670.

Query plan: I have checked the query and the query for the table of polygons take too long.


texts: Distance in km between coordinates, Google Earth Engine

texts:
I have coordinates for a bounding box:
1: [-109.41397027896355,22.86994450876419]
2: [-109.41397027896355,28.00000232204317]
3: [-115.75044209344128,28.00000232204317]
4: [-115.75044209344128,22.86994450876419]

And I want to find the distance between coordinates in order to get the height and width of the bounding box in kilometers. Are there any specific functions that do this in Google Earth Engine?


texts: Extract geojson 'properties' dynamically with mapbox gl

texts:
I have an interactive map based on MapBox WebGL. I have an 'on click' function which is returning all information on data layers at that location (lovely).
For each object I am trying to pull out just 'layer.id' (layer title) and 'layer.properties'(held as an object) and display that only.

I can't change the geojson so i need to be able to access the properties without knowing what they keys are.

I figured it would just be a case of looping through as an object array but i can't get it right.

Currently this partially works in the console;

map.on('click', function (e) {
        map.featuresAt(e.point, {radius: 5, layer: lyrlist}, function (err, features) {
        if (err) throw err;
        //console.log(features)

        var keys = Object.keys(features);
        var val = "";

        for (var i = 0; i < keys.length; i++) {
            val = features[keys[i]];

    //console.log(val.layer.id,val.properties); shows layer title and property object in console.

        }
        //document.getElementById('features').innerHTML = JSON.stringify(features, null, 2); Shows all info in panel.
    });
});

link to demo site: http://www.activenaturalist.org.uk/lcren/imap/mb.html


texts: Color gradient to show direction in polygons (filling, color ramp, gradients)

texts:
My data is 2 point-layers: 1 layer origins and 1 layer destinations. I connected these pairs of origins and destinations: firts with lines. As I know that in qgis, there is still no gradient fill for lines, I did a very small buffer on mi lines. And now, I have polygons, wich represent my flows, and I would like to represent them with a gradient fill wich show the direction of each flow (it means: from the origin point to the destination one). I have an idea but I don't achieve to go further.

Do you know how to use options of the gradients filling?

I have the idea to use the two Reference Points to represent the orientation of the flow (from origin to destination).

I would like to use 4 table: - abscissa of my origin points (as x of reference point 1) - ordinate of my origin points (as y of reference point 1) - abscissa of my destination points (as x of reference point 2) - ordinate of my destination points (as y of reference point 1) But it doesn't work.


texts: gdal warp reproject creates artifact in raster

texts:
I’m trying to re-project a tif-file using gdal warp:
fn = "\dtm.tif"
out = "\dtm_25832.tif"
gdal.Warp(out, fn, dstSRS="EPSG:25832")

The output file is re-projected, but it creates artifacts in the output file. Is there a way to avoid creating artifacts? Artifacts are more visible on hillshade-dem.

from osgeo import gdal
import os
fn = 'DEM.tif'

# open dataset
rds = gdal.Open(fn)
ds = gdal.Open(fn)
ds.GetMetadata()
{'AREA_OR_POINT': 'Area',
'TIFFTAG_RESOLUTIONUNIT': '2 (pixels/inch)',
'TIFFTAG_XRESOLUTION': '72',
'TIFFTAG_YRESOLUTION': '72'}
rds.GetDescription()
'C:\\Python\\Test\\DEM.tif'
img_width,img_height=rds.RasterXSize,rds.RasterYSize
img_width,img_height
(15010, 30010)
ds.GetGeoTransform()
(-39575.0, 1.0, 0.0, 6726005.0, 0.0, -1.0)
ds.GetProjection()

'PROJCS["ETRS89 / UTM zone 33N",GEOGCS["ETRS89",DATUM["European_Terrestrial_Reference_System_1989",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6258"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4258"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",15],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","25833"]]'


texts: List of available transformations between two coordinates system in epsg using geotools

texts:
How we can have the list of available transformation between two coordinate systems using geotools? e.g. the list provided of epsg.io between EPSG:32066 and EPSG:4326

https://epsg.io/32066

Java snippet using geotools would be very helpful.


texts: Is it possible for QGIS to recognise GeoPackage Schema and/or Related Tables extension implementation?

texts:
I have set up a basic geo package example with the schema extension, specifically for data constraints and enums. When I bring this into QGIS v3.16, I cannot see my extension work implemented anywhere (i.e. enum pick-list not showing when I add a new feature).
Is this possible at all in QGIS yet? Or do I just need to set up a "Value Map" widget in the Attributes Form properties window for now?


texts: Odd behavior in a QGIS plugin: my function is triggered twice

texts:
I have a plugin with a function (triggered by a button) that selects a specific layer by name, then shows its attribute table. When I close the plugin and open it again, the attribute table pops up twice, and repeating this keeps adding another table popup. What's wrong?

Code:

def showTable(self):
    for lyrs in QgsMapLayerRegistry.instance().mapLayers().values():
        if lyrs.name() == "TestFeatures":
            lyrTest = lyrs
    self.iface.showAttributeTable(lyrTest)


texts: Is there a way to view KML natively in ArcMap?

texts:
Answering this question: How do I import/export KML to ArcGIS Desktop?  raised another one.  There are loads of tools to import/export kml/kmz.  Are there tools available for native viewing of KML in ArcMap?

I'm curious if anyone has built a custom featurelayer(s) that reads kml, similar to the way ZigGIS uses PostGIS as a data source.

Or maybe something like inmemoryworkspacefactory, which probably doesn't count as native, but would have some advantages:
http://geobabble.wordpress.com/2007/07/10/consuming-georss-in-arcmap-with-inmemoryworkspacefactory/


texts: How do I know if my ArcGIS 9.3.1 has a GIS server on my computer?


texts: Can't Link grass to remote server postgis topology schema for import or export vector data

texts:
I have installed postgre/postgis on a remote server in which i want add topology data from Grass Gis to postgre topology schema.

Grass gis 7.2.1 is installed on another PC on the network which has access to the postgis database through Qgis 2.8.11.

i was wondering how can i add postgis databse connection to the Grass gis 7.2.1 so i can use that connection for importing and exporting vector layer?

Without a configured connection i get this error while trying to define an output format:

In all videos i searched i see examples with already linked to postgis, and i don't know how to link it to my database.

This is the log in log i get:

"(Mon Aug 21 10:35:49 2017)
db.login --overwrite driver=pg database=postgres user=postgres password=xxxx host=localhost port=5433
WARNING: DB connection  already exists and will be overwritten
The password was stored in file (C:\Users\Administrator\AppData\Roaming\GRASS7\dblogin)
(Mon Aug 21 10:35:49 2017) Command finished (0 sec)"


texts: What is a good, small GPS data logger that stores accuracy information?

texts:
I'm looking for a small GPS data logger for geotagging photos, and I've had trouble finding one that meets my requirements:

Small, key-fob size would be ideal.  It should also store its tracks on on-board storage.

I will be carrying it around on my person for days at a time.
Good battery life (16 hours per charge and/or on-the-go USB charging or replaceable batteries)

I can't be fussing with it or stopping to charge during the day.
Records accuracy information (HDOP, VDOP, PDOP, satellites in view, etc)

I want to be able automatically filter out bad data during post-processing.

I've been able to find plenty of units that meet requirements 1 & 2 (see below), but 3 seems like more of a specialist thing and I have had a hard time finding information about what supports it.

I've already tried two loggers, but they've had issues:

Holux M-241: Runs on Alkaline AAs (only high quality), but you need to change them every few hours or it gets flaky (often resulting in data corruption).  It was apparently hackable with earlier firmware revisions to get accuracy info, but the newer versions (which I have) clear the nonstandard config on startup.

Canmore GP-102+: Nice form factor and very good battery life.  However, some of the tracks it recorded had really inaccurate sections (probably when the signal was poor) and I had to spend a lot of time editing out the nonsense by hand.  Stores tracks in .fit files that are accessed like a flash drive.  GPSBabel does not appear to find accuracy info in them when I convert them to GPX.


texts: How can I get GeoKettle (or any other Java based ETL) to run on MacOS Sierra?

texts:
Has anyone managed to get GeoKettle or Spatial extension for Talend running on MacOS Sierra? If so, how did you do it?

These packages have not been updated since 2015 and Apple has more or less discontinued support for Java in Sierra. I have installed both Oracles and Apples versions of Java, but to no avail.

I am looking for a spatial ETL tool (extract, transform, load) that is up to date and won´t break the bank, or my good mood. Most of what we need it for is pretty simple stuff like renaming, moving around and processing data in the tables and so on. We are designers and not software engineers, so writing our own scripts is to time consuming.

(Edited for clarity per request)


texts: How to import OSM data into postgis and add Elevation Information from raster

texts:
In my usecase I want to use OSM data for my routing server, but first, I'd like to add elevation data. And I'd prefer to do that in postgis as I'd need to do further analysis of the data using postgis.

So my question is which tools I should use for that task?

I guess I would use gdal to somehow import the elevation information from a raster-source? But how would a proper postgis sql command for that look?

EDIT:
My raster file is a .vrt file that combines about 3TB of data. So loading the whole raster into a DB will now work. I'd need a solution to go through each geometry in my DB and convert it from 2D to 3D.


texts: Adding multiple features from GDB files using PyQGIS

texts:
I tried to add all features from a GDB with PyQGIS.
    layer = QgsVectorLayer(r"path.gdb")
    QgsProject.instance().addMapLayer(layer)

This source code, only add one feature from GDB, but the GDB contains a lot of features.
Any idea, how can I do that?


texts: How to replace accented characters in the CQL filter, using GeoServer?

texts:
I have  GIS API based on OpenLayers and GeoServer. I want to perform an Ajax request in order to get all the features that have a given keyword.

For example, I could use the CQL filter:

name LIKE '%Setúbal%'

And I will get the feature corresponding to the portuguese district of "Setúbal".

Now, I wanted to simplify a bit and allow myself to write "setúbal".

I use the following CQL filter:

strToLowerCase(name) LIKE '%setúbal%'

It works fine, because the CQL function "strToLowerCase" successfully converts the text of the feature to lower case.

--- UNTIL HERE, EVERYTHING WORKS FINE... THE PROBLEM COMES NEXT ---

As you may have noticed, I want to simplify further by only writing "setubal" without the accented "u" ("ú").

I looked into the Filter Function Reference of Geoserver and I found the "strReplace" method, which seemed to be exactly what I needed to replace the accented letters of the feature's attribute.

So, I tried these CQL filters:

// Replace the accented "u" and convert to lower case
strToLowerCase(strReplace(name, ú, u, true)) LIKE '%setubal%'

// Only replace the accented "u"
strReplace(name, ú, u, true) LIKE '%setubal%'

But neither of them work and I get this parser error:

<?xml version="1.0" encoding="UTF-8"?>
<ServiceExceptionReport xmlns="http://www.opengis.net/ogc" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" version="1.2.0" xsi:schemaLocation="http://www.opengis.net/ogc http://schemas.opengis.net/wfs/1.0.0/OGC-exception.xsd">
    <ServiceException>Could not parse CQL filter list. Lexical error at line 1, column 18. Encountered: "\u00fa" (250), after : "" Parsing : strReplace(name, ú, u, true) LIKE '%setubal%'.</ServiceException>
</ServiceExceptionReport>

I thought that the problem might be in using an accented letter in a CQL filter, and tried this:

strReplace(name, e, u, true) LIKE '%setubal%'

And even so, I get a parser error (but different than the previous one):

<?xml version="1.0" encoding="UTF-8"?>
<ServiceExceptionReport xmlns="http://www.opengis.net/ogc" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" version="1.2.0" xsi:schemaLocation="http://www.opengis.net/ogc http://schemas.opengis.net/wfs/1.0.0/OGC-exception.xsd">
    <ServiceException>java.lang.RuntimeException: java.lang.NullPointerException: PropertyDescriptor is null - did you request a property that does not exist? java.lang.NullPointerException: PropertyDescriptor is null - did you request a property that does not exist? PropertyDescriptor is null - did you request a property that does not exist?</ServiceException>
</ServiceExceptionReport>

Am I using the "strReplace" function in a wrong way? Could anybody help with this issue, please? I want to replace accented characters in the CQL filter.


texts: Find densest point in satellite agriculture map

texts:
I have a .tiff image file with 3 bands.

I want to write a program, (python opencv?), that will find the first N densest points in the region. Densest, meaning the point with the most of its own colored points near it. For example, finding the densest points of red would look like:

I understand I should preprocess the image with filters. Is there any other methods I should apply to improve performance?

What algorithm would be most effective to find the densest points in this image?

Is it possible to get the set of neighbors for each dense point?

Is it possible to obtain the polygon shape from the densest point and its set of neighbor points?


texts:
In a typical GIS setup what are the different components that need to be there?

For example:

there needs to be a client program that can show the tiles and layers
there needs to be some GIS server that can query that layer database and return results that can be overlayed on MAP on client
Do we need a tile server for serving tiles or is it also job of GIS server?
There needs to be a database that saves all the GIS data
We can have shape files locally on the client machine

What other components can be there?


texts: Tool to search for repeated codes in ArcGIS 10.1

texts:
I am using ArcGIS 10.1.

I have a shapefile that has a field called "Code".

I want to make a Modelbuilder model to check if there are "duplicate" codes.

Is there a tool that could do this?


texts: WMTS : Convert geolocation (lat, long) to tile index, at a given zoom level?

texts:
I wanted to know how to get the indexes (x,y) of a WMTS tile for a given geolocation (latitude, longitude) and zoom level.

For exemple, I have a POI located at (48.675, 2.7), I want to get the corresponding open-street-map tile for the zoom 10.

Can I do the math ? Do I need a webservice ?
Precision : I have to do this programmatically.


texts: Union MVT features based on attribute / id in OpenLayers render

texts:
Following Paul Ramsey's article about ST_Subdivide I was able to split the Natural Earth boundaries dataset. Generating tiles which intersect the boundaries of Canada and Russia on the fly is now significantly improved.
How would I go about styling these tiles in OpenLayers?

We tried to ST_Union the MVT geometries but this process requires validation due to self intersections. The union geometries then must be dumped since geometry collections are not valid for MVT geometries. Even with all the steps in place there are still kinks due to PostGIS not always adhering to the winding order of geometries.
Ideally I would like to remove boundaries between features with the same ID inside the MVT as generated by the ST_MVT method.
Any pointer whether this is possible?
I am aware that mapbox-gl provides recipes to union features from MVT tiles based on an attribute but look for a way to achieve similar with ol/layer/VectorTile


texts: Display tabs don't appear when files and folders are selected in QGIS?

texts:
My setup: QGIS 2.8.1 Wien on OS X 10.10.2

I'm a complete beginner at GIS and I'm having trouble on getting the Param/Metadata/Preview/Attributes tabs to appear in the Display Window. When I select a folder or .shp file in the browser, nothing happens. The documentation I'm using shows this as what is supposed to appear when selecting a file or folder:

I'm getting this:

I've been able to find Metadata, Attributes, and Preview through other methods (Metadata by right clicking on the file, Preview by dragging file to Display Window, and Attributes by dragging file to Display Window and clicking the Open Attribute Table button), but I haven't been able to view Param at all.


texts: Showing route between two points in OpenLayers 3

texts:
How will I show the route between two points?

I have created points with below code.

var iconFeature = new ol.Feature({
          geometry: new ol.geom.Point([8637791.36, 1456487.82]),
          name: 'Police Station',
          population: 4000,
          rainfall: 500
        });
var vectorSource = new ol.source.Vector({
  features: [iconFeature]
});

var vectorLayer = new ol.layer.Vector({
  source: vectorSource,
  title: 'point_layer'
});


texts: Where to download QA bands (cfmask bands) for surface reflectance Landsat7

texts:
I downloaded scenes from both Landsat 7 & 8 reflectance. However, QA bands are included in Landsat 8 but not 7. I selected "all data" when downloading. Anyone has a clue where the QA bands are for Landsat 7 (further more, Landsat 4-5?)

QA for reflectance, USGS
This link says QA bands are available for all reflectance bands. Just couldn't find them.


texts:
I have list of 50 lists of pairs of SpatialLinesDataFrames. That means I have something like:

[[1]]

     [,1]    [,2]    [,3]
[1,] spldf1  spldf3  spldf5
[2,] spldf2  spldf4  spldf6

[[2]]

     [,1]    [,2]    [,3]
[1,] spldf1  spldf3  spldf5
[2,] spldf2  spldf4  spldf6

.
.
.
[[50]]

     [,1]    [,2]    [,3]
[1,] spldf1  spldf3  spldf5
[2,] spldf2  spldf4  spldf6

I want a list of 50 elements of pairs with the length of every pair. I did:

for(i in 1:50){
  for(j in 1:3){
      lengthpair[[i]]=sum(length(pairs[[i]][[1,j]]),length(pairs[[i]][[2,j]]))
  }
}

lengthpair is a list to store the sum of length of pairs every time.
pairs is the list with the pairs of SpatialLineDataFrames.

From this all I get is the sum of the length of last pair, of last run of my loop. What can I change to my code to do it?


texts: Given a handheld GPS with an accuracy of roughly 2 meters (verified through field tests) and an observation with an HDOP of 1.2, the accuracy of that one observation should be within 2 meters X 1.2 of actual horizontal position, correct? Where can I find documentation of how this works? That is documentation that shows reported device accuracy multiplied by the HDOP yields horizontal accuracy.


texts: Natural Earth doesn't render well in Robinson World projection

texts:
QGIS Pisa 2.10.1 Mac Pro, Mac OS

I don't seem to be able to create a simple World Map using the Natural Earth Dataset and the Robinson World projection, unless I zoom in. What could be the problem?


texts: How to queue FME output on the same workspace?

texts:
I have an .csv  coordinate file with the following header

num       xc    yc      xmin     ymin    xmax     ymax

I want to tranform  the x/y coordinates on each line of this file to lat / long (LL84)

I am using the following workspace

The output file is a csv file with the lat/long.
The problem I have is FME writing the lat/long to the output file in a different line for each pair of coordinates xc/yc , xmin/ymin, xmax/ymax and not in the same line as it is in the input csv. So the results look something like this,

Is it possible to get it to work  or I have to output three different csv files and merge them ?


texts: Adding buffers around each city of user defined region using geopandas?

texts:
How do I obtain a buffer around each city on the overlay map of cities on the United States area?

I have used the following piece of code to work this out but this does not work out it seems.

I am using geopandas in python to solve this:

 import pandas as pd
 import geopandas
 df = pd.DataFrame(
        {'City': ['Fargo', 'Orange', 'Jersey City'],
         'State': ['Nevada', 'California', 'New Jersey'],
         'Latitude': [46.877186 , 33.787914 , 40.728157 ],
         'Longitude': [-96.789803, -117.853104, -74.077644]})
    df['Coordinates']  = list(zip(df.Longitude, df.Latitude))
    df['Coordinates'] = df['Coordinates'].apply(Point)
    gdf = geopandas.GeoDataFrame(df, geometry='Coordinates')
    print(gdf.head())

world = geopandas.read_file(geopandas.datasets.get_path('naturalearth_lowres'))
#nybb=  geopandas.read_file(geopandas.datasets.get_path('nybb'))
# We restrict to United States Of America.
us_name = world[world.name=='United States']
ax = world[world.name == 'United States'].plot(
    color='white', edgecolor='black',figsize=(20,10))
us_states.plot(ax=ax,figsize=(20,10))

# We can now plot our GeoDataFrame.
gdf.plot(ax=ax, color='red')
ax.set(xlim=(-140,-50),ylim=(0,60))

gdf['Coordinates']=gdf.buffer(100)
gdf.plot()
gplot=geopandas.overlay(us_name,gdf,how='identity')
gplot.plot(edgecolor = 'k',alpha=0.5,cmap='tab10')

But what I want is buffers around each city on the map of the United States but what I get is just the map of US!

I am very new in this.


texts: Combining two shapefiles using common attribute

texts:
I have two shapefiles "sf1" and "sf2" for same area.

"sf1" is population data for year 2015 and has 6 variables.

"sf2" is population data for year 2010 and has 5 variables.

Both have one common attribute "GEOID". I want to get a combined layer with 8 variables (6 from 2015 and the 2 (variable and value) from 2010).
How can I combined the two shapefiles using the common attribute "GEOID"? I tried merging them.
sf3 <- merge(sf1, sf2)

It merges the two layers to give me 11 variables. This is not what I want. I want to join the attribute tables for the two layers and get one of the common attributes.


texts: convert SDE table to polygon feature class with Arcobjects

texts:
I figure there's a few ways to convert a table to polygon feature class, but what's the fatest and cleanest way to do it?
My table has lat, long, azimuth, and beam width and based on all this info I can create polygons.

Couple of methods I am thinking of right now:

Could I just add an OID field and
'shape' field to my existing table
and create polygons
Create an empty polygon fc, generate polygons and then copy all
table rows to fc.

The problem is, I don't really want attribute redundancy by having data in both the table and polygon fc.
Thanks!


texts: Inserting variable into layer name in PyQGIS?

texts:
I want to insert current date into name of new layer.
My current code to save layer with hardcoded name is like this:

_writer = QgsVectorFileWriter.writeAsVectorFormat(layer,r"D:\Maps\layer.shp","utf-8",None,"ESRI Shapefile")

I can get current date using:

cur_date_apndx = time.strftime("%d_%m_%Y")

statement, but I can not insert cur_date_apndx into the code above. I think "r" in the 2nd parameter of writeAsVectorFormat is the reason. If I remove "r" layer is not saved.

Could anyone help?


texts: GeoTIFF created using ArcGIS for Desktop not displaying in Google Earth or Photoshop?

texts:
I originally had an IMAGINE Land Cover raster file that I opened in ArcMap.

I used the Extract Feature from Raster to extract only the Developed class and saved it (still IMAGINE).

Then I used the Copy Raster to convert the IMAGINE to GeoTIFF, but when I try to view the GeoTiff in Google Earth or Photoshop, it doesn't display anything.

However, in ENVI, it displayed successfully.

What am I doing wrong?


texts: Showing vertices outside editing session in ArcMap?

texts: Does anybody know, if it is possible to see features' vertices (or all) when you are not in editing mode?


texts: Averaging three raster images (of same month obtained from three data centres) for approx 150 months using QGIS

texts:
How do I average the three raster images obtained from three data centre for each month of 15 years approx?
I tried to do it using Graphical Model in QGIS but it shows an error.
Error Log in QGIS
Processing algorithm 1/1…
Algorithm GRACE Combined Average starting…
Input parameters:
{'CSRLayer': 'G:/My '
'Drive/My_Data/ASHISH-Unsorted/PhD/Literature/Proposal/River '
'Basin/Data/GRACE/CSR/GRD-3_2003001-2003031_GRAC_UTCSR_BA01_0600_LND_v03.tif',
'GFZLayer': 'G:/My '
'Drive/My_Data/ASHISH-Unsorted/PhD/Literature/Proposal/River '
'Basin/Data/GRACE/GFZ/GRD-3_2003001-2003031_GRAC_GFZOP_BA01_0600_LND_v03.tif',
'JPLLayer': 'G:/My '
'Drive/My_Data/ASHISH-Unsorted/PhD/Literature/Proposal/River '
'Basin/Data/GRACE/JPL/GRD-3_2003001-2003031_GRAC_JPLEM_BA01_0600_LND_v03.tif',
'qgis:rastercalculator_1:GRACE Combined Average': <QgsProcessingOutputLayerDefinition {'sink':G:/My Drive/My_Data/ASHISH-Unsorted/PhD/Literature/Proposal/River Basin/Data/GRACE/Combined_Average/CA_GRD-3_2003001-2003031_GRAC_JPLEM_BA01_0600_LND_v03.tif, 'createOptions': {}}>}

An error occurred while performing the calculation
Error encountered while running Raster calculator
Algorithm GRACE Combined Average failed…
Execution failed after 0.64 seconds
Batch execution completed in 0.74 seconds
1 executions failed. See log for further details.


texts: Calculate the area of a polygon on the surface of the Earth

texts:
I'm trying to write a function that will accept an array of coordinates and calculate the area of the polygon.
function getArea(lat_long: Array<LAT_LONG>){

  lat_long.push(lat_long[0])
  const n = lat_long.length

  const cart = lat_long.map((point) => convertToCartesian(point));

  let sum1 = 0;
  let sum2 = 0;

  cart.map((point, i, cart) => {
    if(cart[i + 1]) {
      const {X: lat1, Y: long1} = point;
      const {X: lat2, Y: long2} = cart[i + 1]
      sum1 += <number>lat1 * <number>long2
    sum2 += <number>long1 * <number>lat2
   }
  })
  return Math.abs(sum1 -sum2) / 2;
}

I'm using the shoelace algorithm to successfully calculate the area of the polygon. I convert the lat/long coordinates into their Cartesian formats before running the algorithm.
Input of an array of coordinates:
[
     -30.755695508016164,139.58931179712357,
    -31.026745675780326,139.67720242212357,
    -30.929020409557967,139.9243948049360,
    -30.563139898321857,140.02052517602982,
    -30.500447229869152,139.79530544946732,
    ])

Convert the above to Cartesian:
const convertToCartesian = ({ lat, long }: LAT_LONG) => {
  const { a, feSq } = WGS84;

  const sinLat = Math.sin(convertRad(<number>lat));
  const sinLong = Math.sin(convertRad(<number>long));
  const cosLat = Math.cos(convertRad(<number>lat));
  const cosLong = Math.cos(convertRad(<number>long));

  const h = 0;
  const v = a / Math.sqrt(1 - feSq * sinLat * sinLat);

  const X = (v + h) * cosLat * cosLong;
  const Y = (v + h) * cosLat * sinLong;
  const Z = (v * (1 - feSq) + h) * sinLat;

  return { X, Y, Z };
};

Output:
[
      {
        X: -4177052.98836046,
        Y: 3556294.7743405653,
        Z: -3242647.5585449743
      },
      {
        X: -4170740.6233994444,
        Y: 3539899.576439874,
        Z: -3268435.3161851424
      },
      {
        X: -4190241.3515973534,
        Y: 3525463.391770664,
        Z: -3259145.5543680554
      },
      {
        X: -4212041.526141104,
        Y: 3531751.3820798383,
        Z: -3224283.8117593047
      },
      {
        X: -4200822.826425785,
        Y: 3550559.9915549974,
        Z: -3218297.530026407
      },
      {
        X: -4177052.98836046,
        Y: 3556294.7743405653,
        Z: -3242647.5585449743
      }
    ]

Results in 774321690.7109375 from my function.
I've compared this to the result from https://www.mathopenref.com/coordpolygonareacalc.html which outputs: 774321690.7119141.
However, when I compare this to a tool such as Google Maps it is incorrect (1513405095.71 m² is the answer). I should be applying the spherical properties of the earth to get an accurate answer but I'm not sure how. I've tried S = [θ - (n - 2)π]R² where θ is the sum of the angles, n is the number of corners.


texts:
I am using ArcMap 10.2 and have connected to my Oracle database, and can see the tables in Catalog. However, when I try to add a layer I get "Could not add the specified data object to the map. Underlying DBMS error[ORA-00936: missing expression]"

I am able to add this Oracle layer to other GUIs such as FME Data Inspector and QGIS. I have added both spatial and non-spatial indexes and have a primary key.

I have been able to connect to and show layers in other Oracle databases in our network, just not from this database. I connected to this database as a user that has select, update, delete, and insert permissions.

Any suggestions?


texts: ArcGIS API for JavaScript - error in Dojo parser.parse()

texts:
I started creating a Dijit based on a template HTML file. I started by referencing Dojo CSS and JS directly, and all was fine.

Later when I wanted to add a map into my dijit, first - changing nothing else - I switched my Dojo CSS links from googleapis.com to js.arcgis.com - and all was still fine.

But then I changed the Dojo JS script reference from googleapis.com to the ArcGIS API for Javascript. Doing this causes an error.

All else staying the same, this worked fine:

    <script src="http://ajax.googleapis.com/ajax/libs/dojo/1.9.1/dojo/dojo.js"></script>

    <script>
        require([
            "dojo/parser",
            "dijit/layout/BorderContainer",
            "dijit/layout/TabContainer",
            "dijit/layout/ContentPane",
            "dojox/layout/ExpandoPane",
            "dojo/domReady!"],
            function (parser) {
                parser.parse();
            });
    </script>

But simply changing out that first script tag with:

    <script src="http://js.arcgis.com/3.8/"></script>

I get: dojo/parser::parse() error TypeError: Unable to get property 'toString' of undefined or null reference when I run the page.

I can provide more detail, but I don't want to put you to sleep. Thanks!

Update: upon further toying, I see this doesn't happen every time the app starts. It happens ~3 starts out of 10 on one of my dev workstations, but every start on my other one. The first is IE10 with fast internet, the second is IE11 with moderate, whatever that's worth.

Chrome gives a little more info:


texts: Combining points that are close in space and time using R

texts:
I have table of data that contain 5 columns:

Unique_ID, Date_time, Lat, Long, Incident_Count

The sensors that supply this data often take what should be one occurrence, and split it up into several. eg I have three rows that look like this:

Unique_ID, Date_time, Lat, Long, Incident Count
xxx, 12/4/19 10:01, xxx, xxx, 2
xxx, 12/4/19 10:02, xxx, xxx, 4
xxx, 12/4/19 10:02,xxx, xxx, 3

and I would like to be able to turn that into:

xxx, 12/4/19 10:01, xxx, xxx, 9

Here's a piece of R script I've used in the past that will take two tables with lat long and time and tell me what from the second table is close in space and time to the first. What is different in this instance is that I want to go off of one table, and remove rows as they match with other rows before them in the data.

time_window = 72 # hours
dist_window = .3048 # km

incidents_list = list()

for (i in 1:nrow(Table1)) {
  if (i %% 100 == 0) {
    print(i)
  }
  Table1_incident = Table1[i,]
  incidents_close_in_time = subset(Table2, DATE_TIME >= Table1_incident$DATE_TIME - hours(time_window)
                             & DATE_TIME <= Table1_incident$DATE_TIME + hours(time_window))
  dist_matrix = spDistsN1(matrix(c(incidents_close_in_time$long, incidents_close_in_time$lat), ncol = 2),
                          matrix(c(Table1_incident$long, Table1_incident$lat), ncol = 2),
                          longlat = T)
 incidents_list[[i]] = (incidents_close_in_time$UNIQUE_ID[dist_matrix <= dist_window])
}
names(incidents_list) <- c(Table1$UNIQUE_ID)
Table1_incidents_df <- map_df(incidents_list, ~as.data.frame(.x), .id = "id")


texts: Choosing correct proj4string in R for a set of points

texts:
I was given a set of points which I eventually will overlay with a polygon shapefile in R. The issue is that the person who provided me with the points did not give them to in long/lat coordinates but the points were already projected (I will not be able to ask them to give me the points in a better format for a while). I was given the following information about the projection of the coordinates:

Top: 583309.373443 m
Bottom: -1996690.626557 m
Left: 221989.400576 m
Right: 3572989.400576 m

Projected Coordinate System: World_Sinusoidal
Projection: Sinusoidal
False_Easting: 0.00000000
False_Northing: 0.00000000
Central_Meridian: -76.00000000
Linear Unit: Meter

Here is the head of the csv I was given:

   [,1]     [,2]
1 2413489 -1995191
2 2416489 -1995191
3 2419489 -1995191
4 2422489 -1995191
5 2425489 -1995191
6 2428489 -1995191

I need to assign the correct projection now to read these points in correctly. Here is the pro4string and SpatialPoints code I have so far:

sinu_prj <- "+proj=sinu +x_0 +y_0 +units=m +no_defs"
brazil_fish_pts <- SpatialPoints(brazil_fish_grid2,  proj4string=CRS(as.character(sinu_prj)), bbox = NULL)

I know that the proj4string code is not correct but I don't know how to fix it? Can anyone be of help?


texts:
I tried to install GeoNetwork on a local Tomcat server. I'm using Tomcat 9.0.14 with GeoNetwork version 3.4.4.
I use Ubuntu version 18.04 and Firefox as a web browser.

Tomcat is up and running.
I installed the GeoNetwork application with a .war file copied in the webapps folder. When I deploy it I don't have any error message.

In the localhost:8080/manager page I can see GeoNetwork installed with the status "running=true".

However,when I try to access it via localhost:8080/geonetwork I receive an error 404.

In the catalina.out log file I don't see any error message.

Any idea of the problem that could cause this? I read about every post I could find but I don't find the answer.

Edit:
One of the problem was missing permission for the folder geonetwork (execute I think). Now I have an empty yellow page instead of error 404 and the geonetwork log speak to me, which can be helpful.

  2018-12-19 09:11:10,853 FATAL [jeeves.engine] - Raised exception during init
  2018-12-19 09:11:10,853 FATAL [jeeves.engine] -    Exception : org.springframework.context.ApplicationContextException: I/O error parsing bean definition source for Root WebApplicationContext; nested exception is java.nio.channels.ClosedByInterruptException
  2018-12-19 09:11:10,853 FATAL [jeeves.engine] -    Message   : I/O error parsing bean definition source for Root WebApplicationContext; nested exception is java.nio.channels.ClosedByInterruptException
  2018-12-19 09:11:10,856 FATAL [jeeves.engine] -    Stack     : org.springframework.context.ApplicationContextException: I/O error parsing bean definition source for Root WebApplicationContext; nested exception is java.nio.channels.ClosedByInterruptException

These are the first lines. Then it fails at creating the db

  2018-12-19 09:12:55,904 ERROR [org.hibernate.tool.hbm2ddl.SchemaUpdate] - HHH000388: Unsuccessful: create table Address (id integer not null, address varchar(255), city varchar(255), country varchar(255), state varchar(255), zip varchar(16), primary key (id))

And finally it repeats that the database is read only. I can't find the file geonetwork.h2.db which should at least exist?!
Should I create it myself first to that geonetwork can write in it?

Edit 2 : Somehow it managed to create the geonetwork.h2.db in the /webapps/geonetwork/WEB-INF/config-db folder. After giving full permission it finally works.


texts: Conversion of .shp to EPSG:4326 has slight offset

texts:
I just started this project where I need to project the Dutch municipalities on a map. Luckily this seems very easy with the help of GeoJSON. However, I'm having a slight issue where my conversion to EPSG:4326 seems slightly off and I don't know why or how to correct it (The original .shp-file is in the Dutch Rijksdriehoekscoördinaten system).

My steps:

I downloaded the most recent shapes from cbs.nl: http://download.cbs.nl/regionale-kaarten/shape-2015-versie-1-0.zip

I used OSGeo4W to transform the Gem_2015.shp file to GeoJSON with the option -t_srs EPSG:4326 to convert to lat/lon:

ogr2ogr -f GeoJSON gem.geojson Gem_2015.shp -t_srs EPSG:4326

This results in a nice geojson file. But if we take a closer look at one entry, it seems that it has a slight offset:
https://gist.github.com/Waterstraal/95f73c22662775e7169d079f222869cb/60ae6e504efcd929bfebf8ae23fdb173038a7604

Am I using a wrong projection, or am I doing something else wrong?


texts:
I have recently installed QGIS and trying to work with it but I don't see this option "Get More" at all. As suggested by some previous source, I also checked whether I am behind a proxy but that isn't the case.

What could be the cause?


texts: ArcGIS Pro tool is failing. "Set LAS Class Codes Using Features"

texts:
I'm using ArcGIS Pro and I have a nice LAS Dataset. It came to me with the points already classified. The classifications are off, many buildings are classified as vegetation, for example.

I'm trying to use a building footprints layer to re-classify the points. For some reason the tool keeps failing. It doesn't seem to matter if my LAS Dataset is large or 1 square. It fails every time. I do get small results from it once in a while.

Any ideas why the tool is failing?


texts: What is QgsCoordinateTransformContext and how to set it up for distance calculation in QGIS?

texts:
When setting source coordinate system for QgsDistanceArea with setSourceCrs() you now (QGIS 3) need to specify context parameter. Can someone explain what is this context parameter, and how to set it up for basic distance calculation (point to point, and they are in different crs)? The QgsDistanceArea doc only has the basic info so far.
This parameter is non-optional, so not providing it results in an error, and none/'' parameters are not accepted..

    TypeError: QgsDistanceArea.setSourceCrs(): not enough arguments


texts: No Data was found error after applying a MaxEnt classifier (GEE)

texts:
I am trying to redo this tutorial: https://developers.google.com/earth-engine/classification, using the MaxEnt classifier with my data.

Code: https://code.earthengine.google.com/0ce8a254cd46d7b8d0384d9fdd2d1675

//*******************************************************************************************
// Set Study Area
var countries = ['HU', 'BU', 'IT', 'LO', 'AU', 'EZ'];
var geometry = ee.FeatureCollection('USDOS/LSIB/2013').filter(ee.Filter.inList('cc', countries));
Map.addLayer(geometry.draw({color: 'ffffff', strokeWidth: 3}), {},'Study Area');
Map.centerObject(geometry, 5);

//*******************************************************************************************
// Import Data

var points = ee.FeatureCollection("users/farhatlokmen/ecoModeling/falcoCherrug");
var current245 = ee.Image("users/farhatlokmen/ecoModeling/current245").clip(geometry);

//*******************************************************************************************
// MaxEnt

var bands = ['b1', 'b5', 'b6', 'b12', 'b13', 'b14'];
var label = 'SPECIESCODE';

// Sampling
var sampleData = current245.select(bands).sampleRegions({
  collection: points,
  properties: [label],
  scale: 500
});

// Build a classifier
var classifier = ee.Classifier.gmoMaxEnt(0,  0.00001, 0.00001, 0, 50)
                            .setOutputMode('PROBABILITY')
                            .train(sampleData, label, bands); // training

// Classify the image with the same bands used for training.
var classified = current245.select(bands).classify(classifier);

// Display the inputs and the results.
Map.addLayer(classified, {min: 0, max: 1, palette: ['black', 'red']}, 'classification');

Error Message: classification: Layer error: No data was found in classifier training input.


texts: Can I download the current version of QGIS if I am running Mac 10.8.2?

texts: I am taking an introductory course in GIS at my university and my teacher recommended QGIS for Mac users. Unfortunately, I haven't been able to download the software. I just upgraded to Mac OS 10.8.2 recently, and I am not sure if Lisboa can work on Mountain Lion. Please let me know what I need to do in order to run QGIS on my laptop.


texts:
I want to change data type of a field to type esriFieldTypeOID in a layer. So i copy existing field into a new field than delete existing field first and than change datatype for copied field and than add copied filed with new data type as new field. Please do note that currently that field holds values in it. But below code throws an error saying "A syntax error occurred trying to evaluate a query string (Exception from HRESULT: 0x80040203)" as soon as DeleteField function is called.

 IDocument doc = ArcMap.Application.Document;
        IMxDocument mxDoc = doc as IMxDocument;
        IMap map = mxDoc.FocusMap;

        ILayer layer;
        for (int i = 0; i < map.LayerCount;i++ )
        {
            layer = map.Layer[i];
            IFeatureClass featureClass = (layer as IFeatureLayer2).FeatureClass;

            IFields2 fields = featureClass.Fields as IFields2;
            IField keyField = null;
            IField tempKeyField = null;
            int fieldIndex = fields.FindField("layerField");
            if (fieldIndex >= 0)
            {
                keyField = fields.Field[fieldIndex] as IField2;
                tempKeyField = keyField;
                featureClass.DeleteField(keyField);
                IFieldEdit2 editedKeyField = (IFieldEdit2)tempKeyField;
                editedKeyField.Type_2 = esriFieldType.esriFieldTypeOID;
                featureClass.AddField(editedKeyField);
            }

        }


texts:
I have been working with two layers and I am trying to join them and turn them into a raster data. The first layer is a watershed layer which I have assigned each shed a location from the center point (green dot) for instance shed 1,2,3, and so on. The second layer is distance rings which have a certain distance from the center point. I joined them using the Union tool. After, I use the Calculate Field tool to assign probability to each polygon. For intance if a polygon is located on the first water shed and is on the 1.3 mile distance ring it gets a probability of 7.042 and so on.

Now what I am having trouble with is working with the Summary statistics tool and the Join tool. By using the summary statistics tool I am trying to find the total area that each colored section has. For instance to find the total area of the yellow polygons or the total area of the orange polygons and so on. I am doing this because I want to be able to divide a single polygon that is located on the same colored section by that colors area. For instance the area of the top left yellow colored polygon by the entire area of the yellow colored section. I figured the summary statistics tool would be able to give me that information.

Then by using the Join Field tool I would like to be able to join that area with the statistics that I assigned it before with the Calculate Field tool.

How do I get the two to join?

I have provided two images of my watershed and distance rings united and the modelbuilder model


texts:
How to set NULL value by arcpy.da.InsertCursor on date column?

I try insert "NULL", None, '', 0

dictIN= [{'field1':"test1",'date1':'01-02-2015','date2':''},\
 {'field1':"test2",'date1':'01-02-2015','date2':'01-02-2015'},\
 {'field1':"test3",'date1':'','date2':''},\
 {'field1':"test4",'date1':'01-02-2015','date2':''}]

def setNULLIf(field):
    if field == '':
        return None #try NULL,'',0
    else: return field

cursor = arcpy.da.InsertCursor(GDB,("field1", "date1", "date2"))
for x in dictIN:
    cursor.insertRow((x['field1'], setNULLIf(x['date1']), setNULLIf(x['date2'])))

Is it possible by arcpy.da ?
I use ArcGIS 10.2. I try insert data to GDB file.


texts: Reducing a multi-band image and keeping the bands

texts:
I would like to calculate the area of all non-zero pixels on multiple bands of an image. My idea was to use

var areaImage = image.multiply(ee.Image.pixelArea());

To do that I will need an image where pixel values equal 1, so I was trying to use

image.reduce(ee.Reducer.anyNonZero());

However, since I am using a multi-band image, I receive a flattened image with only 1 band.

How could I assign all non-zero pixels a value of 1 across all bands without losing the bands? Or is there a smarter way to go about it?

Here is the code I tried:

//----load MCD64A1 Version 6 Burned Area and create a multi-band image

var datasetAll = ee.ImageCollection(MODISglobal500)
                  .filter(ee.Filter.date('2000-01-01', '2019-12-31'))
                  .filterBounds(border);

var DSall1B = datasetAll.select('BurnDate')

var DSallImage = DSall1B.toBands()

//----convert pixels to value = 1

var ImageValue1= DSallImage.reduce(ee.Reducer.allNonZero())

//----create an image with pixels value = area

var ImageValueArea = ImageValue1.multiply(ee.Image.pixelArea())


texts: Reduce the height of the selected cesium entity to zero on mouse click in cesium

texts:
How do i need to reduce the height of the selected entity to zero on mouse click in cesium. Below is the code which i had tried
var handler = new Cesium.ScreenSpaceEventHandler(viewer.scene.canvas);
handler.setInputAction(function(click){
    var pickedObject = viewer.scene.pick(click.position);
    if (Cesium.defined(pickedObject)) {
        var entityId = pickedObject.id._id;
        pickedObject.id._polygon._extrudedHeight._value = 0
    }
},Cesium.ScreenSpaceEventType.LEFT_CLICK);


texts: Get an XML response with a Java HTTP POST request on GeoNetwork

texts:
I try to complete automatically my GeoNetwork catalog with a Java client. My metadata respect the ISO 19115 schema and the interactions are performed trough CSW.

To insert metadata, I send HTTP POST request. The content of my request respects the CSW XML request format.  No matter the library (io, apache commons, geotools), I always have the same HTML response as follows:

<!DOCTYPE html
  SYSTEM "html">

      <meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
      <title>My GeoNetwork catalogue - My organization</title>
      <meta charset="utf-8">
      <meta name="viewport" content="initial-scale=1.0, user-scalable=no">
      <meta name="apple-mobile-web-app-capable" content="yes">
      <meta name="description" content="">
      <meta name="keywords" content="">
      <link rel="icon" sizes="16x16 32x32 48x48" type="image/png" href="../../i$
      <link href="rss.search?sortBy=changeDate" rel="alternate" type="applicati$
      <link href="portal.opensearch" rel="search" type="application/opensearchd$
      <link href="/geonetwork/static/gn_login_default.css" rel="stylesheet" med$
      <link href="/geonetwork/static/bootstrap-table.min.css" rel="stylesheet" $
      <link href="/geonetwork/static/ng-skos.css" rel="stylesheet" media="scree$
      <link href="/geonetwork/static/srv_custom_style.css" rel="stylesheet" med$

      <div class="navbar navbar-default gn-top-bar" data-ng-hide="layout.hideTo$


texts: How to add GPS UI on CARTO public_map and/or embed_map?

texts:
I want to add GPS UI for getting current location on CARTO public_map and/or embed_map, when they are accessed from mobile phones.
I know I can do it coding in JavaScript and carto.js, but I want to do it without coding.

Is it possible?


texts:
I am working in a navigation company, and they have purchased server license for hundreds of computers, but ArcObjects is not installed in any of the system.

So I need a loose copy of ArcObjects SDK for learning purpose.

I am not able to find it on Esri website, there is no download option for Arcobjects SDK. Instead ArcGIS RunTime SDK is available to download.

Is it available on internet freely?


texts: Renaming the value of a proportional symbols layer in Table of Contents in ArcMap

texts:
In the case of a - for example - graduated colors or graduated symbols layer, it is easy to rename the value in the Table of Contents. I always click the value name and press F2.

However, pressing F2 on the value name of a proportional symbols layer in ArcMap doesn't give you the opportunity to rename the value. Is there any other way to do so?

So, I'm trying to change the name of the 'PURBAN'-value in the second layer.


texts: Getting non overlapping bbox in OpenLayers?

texts: I have tree bbox as follow: bbox 1, bbox 2, and bbox 3. Now I want to get Final bbox from bbox 1 in OpenLayers 5. How?


texts:
I'm trying to merge 2 geodataframes and combine all entities.
So I use :

cloudNshadow = gpd.overlay(cloud, shadow, how='union')
cloudNshadow = cloudNshadow.unary_union()

But it returns me :

cloudNshadow = cloudNshadow.unary_union()
TypeError: 'MultiPolygon' object is not callable

So I insert a gdf.explode() :

cloudNshadow = gpd.overlay(cloud, shadow, how='union')
cloudNshadow = cloudNshadow.explode()
cloudNshadow = cloudNshadow.unary_union()

But same error...
Explode() isn't supposed to explode multi-entities ??


texts:
According to the specifications, GSV sentence contains information about satellite. But in this sentence (extracted from my Sierra receiver):
$GPGSV,3,3,11,04,,,,08,04,262,,10,06,151,*75

If I understand right:

3-> Number of sentences
3-> Index (3rd sentence)
11 -> Satellites being tracked
04 -> PRN1
,, -> Elevation Null
,, -> Azimuth Null
,,    -> SNR Null
08 -> PRN2
04 -> Elevation ?
262 ->Azimuth
,, -> SNR Null
10 -> PRN3
06 -> Elevation?
151 -> Azimuth

Where is SNR for PRN3? If it is null...why are the commas missing??
*75 -> Checksum
What do you think?


texts:
I get an error message if I try to install some symbols in QGIS with QGIS Resource Sharing Plugin. My System is Windows 10 and QGIS 2.18.19

Could someone give a hint to solve that Problem?


texts: How can I code a 'zoom to layer' applied to a vector layer when I build a QGIS standalone Python application?

texts:
I am trying to code a standalone QGIS application that is able to print a composed map that zooms in the same layer every time I call it. The feature that have to be zoomed in change its location when is regenerated so my map application has to zoom in to a different place when I call it.
I have been preparing an Atlas composer but when I call the print Python function it keeps the same area, it does not zoom to the chosen Atlas Layer.
Another possibility is to zoom in to my layer every time I call the mapper, but I can't achieve it.
I am adapting a code I found in the following Git repo without success: https://github.com/rjspiers/qgis-standalone-map-export

from qgis.core import *
import sys

print 'Hello from QGIS'

from qgis.core import (QgsProject, QgsComposition, QgsApplication, QgsProviderRegistry)
from qgis.gui import QgsMapCanvas, QgsLayerTreeMapCanvasBridge
from PyQt4.QtCore import QFileInfo, QObject
from PyQt4.QtXml import QDomDocument

gui_flag = True
app = QgsApplication(sys.argv, gui_flag)

# Make sure QGIS_PREFIX_PATH is set in your env if needed!
app.initQgis()

# project.qgs
project_path = r'D:\Fmesalessimulation\Simulator\qgis-standalone-map-export-master\smartfee_map.qgs'  # absolute filepath example: 'C:/folder/subfolder/project.qgs'

# template.qpt
template_path = r'D:\Fmesalessimulation\Simulator\qgis-standalone-map-export-master\smartfee_maptemplate.qpt'

#LayerIds
unused = 'BOUGHT_PAYS_DIFF20170130111140143'
sailedthrough = 'BOUGHT_UNITS20170130111140159'
optimalunits = 'PAYS_UNITS20170130111140165'
boundingbox = 'BOUNDING_BOX20170130111140161'

def make_pdf():
    canvas = QgsMapCanvas()
    # Load our project
    QgsProject.instance().read(QFileInfo(project_path))

    # Set canvas extent
    # canvas.setExtent(QgsRectangle(xmin, ymin, xmax, ymax))

    # Load layers here if they are not already in the project
    for layer in QgsMapLayerRegistry.instance().mapLayers().values():
        if layer.name() == 'BOUNDING_BOX':
            lyr = layer
            break
    # lyr.setSubsetString(filter_expression)
    #################
    list_layers = []
    for f in lyr.getFeatures():
        if f['unit_name'] != '':
            list_layers.append(f.id())
    print list_layers
    lyr.setSelectedFeatures(list_layers)
    box = lyr.boundingBoxOfSelected()
    lyr.setSelectedFeatures([])
    canvas.setExtent(box)
    canvas.refresh()
    canvas.show()
    #################
    # bridge used in standalone script: http://docs.qgis.org/testing/en/docs/pyqgis_developer_cookbook/loadproject.html
    bridge = QgsLayerTreeMapCanvasBridge(
        QgsProject.instance().layerTreeRoot(), canvas)
    bridge.setCanvasLayers()

    # Check data providers load ok
    print('Provider List')
    print(QgsProviderRegistry.instance().providerList())

    # Print layer validity - I've had some trouble with some wms layers
    for layer in QgsMapLayerRegistry.instance().mapLayers().values():
        print(layer.name())
        print(layer.isValid())

    template_file = file(template_path)
    template_content = template_file.read()
    template_file.close()
    document = QDomDocument()
    document.setContent(template_content)
    composition = QgsComposition(canvas.mapSettings())
    # You can use this to replace any string like this [key]
    # in the template with a new value. e.g. to replace
    # [date] pass a map like this {'date': '1 Jan 2012'}
    substitution_map = {'title': 'PyQGIS test'}
    # substitution_map = ''
    composition.loadFromTemplate(document, substitution_map)
    # You must set the id in the template
    map_item = composition.getComposerItemById('Main Map')
    map_item.setMapCanvas(canvas)
    map_item.zoomToExtent(canvas.extent())
    # You must set the id in the template
    legend_item = composition.getComposerItemById('Legend')
    legend_item.updateLegend()
    composition.refreshItems()
    composition.exportAsPDF(r'D:\Fmesalessimulation\Simulator\OUTPUT\JPG\PaysComparison.pdf')
    QgsProject.instance().clear()

make_pdf()


texts: Get multiple values from summary tables in Model builder

texts:
I working with ModelBuilder and  I have two summary statistics tools,  with 'case fields' that have multiple sum values, I want to use these sum values from one table and divide the sum values from the other table to get percentages,I have added a field in one table to store the percentages.
I am using Iterate Field Values tool to get the output values that I can input into Calculate field to get the percentage, it works but only for the first record, obviously because it is only getting one value from the summary tables, I tried Collect values tool to get multiple values, but I get an error 000539, I using the code from this discussion:
How to extract multiple values out of a table in ModelBuilder

I have more than than 2 values so I guess List and not index is required.

The whole model runs I get all the results, except the error with collect Value tool. Here is snap shot of the model,
I will separate the iterator tool and add it as submodel once it works.

Here are the input tables and I want to calculate PctCtyPop using calculate field by dividing  the 2 SUM_POPULATION,
 it giving me the right value for the 1st record 5505/20202 but for the second instead of getting 118774/131506, it is giving some thing else.


texts: Populating Attribute field using Python code without copying over null rows?

texts:
I am trying to copy the values from one column within my attribute table (!roomName!) over to another column (!spaceDiscriminator!). I understand how to do an easy field calc over for all the values but the problem is, I only want the rows in !roomName! that have values to copy over. If !roomName! contains a Null row then I do not want it copying over, instead I want the original value in !spaceDiscriminator! to remain and not get edited. I cannot figure out how to do this adequately. I have attached my code below:

In the first if statement I have tried many alternative scenarios. If x == None: seems to be the one that accurately selects the Null fields I don't want copying over. return None still edits the !spaceDiscriminator! field however. I've tried to use continue there but it throws and error. Is there any simple statement I can use there that basically says, "if the value is none/Null continue to the next row and leave original !spaceDiscriminator! value".


texts: Calculating 3D surface area base on slope in ArcGIS Desktop?

texts: How can I calculate 3D surface area base on slope in GIS for each morphological units in watershed?


texts: How can I retrieve private tag GDAL_NODATA from a TIFF using libtiff?

texts: In an earlier post I reported about my attempts to read TIFF with data stored in tiles using python and libtiff. First I tried to use libtiff's function TIFFReadEncodedTile but finally I tried function TIFFReadTile. I have been able to make progress now. The question I have today: how can I retrieve the value of private / custom tag GDAL_NODATA. Both GDAL and ArcGIS make use of this tag. I've seen examples of how to do this in C++ and C#. I still wonder how to do this in Python though. I guess it would be possible by means of the class TIFFExtender and the function add_tags - both defined in module libtiff_ctypes. Could anybody provide an example?


texts: ArcMap Field Calculator - find and add last number

texts:
I'm trying to get a part of a string (numbers after letter) in a new table (field Description) using field calculator, but I need it only from a field next to it (field MUID) containing the first letter "P" while the rest, that contains something else will be null.
That means if MUID = P14 then Desc = 14 or MUID = P1 then Desc = 1 and if MUID = HT9 then Desc = 0.
Is there a way to write it in the Python parser?
Will need it in the future to implement it via ModelBuilder.
I tried to use this, but it keeps getting Error signs


texts: confused about distances on the wgs-84 ellipsoid

texts:
As far as I understand given a certain latitude the length of a degree of longitude is constant.

So using any distance calculator (I use Mathematica here) this equals:
(circumference of the earth)

GeoDistance[{0, 0}, {0, 1}]*360 = 40075km
GeoDistance[{0, 0}, {0, 90}]*4 = 40075km

However, if we change the latitude to 15° for example, this happens:

GeoDistance[{15, 0}, {15, 1}]*360 = 38718.1km
GeoDistance[{15, 0}, {15, 90}]*4 = 38372.7km

Are these roundoff errors? Which one is correct?


texts: Nominatim installation error: Out of memory for node cache dense index

texts:
I am installing Nominatim on my Ubuntu server 14.04 through this link and trying to load OSM data into my postgreSQL database. For loading into my database I am using this command as given in the link:

~/Nominatim/utils/setup.php   --osm-file data/latest.osm.pbf   --all --osm2pgsql-cache 1024 2>&1   | tee setup.log

But it is throwing an error

Using projection SRS 4326 (Latlong)
NOTICE:  table "place" does not exist, skipping
NOTICE:  type "keyvalue" does not exist, skipping
NOTICE:  type "wordscore" does not exist, skipping
NOTICE:  type "stringlanguagetype" does not exist, skipping
NOTICE:  type "keyvaluetype" does not exist, skipping
NOTICE:  function get_connected_ways(pg_catalog.int4[]) does not exist,  skipping
Allocating memory for dense node cache
Out of memory for node cache dense index, try using "--cache-strategy  sparse" instead
Error occurred, cleaning up
ERROR: Error executing external command: /home/nominatim/Nominatim /osm2pgsql/osm2pgsql -lsc -O gazetteer --hstore -C 256 -P 5432 -d nominatim /home/nominatim/data/latest.osm.pbf
Error executing external command: /home/nominatim/Nominatim/osm2pgsql/osm2pgsql -lsc -O gazetteer --hstore -C 256 -P 5432 -d nominatim /home/nominatim/data/latest.osm.pbf

My ubuntu has 7GB of RAM and I am not sure why this error is coming. I ran the above command by reducing the cache size to 512 and then 256 but it is running into same error. I am not sure how to use --cache-strategy  sparse in the given command.

Screenshot of my Error:

Could anybody help with the setup?


texts: Osmosis failed: ERROR: relation "pend_ways" does not exist

texts:
I tried to generate osm update using osmosis and osmium(omsium only for setting configuration.txt from existing russia-latest.osm.pbf).

I use osmosis like that:

osm2pgsql -append --slim --cache 10000 --proj 3857  --append -H  local -d osm -U osm_owner -W -P 5432  --input-reader xml -k change.osc

There is reuslt:

Password:
Allocating memory for dense node cache
Allocating dense node cache in one big chunk
Allocating memory for sparse node cache
Sharing dense sparse
Node-cache: cache=10000MB, maxblocks=160000*65536, allocation method=11
Mid: pgsql, cache=10000
node cache: stored: 0(-nan%), storage efficiency: -nan% (dense blocks: 0, sparse nodes: 0), hit rate: -nan%
Osm2pgsql failed due to ERROR: PREPARE mark_ways_by_node(int8) AS select id from pend_ways WHERE nodes && ARRAY[$1];
PREPARE mark_ways_by_rel(int8) AS select id from pend_ways WHERE id IN (SELECT unnest(parts[way_off+1:rel_off]) FROM pend_rels WHERE id = $1);
 failed: ERROR:  relation "pend_ways" does not exist
LINE 1: PREPARE mark_ways_by_node(int8) AS select id from pend_ways .

my configuration.txt

baseUrl=http://download.geofabrik.de/russia-updates
maxInterval=60

my state.txt

#Thu Nov 21 16:44:37 MSK 2019
sequenceNumber=1473
timestamp=2019-11-20T21\:18\:03Z

What is this? It is a bag? I didn't find  pend_ways in osm2pgsql


texts: Determining which polygon(s) each point lies in using QGIS?

texts:
I am using QGIS software. I have shapefile layer of points and a shapefile layer of polygons(many of them overlapping). I need to determine which polygon(s) each point lies in.

I tried "Join attributes by location" with the solution posted here:(How to find the polygon inside which a point lies?). My question is different because many of my points will be inside multiple polygons. This only gives one polygon per point. My solution will have many polygons for each point.

For example, I would like to know the identity of every blue circle that lies underneath each black dot.

How can I complete this analysis?


texts:
I'm trying to clip and resample a raster to match the dimensions and resolution of another raster. Here's my gdalwarp command:

    gdalwarp  -t_srs EPSG:4326 -te 102.3375307079206 14.6874037407838 107.6323189079206 10.350077240783799 -tr 0.0008332999999999999 0.0008333000000000003 -r bilinear input.tif output.tif

The extents and resolution are grabbed from the other raster and are correct. However, I get the following error:

    ERROR 1: Attempt to create 6354x-5204 dataset is illegal,sizes must be larger than zero.`

The problem is the height of the raster -5204. I've tried supplying the x and y resolution -tr with negative and positive arguments with no avail. Have also tried putting (abs) around all variables to avoid negative values. Here is the gdalinfo call for my input.tif file

    Driver: GTiff/GeoTIFF
    Files: input.tif
    Size is 6187, 8475
    Coordinate System is:
    GEOGCS["WGS 84",
DATUM["WGS_1984",
    SPHEROID["WGS 84",6378137,298.257223563,
        AUTHORITY["EPSG","7030"]],
    AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433],
AUTHORITY["EPSG","4326"]]
    Origin = (97.277217688999997,23.425182730000000)
    Pixel Size = (0.002105966899941,-0.002105966899941)
    Metadata:
      AREA_OR_POINT=Area
    Image Structure Metadata:
      INTERLEAVE=BAND
    Corner Coordinates:
    Upper Left  (  97.2772177,  23.4251827) ( 97d16'37.98"E, 23d25'30.66"N)
    Lower Left  (  97.2772177,   5.5771133) ( 97d16'37.98"E,  5d34'37.61"N)
    Upper Right ( 110.3068349,  23.4251827) (110d18'24.61"E, 23d25'30.66"N)
    Lower Right ( 110.3068349,   5.5771133) (110d18'24.61"E,  5d34'37.61"N)
    Center      ( 103.7920263,  14.5011480) (103d47'31.29"E, 14d30' 4.13"N)
    Band 1 Block=6187x1 Type=Int16, ColorInterp=Gray

Any ideas how I can resolve this problem?


texts:
Using ArcMap, I am trying to add a series of co-ordinates to a map to show the distribution of shelf seas I am analysing. However, when I go through the process of adding XY data to the basemap, it will only plot 0, 0 (which isn't even one of my data points) but none of the actual 40 points.

The data is on an xls file, the lat and long is as decimals, with no spaces, and with - for south and west data.
I have tried projecting in a number of different coordinate systems (WGS 84, etc) but none produce any different result.

Also, despite there being data in the file, the attribute table is always empty (which is probably the issue, but I can't work out how to fill it).


texts:
We need to do are very simple administrative boundary drawings over a (insert fav map here) via some website.

Our requirements are:

click to create points which then form a polygon
drag existing points to edit a polygon

I know Google Maps and others can do all of this but we are hoping to leverage an existing site that can already offer this.

We would also love to 'load' some WKT data, so we can edit that (if required), where load means paste into a textbox.

Anyone know of any sites out there with this functionality?


texts: Comparing two spatial datasets statistically

texts:
I'm new to this forum. I'm working on the following problem -

I have two datasets. Let's call it dataset A and B. Both of them will contain a set of points like

A1,LonA1,LatA1
A2,LonA2,LatA2
..
..

and

B1,LonB1,LatB1
B2,LonB2,LatB2
..
..

I also have link data for both A and B like A1->A2, A1->A3. So, I can essentially create a graph. I want to statistically compute how similar the two datasets (graphs) are. I want to say something dataset A is X% similar to dataset B.

I'm using coupled node edge scoring to compute similarity matrix and then use Hungarian algorithm to compute similarity between the graphs created. I was wondering if there is something simpler that I could do to do this. Thanks!


texts: Data Driven Pages Navigation in a Script

texts:
Is there a method to simply progress from one data driven page to the next in a script?  Essentially, like pressing the "Next Page" button on the toolbar.  I have started with:

import arcpy
mxd = arcpy,mapping.mapdocument ("Current") # or should I reference the mxd directly?
ddp = mxd.datadrivenpages
cp  = ddp.currentPageID
np  = cp + 1              #Advance to next PageID?
arcpy.RefreshActiveView()

But cannot seem to find a way to advance to the next page, most searches yield methods to only export the pages.

The end goal is to build a tool that advances to the next data driven page and selects all features in a layer that has a page definition set to match the index.  This will facilitate attribute editing for the selected features.

The image below shows the selection model which works just fine, I just need to be able to script something to move forward a page and then run the Model (or script it all) to simplify the process.


texts: Using threading module with Arcpy: Manager/Worker architecture

texts:
I asked this question on Stack Overflow on how to setup a manager/worker scenario.

But I am getting some error when trying to use threading on my Python install with ArcGIS 10 (just the python install from ArcGIS 10.0).

Is anyone using threading successfully in arcpy.

My end goal is to be able to poll an arcgisjobs directory on ArcGIS Server 10.0 and call worker processes for each new folder.  Each new process should be called and set to go and return control back to the manager to call more processes before previous ones are finished.


texts: What is the difference between using a PolyLine layer and a GeoJSON layer in leaflet?

texts:
This question arises because on my project I will have a backend that consists of gpx tracks. Then, my frontend queries it (based on time, location filters, etc.) and I will either get the corresponding tracks and show them on a leaflet map or some specific points to also show on that map.

Then, I need to be able to style those tracks and be able to interact with them. Is it easier to achieve this through PolyLines or GeoJSON layers?


texts: Reset particular layer in Google Earth Engine

texts:
I have two layers in my GEE and I want to reset one of them base on my selector. I have the code as following:

Map.layers().reset(['UHI']);
Map.addLayer(UHIAdded, {opacity: 0.5}, 'layer1');
Map.addLayer(uhiday, {min: -2.0, max: 3.0, palette: palette, opacity: 0.5}, 'UHI');

And there's the error:

  Cannot add an object of type string as a layer.

If I get rid of the [], it will give me the error:

  The argument provided to ui.List.reset() must be an array.

How do I specify the layer that I want to reset?


texts:
I have a variogram and i'm trying to fit an exponential variogram Variogram

and i try to estimate the parameter of this variogram through the data:

    vk=variogram(z~1, data=data)
    nugget_ini=.25*median(c(rep(vk$gamma[1],vk$np[1]),rep(vk$gamma[2],vk$np[2])))
    model=c("Exp")     sill_ini=median(c(rep(vk$gamma[length(vk$gamma)-1],vk$np[length(vk$gamma)-1]),
                            rep(vk$gamma[length(vk$gamma)-2],vk$np[length(vk$gamma)-2]),
                            rep(vk$gamma[length(vk$gamma)-3],vk$np[length(vk$gamma)-3]),
                           rep(vk$gamma[length(vk$gamma)-4],vk$np[length(vk$gamma)-4])))-nugget_ini
        range_ini=vk$dist[min(which(vk$gamma>.95*(sill_ini+nugget_ini)))]/3
        vm=vgm(sill_ini,model[1],range_ini[1],nugget_ini[1])
        vk.fit = try(fit.variogram(vk, vm, fit.sills=T, fit.method = 7, debug.level = 1))
        plot(vk, model = vk.fit, main = 'gstat')

but i get the warning message :

  Warning message: In fit.variogram(vk, vm, fit.sills = T, fit.method =
  7, debug.level = 1) :   No convergence after 200 iterations: try
  different initial values?

and if i try to plot the vk.fit ( the exponential variogram) it doesn't appear

data:

    x          y         z
46  1.200000 -0.1000000 -1.825398
47  1.600000 -0.1000000 -2.225398
117 1.373737 -0.2244898 -2.083230
203 1.121212 -0.5918367 -1.898176
289 1.272727 -0.6734694 -2.028034
345 1.070707 -0.5510204 -1.853502
358 1.424242 -0.6734694 -2.179549
390 1.222222 -0.5510204 -2.005017
432 1.272727 -0.3061224 -2.020537
534 1.070707 -0.1428571 -1.728554
551 1.474747 -0.2244898 -2.184240
609 1.373737 -0.4693878 -2.158198
650 1.424242 -0.6326531 -2.192044
771 1.222222 -0.4693878 -2.006683


texts: Remove a legend entry if the layer is not in the composer map

texts:
I am working on a QGIS plugin in python. In this plugin I load two layers. These layers are polygons which represent administrative boundaries. My plugin load also several print layout in a Qgis project. In these print layout there is a legend. I would like that in the legend of the composer it is displayed only the layers that are visible in the map of the composer. The problem is that my administrative boundaries are represented as a polygon and therefore the layer is always present in my map while the administrative boundaries are not visible in the map if I use the option "Show features only inside the linked map". How can I do not to display the layers of the administrative boundaries in the legend if they are not visible in the extent of the map ?
I tried this code.
project = QgsProject.instance()
root = QgsProject.instance().layerTreeRoot()
groupe_admin = root.addGroup("Limites administratives")

dep_layer = QgsVectorLayer('C:/dep.shp', 'DEPARTEMENT', 'ogr')
reg_layer = QgsVectorLayer('C:/reg.shp', 'REGIONS', 'ogr')
QgsProject.instance().addMapLayer(dep_layer, False)
QgsProject.instance().addMapLayer(reg_layer, False)
groupe_admin.addLayer(dep_layer)
groupe_admin.addLayer(reg_layer)
processing.run("qgis:polygonstolines", {'INPUT': dep_layer, 'OUTPUT': 'C:/dep_lignes.shp'})
processing.run("qgis:polygonstolines", {'INPUT': reg_layer, 'OUTPUT': 'C:/reg_lignes.shp'})
depL_layer = QgsVectorLayer('C:/dep_lignes.shp', 'DEP', 'ogr')
regL_layer = QgsVectorLayer('C:/reg_lignes.shp', 'REG', 'ogr')
QgsProject.instance().addMapLayer(depL_layer, False)
QgsProject.instance().addMapLayer(regL_layer, False)

layout = QgsPrintLayout(project)
document = QDomDocument()
template_file = open('C:/Users/Thysma/Documents/Stage/Composeurs/A4 PT.qpt')
template_content = template_file.read()
template_file.close()
document.setContent(template_content)
layout.loadFromTemplate(document, QgsReadWriteContext())
project.layoutManager().addLayout(layout)
map = layout.itemById("map")

processing.run("qgis:extenttolayer", {'INPUT': map.extent(), 'OUTPUT': 'C:/emprise_carte.shp'})
layout_map_layer = QgsVectorLayer('C:/emprise_carte.shp', 'Emprise Carte', 'ogr')
processing.run("qgis:selectbylocation", {'INPUT': depL_layer,
                                         'PREDICATE': 0,
                                         'INTERSECT': layout_map_layer,
                                         'METHOD': 0})

depL_feat = depL_layer.selectedFeatureCount()
if depL_feat == 0:
    depL_layer.removeSelection()
    depL_feat = 0
    legend = [i for i in layout.items() if isinstance(i, QgsLayoutItemLegend)][0]
    legend.setAutoUpdateModel(False)
    rootGroup_admin =
    legend.model().rootGroup().findGroup('Limitesadministratives')
    rootGroup_admin.removeLayer(dep_layer)
    legend.adjustBoxSize()
    layout.refresh()


texts: I have a polygon on my map and I want to resize it in that way I add  X km to all of the vertexes. Any idea?


texts:
I try to use the build-in Python-console of QGIS to create voronoi skeleton lines from a set of polygons. Therefore I iterate over the features of my input layer (the polygons) and then I use the runalg command for v.voronoi. However, I reckon I mess things up by defining the output-file, because I keep getting the same error-message:

>>> mapLines(layer, coord_ext)
Layer failed to load!
Traceback (most recent call last):
  File "<input>", line 1, in <module>
  File "c:/users/.../tmp1qh0ii.py", line 43, in mapLines
    cent_feat = centroid_layer_obj.getFeatures()
AttributeError: 'NoneType' object has no attribute 'getFeatures'

I use the following code:

import sys
from PyQt4.QtCore import *
import processing

#Retrieve a layer
layer = processing.getObject("voetpaden_test_gebied")
#Retrieve layer extent
ext = layer.extent()
coord_ext = "%s, %s, %s, %s" %(ext.xMinimum(), ext.xMaximum(), ext.yMinimum(), ext.yMaximum())
#output folder
output_folder = r"F:\My Documents\test"

#create layer for centroid lines
line_layer = QgsVectorLayer('LineString?crs=epsg:28992', 'line', 'memory')
line_prov = line_layer.dataProvider()

#This function should create a centroid line per feature of layer
def mapLines(one_layer, coord_ext):
    for f in one_layer.getFeatures():
        #create a temporary polygon layer for each feature
        poly_layer = QgsVectorLayer('Polygon?crs=epsg:28992', 'polygon', 'memory')
        poly_prov = poly_layer.dataProvider()
        poly_prov.addFeatures([f])
        poly_layer.updateExtents()
        #place temporary layer on map canvas
        QgsMapLayerRegistry.instance().addMapLayers([poly_layer])

        #create a centroid line for the temporary layer
        centroid_line = processing.runalg('grass7:v.voronoi.skeleton', poly_layer, -1.0, coord_ext, -1.0, 0.0, 2, output_folder + "/temp.shp")
        #open the centroid line layer
        centroid_layer = iface.addVectorLayer('temp.shp', 'centroid_temp', 'ogr')
        if not centroid_layer:
            print("Layer failed to load!")

        #get the centroid layer
        centroid_layer_obj = processing.getObject("centroid_temp")
        #get the feature of the centroid layer
        cent_feat = centroid_layer_obj.getFeatures()
        line_prov.addFeatures([cent_feat])
        line_layer.updateExtents()
        #remove the temporary layer
        QgsProject.instance().removeMapLayer(poly_layer)
    #Add the new layer with centroid lines to the map canvas
    QgsMapLayerRegistry.instance().addMapLayers([line_layer])

Unfortunately, I do not see any saved object in the folder I specified (F:\My Documents\test). I figured that is why I get the NoneType AttributeError.


texts:
I cannot find the heatmap tool in QGIS 3.0.  I read that it has been moved from "Raster" to "Processing" but I cannot find it in processing either.

Update:
I do not see "Interpolation" under "processing"


texts: Why use C# rather than Java API of ArcGIS Engine (or vice versa)?

texts:
I'm going to lead a development of an ArcGIS Engine application and facing the choice mentioned in the title. We have all licenses to develop our application in C#, so Microsoft's licences are out of the equation.

I've heard (from colleagues and from the internet) that for some unknown reasons it's better to use C# instead of Java. Is that true?

The ArcGIS Engine website says that its API is available in both C# and Java. I know that ArcObjects are written in COM manner in C++, so they still need to be wrapped and it shouldn't matter whether it's in C# or Java. I also am aware of the fact that C# API has been available longer than the Java one, but does it make it more stable or preferable?

We are starting our application from scratch and it doesn't really matter for us whether it'll be written in C# or Java. Are there any reasonable reasons why one should you one or the other?


texts: Does Calculate Value run only once in Batch mode?

texts: I wrote a model with a model-only tool (Calculate Value) as the first step. While running the model in single time, it works fine. But while run in batch mode, it yields silly results. I checked the log (messages) of the batch processing and found that the Calculate Value ran only once while the rest of the model run multiple times. I confirmed Calculate Value running only once is the cause of the silly results. Does this mean the model-only tool (Calculate Value) can only run once in batch mode? Is there any way to fix it?


texts: Not all features available in OSM exports in TileMill?

texts:
I'm using TileMill and OSM Bright to style an OpenStreetMap for hiking. The problem is that OSM Bright does not import all the OSM features for styling.

I believe it's possible to add features in the osm-bright.osm2pgsql.mml file, but I'm not certain how. Does anyone know how to add, for example, path labels?


texts:
How do you manage your ArcGIS floating license pool?

I have to field several calls a week, or sometimes a day, from staff wanting to be able to check out Spatial Analyst or ArcInfo but they can't because there are none available. I've a little script I run to get a list of who's got what, and then need to either make phone calls to find out if Sam or Ellen really is still using SA or pass on the list of names to the caller so they can do the leg work. Better than half the time the user has simply forgotten to downgrade their license level, or close ArcMap/Catalog after they've done what they needed to do.

This is monkey work that really should be handled by a computer. I've not found anything that really works for us yet. What do you do?


texts: GLCM measure in an object-based classification?

texts:
I am trying to classify agricultural fields with randomForest model. I will conduct an object-based classification. To do so, in the code by Ned Horning, I provide the boundaries in a shapefile. The code extracts pixels from this area and calculates the measures such as mean of the included pixel values and std of those values. In this case the followings are my predictors:

mean_Red_Band of an object (such as a field)
mean_NDVI_Band of an object
std_Red_Band of an object
and etc.

My observations in this case are the polygons that I provided as shapefile, i.e. 3600 fields.

I have no problem with the code but the concept. My question is: can I include a unique GLCM measure of an object as a predictor in the model?

Here is a detailed explaination of my problem:

Suppose that I have an object image as this: https://imgur.com/a/OlfAzLR

When I calculated the GLCM measures for the object here is different GLCM images for different measures: https://imgur.com/a/GIKuONU

Specifically, a GLCM contrast image is as this: https://imgur.com/a/4i8aWUs

It calculates it with the following code in R:

glcm <- glcm(rasterobject,
                   window = c(9,9), ## window matrix size ##
                   shift = list(c(0,1), c(1,1), c(1,0), c(1,-1)), ## matrix direction ##
                   statistics = c("mean", "variance", "homogeneity", "contrast",
                                  "dissimilarity", "entropy", "second_moment") ##statistics##
      )

What is told about using GLCM measures is that you can use it to understand texture differences in an image. In my example, you can check the glcm_contrast image and understand the texture of the object.

However, I have multiple objects. For example, my second object and its glcm measures as follows:

https://imgur.com/a/1HZTqFu

https://imgur.com/a/gMMklLO

https://imgur.com/a/qG0QyKC

As you can see GLCM measures are calculated for different objects. If we wanted like to see this objects as a matrix, we would see a matrix like this: (It is the a part of the matrix of object2_glcm_contrast https://imgur.com/a/qG0QyKC):

A part of the matrix: https://imgur.com/a/2tzmt81

The GLCM measures are calculated and written to mid-pixel of the matrix by nrow=9 and ncol=9, thus there are multiple glcm values of an object.

My question is can I calculate a unique GLCM measure for an object? If I can so, I would be able to put it as a predictor variable of an observation, which is an object.

A unique GLCM measure could be this, for example:

glcm_contrast_mean <- sum(as.matrix(glcm$glcm_contrast),
                                na.rm = TRUE)/(length(as.matrix(glcm$glcm_contrast))
                                               - sum(is.na(as.matrix(glcm$glcm_contrast))))

And it is:

> glcm_contrast_mean
[1] 2.759012

This value is only for an object whereas, on above, there are multiple values for different parts of an object.

Could this measure (or any other) represent the object's texture?


texts: How would I draw and visualize custom maps based on OSM data?

texts:
I would like to be able to use openstreetmap data for a specific region as a source and "draw" from that the map of infrastructure improvements - mostly roads.

I want it to be self-hosted (not spoiling any data on the real OSM servers), want to have quite big flexibility and detail level in the editing process and it would be nice to render it in a pretty way.

I'm not familiar with lots of GIS software. Actually it's not that important to have OSM data as the source, it could be a satellite image as well. The important things - flexible drawing bound to real coordinates and visualization of the map.

Where could I start?

Hope this question makes sense here :)


texts:
I am new in PyQGIS and am having trouble programming how to copy features between layers.

Goal: from a robust layer (with more fields) copy its features to a 'simple' version layer.

from qgis.core import *

mapLayer = QgsProject.instance().mapLayersByName("MySQL:forexeasy,host=localhost,port=3306 authcfg='z4s291g',tables=avenida_rua|layername=avenida_rua")[0]

# selecionar estrutura dos fields name
# criar um feature
# clonar os
if (mapLayer.isValid()):
    layer = QgsVectorLayer("D:\\workstation\\qgis\\temp\\fid_62.shp", 'lay', 'ogr')

    if(layer.isValid()) :
        print('init')
        features = []
        for feature in layer.getFeatures():
            features.append(feature)

        mapLayer.startEditing()
        dataP = mapLayer.dataProvider()
        dataP.addFeatures(features)
        mapLayer.commitChanges()

Mean while I get this error:

exec(open('C:/Users/LUISMH~1/AppData/Local/Temp/tmp2vbxv1cy.py'.encode('utf-8')).read())
Traceback (most recent call last):
  File "D:\OSGeo4W\apps\Python37\lib\code.py", line 90, in runcode
    exec(code, self.locals)
  File "<input>", line 1, in <module>
  File "<string>", line 3, in <module>
IndexError: list index out of range

I don't know if I would have to create a feature for the destination layer. And pass on the attributes features from the layer to destination feature. if so how?