texts:
I saw many examples of Feature Popup activated when the cursor is moved over a feature for OpenLayers 2.0 but none for OpenLayers 3.0.

Does anyone know how to do the following : http://jorix.github.io/OL-FeaturePopups/examples/feature-popups.html on OpenLayers 3.0 ?


texts: How can I temporarily disable creation of QGIS project backup files (.qgs~) with Python?

texts:
Is it possible to temporarily disable the creation of QGIS project backup files (*.qgs~) with Python?

I'm trying to save QGIS projects in Microsoft SharePoint, where the creation of *.qgs~ files is not allowed. SharePoint has it's own file version handling, why backups are not needed.


texts: I've got QGIS 2.8.1 - Wien on my system and often use the Grass tool r.out.xyz. Two of my colleagues have recently installed QGIS 2.18.12, which shows both Grass and Grass7 processing tools. Grass7 doesn't seem to have r.out.xyz (although I thing it should acoording to the reference manual). The Grass processing tools don't work altogether, even when changing the Grass and Msys folder locations in the processing options. I've installed 2.18.12 myself (next to the 2.8.1 version), and indeed there is no Grass7 r.out.xyz tool and the Grass version does't work. In 2.8.1 the Grass r.out.xyz fortunately still works. Are we doing something wrong?


texts: PyQGIS exception: Could not load source layer for INPUT, xxx.shp not found

texts:
I'm writing a standalone script containing a series of processing tools. I appended system paths to make it consisten with QGIS consolo.

# -*- coding: utf-8 -*-
import sys
import os
sys.path.append('D:/QGIS 3.4/apps/qgis/python')
sys.path.append('D:/QGIS 3.4/apps/qgis/python/plugins')
#sys.path.append('D:QGIS 3.4/apps/qgis/python/qgis')
sys.path.append('D:/QGIS 3.4/bin')
sys.path.append('C:/Users/Administrator/AppData/Roaming/QGIS/QGIS3\\profiles\\default/python')
sys.path.append('C:/Users/Administrator/AppData/Roaming/QGIS/QGIS3\\profiles\\default/python/plugins')

os.environ['QT_QPA_PLATFORM_PLUGIN_PATH'] = 'D:/QGIS 3.4/apps/qt5/plugins/'
os.environ['QT_PLUGIN_PATH'] = r';D:/QGIS 3.4/apps/Qt5/plugins;D:/QGIS 3.4/apps/qgis/qtplugins;D:/QGIS 3.4/apps/qgis/plugins'
os.environ['PATH'] += ';D:/QGIS 3.4/apps/qgis/bin/;D:/QGIS 3.4/apps/qt5/bin/;D:/QGIS 3.4/bin/'

from qgis.core import *

import processing
from processing.core.Processing import Processing

from qgis.analysis import QgsNativeAlgorithms

prefixPath = "D:/QGIS 3.4/apps/qgis/"
QgsApplication.setPrefixPath(prefixPath, True)
qgs = QgsApplication([], False)

# load providers
qgs.initQgis()

Processing.initialize()
QgsApplication.processingRegistry().addProvider(QgsNativeAlgorithms())

areaLayer = QgsVectorLayer(r"E:\data\area_x.shp","area","ogr")
intensityLayer = QgsVectorLayer(r"E:\data\intensity.shp","intensity","ogr")

result1 = processing.run("qgis:clip",
        {'INPUT':intensityLayer,
        'OVERLAY':areaLayer,
        'OUTPUT':r'E:\data\clip1.shp'})

result2 = processing.run("qgis:clip",
        {'INPUT':r"E:\data\intensity.shp",
        'OVERLAY':r"E:\data\area_x.shp",
        'OUTPUT':r'E:\data\clip2.shp'})

qgs.exitQgis()

It works when I use a QgsVectorLayer instance as input parameter ( result1 ). However, When I change it to the layer path ( result2 ), I always get a exception:

Could not load source layer for INPUT: E:\data\intensity.shp not found.

It has troubled me for a week. I need to call a number of processing algorithm so I prefer to write commands like result['OUTPUT'] rather than QgsVectorLayer(result['OUTPUT'], 'lyrName','ogr').


texts:
I have two datasets that come from ArcGis, stored in SQL Server. I have a column called - Shape which is a geometry data type and the value looks like this:

0x110F0000010407000000DCD781517FD76F4140FAED5B2EFF46C190C2F53C9FD76F41A0703D5242FF46C1C8BAB873A1D76F416888632504FF46C1B0BFECF4A3D76F41801D38F7C5FE46C19031770984D76F41A080264AB5FE46C114D0448881D76F4100B37B6AF2FE46C1DCD781517FD76F4140FAED5B2EFF46C101000000020000000001000000FFFFFFFF0000000003

After executing this query for my Polygon dataset: SELECT Shape.STGeometryType() AS GeometryType FROM dbo.POLY

I get Polygon.

I have another dataset that is a POINT dataset which also has the same column Shape.

I use this query to do an Intersect between the two datasets:

update ABC
    set ABC.ID = mb.ID
    from STG ABC
    join POLY mb with (index([SPIDX_POLY_Shape]))
    on ABC.Shape.STIntersects(mb.Shape) = 1

My question is : How can I read this data into R, specially the Shape column, when I connect my SQL Server database. And my second question is which function will help me execute the above spatial intersect query?
I am good with R, but bad with Spatial Analysis. I am not looking for you all to write me a code, but help me understand some concepts.


texts: Setting geometry types importing to PostGIS with OGR

texts:
I'm using ogr2ogr to import data into a PostGIS database. Each dataset I need to import has (theoretically) the same exact data structure, but it appears that they are not the same geometry types. Example, I've already imported the first layer, and there's a layer called "streets". Running this command, I import the first feature class from the first database (I have GDAL configured to read .mdb):

ogr2ogr -f PostgreSQL PG:"host=localhost user=postgres dbname=database port=5432" city1.mdb streets -nln streets -nlt MULTILINESTRING -lco GEOMETRY_NAME=geometry

When I go to append (using the -append flag) the same streets layer from city2.mdb, I get this error (presumably due to conflicting geometry types for the two "streets" layers:

Warning 1: Geometry to be inserted is of type 3D Line String, whereas the layer geometry type is 3D Multi Line String.
Insertion is likely to fail
ERROR 1: INSERT command for new feature failed.
ERROR:  Geometry type (LineString) does not match column type (MultiLineString)

Is there a straightforward way with ogr2ogr to cast geometries to the same type on import? Or a way to handle this by predefining the schema and geometries for each table beforehand?


texts: Transforming from SRID 900913 to 4326 with PostGis

texts:
I have two Postgres tables with different SRIDs: 4326 and 900913. I want to convert the latter to 4326. I first changed the SRID using:

Select UpdateGeometrySRID('table', 'geomcolumn', 4326).

I then transformed the reference data using:

update table set geomcolumn = ST_Transform (geomcolumn, 4326).

This seemed to work. But when I changed the SRS setting of the relevant layer in TileMill to WGS84, I get this error:

  Detected out of bounds geographic extent.

Presumably the transformation didn't work but it's not clear why.


texts: Providing QGIS expression within a plugin

texts:
I have written a custom function in the field calculator in QGIS. How do I distribute this function to other QGIS installations the best way? I can of course instruct the user what to type or copy / paste in, but that feels kind of kludgy.

Is there any way to distribute a custom function like a, or together with or as a part of a plugin? The latter would be great as the function is to be used together with a plugin I am writing.


texts: I have a Polyline ZM shapefile. The M values are in metres. I want to change the M values to kilometres. Seems simple enough, but I can't find a tool to do this. There is only one route in the shapefile.


texts: Problem of resolution - conversion vector (points) to raster QGIS

texts:
I continue the following of the first problem :

How to convert points to a raster in QGIS (rasterize does not show z field)?

My goal is to create outlines watershed at a resolution of 5 min (10 km) :

Initially I have a data as a table (.txt) : 3 columns: the geographical coordinates (x, y) and attribute (z). (See data)
I imported the data in QGIS then I created a vector (shapefile : see map data points) and I converted the vector to a raster (see map raster test): raster> Conversion> rasterize:  vector to raster.

My problem, I do not find a better resolution for raster map.

The data (points) on my table txt or shapefile are distributed at a resolution of 5 minutes, but QGIS does not consider this resolution, I try every time to adjust "pixel size", the best resolution found is : horizontal = 14 = vertical 14. In addition I noticed that some points are outside of pixels (see raster test).

In addition I noticed that the points are not located in the middle of pixels (see raster test).

I explained, once the raster is created, I converted to vector (polygon), and when I add the original data on these polygons, I see that there are some points that are on the outline of watershed.

My goal is to make the contours of watersheds and pixilated data is located within the basin, my goal is to get this map "desired contours).

x   y   z
2.05    42.5958 1
2.13334 42.5958 1
2.05    42.6792 1
2.13334 42.6792 1
1.96667 42.7625 1
2.05    42.7625 1
2.13334 42.7625 1
2.21667 42.7625 1
1.96667 42.8458 1
2.05    42.8458 1
2.13334 42.8458 1
2.21667 42.8458 1
2.3 42.8458 1
2.13334 42.9292 1
2.21667 42.9292 1
2.3 42.9292 1
2.38334 42.9292 1
2.46667 42.9292 1
2.55    42.9292 1
2.13334 43.0125 1
2.21667 43.0125 1
2.3 43.0125 1
2.38334 43.0125 1
2.46667 43.0125 1
2.55    43.0125 1
2.63334 43.0125 1
2.71667 43.0125 1
2.05    43.0958 1
2.13334 43.0958 1
2.21667 43.0958 1
2.3 43.0958 1
2.38334 43.0958 1
2.46667 43.0958 1
2.55    43.0958 1
2.63334 43.0958 1
2.71667 43.0958 1
2.8 43.0958 1
2.88334 43.0958 1
2.05    43.1792 1
2.13334 43.1792 1
2.21667 43.1792 1
2.3 43.1792 1
2.38334 43.1792 1
2.46667 43.1792 1
2.55    43.1792 1
2.63334 43.1792 1
2.71667 43.1792 1
2.8 43.1792 1
2.88334 43.1792 1
1.96667 43.2625 1
2.05    43.2625 1
2.13334 43.2625 1
2.21667 43.2625 1
2.3 43.2625 1
2.38334 43.2625 1
2.46667 43.2625 1
2.55    43.2625 1
2.63334 43.2625 1
2.71667 43.2625 1
2.8 43.2625 1
2.88334 43.2625 1
2.96667 43.2625 1
3.05    43.2625 1
3.13334 43.2625 1
3.21667 43.2625 1
1.88334 43.3458 1
1.96667 43.3458 1
2.05    43.3458 1
2.13334 43.3458 1
2.21667 43.3458 1
2.3 43.3458 1
2.38334 43.3458 1
2.46667 43.3458 1
2.55    43.3458 1
2.63334 43.3458 1
2.71667 43.3458 1
2.8 43.3458 1
2.88334 43.3458 1
2.96667 43.3458 1
3.05    43.3458 1
2.21667 43.4292 1
2.3 43.4292 1
2.55    43.4292 1
2.63334 43.4292 1
2.71667 43.4292 1
2.8 43.4292 1
3.13334 43.3458 2
3.21667 43.3458 2
2.88334 43.4292 2
2.96667 43.4292 2
3.05    43.4292 2
3.13334 43.4292 2
2.71667 43.5125 2
2.8 43.5125 2
2.88334 43.5125 2
2.96667 43.5125 2
3.05    43.5125 2
3.13334 43.5125 2
2.88334 43.5958 2
2.96667 43.5958 2
3.05    43.5958 2
3.13334 43.5958 2
2.96667 43.6792 2
3.05    43.6792 2
3.13334 43.6792 2
3.21667 43.6792 2
3.05    43.7625 2
3.13334 43.7625 2
3.21667 43.7625 2
3.13334 43.8458 2
3.21667 43.8458 2


texts: How to change all values in a field in one operation?

texts: I am setting up a couple of shapefiles with text x, text y, text rotation and text size fields in order to utilize the labeling toolbar where I can move and rotate labels. I would like the values for all the records in the text size field to be 8. Is there a way to change them all at once, so they don't have to be done individually?


texts: Why does no line appear when using Interpolate Line in ArcMap?

texts:
I'm attempting to make an elevation profile from a DEM in ArcMap using the Interpolate Line tool.  I've done this before, and it has worked great--basically I just draw a line with Interpolate Line, then click Profile Graph to get the elevation profile.  However, for some reason, it's not working now: when I draw out my line, nothing shows up on the map.  I don't get a line, and I can't figure out what I'm doing wrong.  I've searched online for an answer, but to no avail.

Any tips?


texts: Converting probability buffers to a probability SURFACE

texts:
I have a point feature, which I have buffered using the estimated location error of each point. Each point has a different sized buffer, representing 95% confidence that the true location is withing the bounds of the buffered area.

I would like to generate a series of individual rasters (one for each of the 208 points) where the value of each cell equals 0.95/(the number of raster cells representing the buffered area). Cells representing smaller buffers have higher values.

Essentially, I think I would like to convert a multipart feature (the buffer polygons) to a SERIES of single part features (rasters representing the buffers).

Then I would like to sum the series of rasters, so that the values of each cell in the output raster represents the degree to which the buffers overlap. The more buffers that overlap, the higher the cell value will be.

Ultimately, I am trying to produce a relative probability surface, where areas where more smaller buffers overlap will have higher values than areas where fewer, or larger buffers overlap.

Any idea how I might do this?


texts: Feature Rich Mapping Service for local deployment (VM) & web deployment

texts:
Since initially posting this, I have stumbled accross solutions for some of my sub-questions. This is the updated question, with fewer requirements.

I would like to create a mapping service using OpenStreetMaps (OSM) data to meet the following requirements:

Routing

offline, for local testing & online (when deployed to the web, hosting it myself)
mode: walking, driving & bicycle
data: route geometry/shape & driving instructions / directions)
criteria: shortest, fastet, etc.

Geocoding & Reverse Geocoding

Ideally, this would all be provided by a single solution, but combining multiple 3rd party projects into a single solution would be the next best thing.

QUESTIONS:

Is there a single solution incorporating the features above?

Alternatively, can someone recommend several solutions for the various features?

Solutions I have considered so far:

GraphHopper: https://graphhopper.com/
Valhalla: https://github.com/valhalla


texts: Feature causing error in Climate linear regression analysis

texts:
In my code, for an unknown reason one of the features (in a FC with a few thousand conflict locations and their dates) is causing a problem leading to this error:
FeatureCollection (Error)
Error in map(ID=00000000000000001ae5): Number.subtract: Parameter 'left' is required.
When I print it, it does't seem to be different.
I know that it cannot be a problem with the FC because the same FC works together with an other Image Collection in a script doing the same thing.
Also when I limit the FC to exclude number 137 (ID=00000000000000001ae5) then the script works. var conflicts = ee.FeatureCollection('users/welp1992/Conflicts25_Date').limit(136)
This is the FC:
https://code.earthengine.google.com/?asset=users/welp1992/Conflicts25_Date
And this the Code:
https://code.earthengine.google.com/?accept_repo=users/welp1992/Master


texts: Number converted to dictionary in Earth Engine

texts:
When reducing an image to number on it's geographical coordinates, the code does return a number, which is what should happen.

var accum = ETo.reduce(ee.Reducer.sum());
var accumNum = accum.reduceRegion(ee.Reducer.sum(),geom,30);
accumNum = accumNum.getNumber('Potential_Evaporation_Rate_surface_6_Hour_Average_sum');
print('accum:', accumNum);

Resulting in the following print:

However, when I try to perform any other basic calculations such as an addition, it prints a weird string of code which doesn't contain the number at all.

print('accum:', 5 + accumNum);

Resulting print from previous operation:

This only concatinates the 5 to the string of code description returned by the operation.

The problem also happens when I try to set a label to the value of accumNum. However, in this case I can't display the number at all even without performing any calculations.

The only thing I need is to get the actual value from the reduced region.
Does anyone have a suggestion for this?


texts: OpenLayers performance issues adding features to vector layers

texts:
I have a fairly standard looking setup (standard in my eyes) but my performance is dreadful. My issue is that any fix stands me up for a fairly large reactor so I was hoping to get some advice.

I have a WMTS basemap, a request is made to my WFS to understand all the layers it may provide. User draws a polygon and it builds a request to my WFS server to retrieve data with a filter of this polygon, this request is fired for each of the potential 20 to 30 layers/feature types held in my WFS.
For each data type I create a new OpenLayers VectorLayer, add styling directly to the features on a per layer basis,  and add layer to the map.

The requests to the WFS all come back quickly but drawing each layer on to the map is really really poor. Depending on the data type the user can be waiting up to 20 seconds per layer. Points are quick, Multi Polygons take a long time to render.

I am dealing with a large number of features and I simply can't restrict what is visible so not expecting it to be lightning but the performance is not sufficient.

I am using OpenLayers 5, 6 was just released and when I have understood the changes I will update but looking at the release notes I am not sure it will increase performance.

Has anyone had any similar issues or success stories?


texts: Error in response to "import arcpy" - RuntimeError: NotInitialized

texts:
I am setting up a new computer and installed ArcGIS 10.4 with Python 2.7.10. I then downloaded PyScripter to use as an editor. When I try import arcpy, I get the following error:

Traceback (most recent call last):
  File "<interactive input>", line 1, in <module>
  File "C:\Program Files (x86)\ArcGIS\Desktop10.4\ArcPy\arcpy\__init__.py", line 22, in <module>
    from arcpy.geoprocessing import gp
  File "C:\Program Files (x86)\ArcGIS\Desktop10.4\ArcPy\arcpy\geoprocessing\__init__.py", line 14, in <module>
    from _base import *
  File "C:\Program Files (x86)\ArcGIS\Desktop10.4\ArcPy\arcpy\geoprocessing\_base.py", line 608, in <module>
    env = GPEnvironments(gp)
   File "C:\Program Files (x86)\ArcGIS\Desktop10.4\ArcPy\arcpy\geoprocessing\_base.py", line 605, in GPEnvironments
    return GPEnvironment(geoprocessor)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.4\ArcPy\arcpy\geoprocessing\_base.py", line 561, in __init__
    self._refresh()
  File "C:\Program Files (x86)\ArcGIS\Desktop10.4\ArcPy\arcpy\geoprocessing\_base.py", line 563, in _refresh
    envset = (set(env for env in self._gp.listEnvironments()))
RuntimeError: NotInitialized

I ensured that Python and PyScripter are both running 32 bit versions, ensured the paths were correct in PyScripter, tried setting the PYTHONPATH under the system environment settings, and I have tried suggestions found here but continue to get the same error. Any more suggestions?


texts:
I opened a polygon shape file to read the coordinate (X, Y, Z) of vertices using the following code. The coordinates of one of the triangles are show below in the picture. My question is how I get X Y Z from the PolygonZ format? I tried .zAt(0), but it works for linestring. The PolygonZ has 3 points and each point has X, Y and Z

layerT = QgsProject.instance().mapLayersByName(tPolygonLayer)[0]
featsT = [feat for feat in layerT.getFeatures()]
geom=featsT[0].geometry()


texts: joining non-spatial pollution model (csv file) to Qgis

texts:
I have csv file with concentration data but without coordinates.
Start of pollution concentration connect with that one point on map.

Is there any chance ?

Start of pollution concentration connect with that one point on map. Is there any chance ?

I am using QGIS and Rstudio.


texts: Is the OpenGeo Suite available via Docker?

texts:
Given that many GIS apps such as GeoServer & PostGIS are available as docker images, has anyone made a OpenGeo Suite Docker Image?

I'm searching for this, since the installation page for OpenGeo, says that it should not be installed on the same system as Qgis, since there will be a conflict.

This looks like a perfect use-case for Docker.

Is the OpenGeo Suite Available as a Docker Image?


texts: Making Polygon that covers the Northern hemisphere in PostGIS?

texts:
In PostGIS, I would expect ST_Geography polygons to use a clockwise winding rule, i.e. everything to the right of the edges of the polygon is inside of the polygon.

Therefore, I expected a polygon going over the equator to the west to encompass the entire northern hemisphere, and a polygon going over the equator to the east to encompass the southern hemisphere.

This is not the case, though:

select st_covers(st_geogfromtext('polygon((0 0, -90 0, -180 0, 90 0, 0 0))'),
       st_geogfromtext('point(0 -5)'))

returns true, so the polygon going to the west contains a point on the southern hemisphere.

In fact, if I change the direction of the polygon to go to the east, it also contains points on the southern hemisphere, and I can't figure out how to make a polygon that contains points on the northern hemisphere.

So where is my misunderstanding, and how do I create a polygon that contains points on the northern hemisphere?


texts: How to create new attribute on a Shapefile based on a preexisting attribute with Fiona library in PYTHON

texts:
I am initiating python with little knowledge on the theme.

Recently, I have been trying to manipulate a data-point Shapefile. Each point has attributes as date acquired and remote sensor applied. Now I have to filter this data according to these two attributes, but I am having real trouble.

I have been able to filter the sensor, but not the date. The date parameter is in string format, as the following format: "yyyy/mm/dd hh:mm:ss".

My first attempt was to split the date attribute into separate attributes like: year, month, day and so on. So later, I could filter each point through GDAL library.

Here is the code which I have have extracted from the this site and tried to apply without success:

   # pay note that elem.GetFieldAsString(7)[5:7] is supposedly the function to return the month from the date string, once it is the eighth attribute from my shapefile.

import sys
import os

try:
    from osgeo import gdal,osr,ogr
except:
    sys.exit("ERROR: GDAl library !")

import fiona
# schema: it is a simple dictionary with geometry and properties as keys
schema = {'geometry': 'LineString','properties': {'test': 'int'}}
# for defining the geometry, you need Shapely
from shapely.geometry import LineString, mapping
# two simples geometries
lines = [LineString([(272830.63,155125.73),(273770.32,155467.75)]),LineString([(273536.47,155914.07),(272033.12,152265.71)])]
with fiona.open('myshp.shp', 'w', 'ESRI Shapefile', schema) as layer:
    for line in lines:
        # filling schema
        elem = {}
        # geometry with mapping function of shapely
        elem['geometry'] = mapping(line)
        # attribute value (the same here)
        elem['properties'] = {'test': 145}
        # writing element in the file
        layer.write(elem)

# 1) Open the original shapefile:
vector_file = "original_path_file.shp"
shapefile_exit = "new_path_file.shp"

shapefile =fiona.open(vector_file)

#read the schema
schema2 = shapefile.schema
print schema2

# 2) As it is a dictionary, it is easy to add new fields/keys in the properties:

schema2['properties']['string']='oi'

# 3) Now we create a new shapefile copying myshp.shp with the new schema :

with fiona.open(shapefile, 'r') as input:
    schema = schema2
    # writing the new shapefile
    with fiona.open(shapefile_exit, 'w', 'ESRI Shapefile', schema) as output:
        for elem in input:
            # add the new attribute value

            geometry = elem.GetGeometryRef()
            elem['properties']['string']= elem.GetFieldAsString(7)[5:7]
            output.write({'properties': elem['properties'],'geometry': mapping(geometry(elem['geometry']))})


texts:
I have an app that is consuming ArcGIS REST Services through the ArcGIS JavaScript API (10.2 or 10.3).

One of the available layers contains several Polylines. These lines belong to different groups, defined through a certain attribute (e.g., groupID).

The user can click on a checkbox (among many) to select which group he wants to see. I am using setLayerDefinitions to filter the existing layer (ArcGISDynamicServiceLayer).

    var queryStatement = "groupid = " + group_id
    // group_id is provided by the selected checkbox
    var layerDefinitions = []
    layerDefinitions[23] = queryStatement
    var layer = mapServiceLayer.getLayer("layer1")
    layer.setLayerDefinitions(layerDefinitions)

This approach only allows to see one group at a time. How to see multiple groups at the same time and give them different colors in the map? I could combine multiple groups in the same queryStatement, but then they would have the same color on the map.

Another possible approach is to create a FeatureLayer on-the-fly, after each checkbox is checked. Would this scale well?


texts: Getting one image by month+year to export from MODIS NDVI using Google Earth Engine JavaScript

texts:
I previusly made a similar question to export one image with 12 bands for several years, Automate exporting one image by year for several years GEE in the previous one the images were from "MODIS/006/MCD64A1" burned area that are monthly, now I would like to do the same but with NDVI so the following code export all the images (2 by month) and I need only one by month using the median of the 2 NDVI images, I reviewed some questions about calculating NDVI by month but I can't figure out how to adapt the code here.
// Define start and end ee.Dates.
var startDate = ee.Date('2000-01-01');
var endDate = ee.Date('2021-01-01');

var col =  ee.ImageCollection("MODIS/006/MOD13A1")
             .filterDate(startDate, endDate )
             .select('NDVI');

// Define a mask to clip the data.
var countries = ee.FeatureCollection('USDOS/LSIB_SIMPLE/2017');

var country = ee.Feature(
  countries
    .filter(ee.Filter.eq('country_na', 'Australia'))
    .first()
);
var ecoregions = ee.FeatureCollection('RESOLVE/ECOREGIONS/2017');

var protectedAreas = ecoregions
  .filter(ee.Filter.and(
    ee.Filter.bounds(country.geometry()),
    ee.Filter.eq('BIOME_NUM', 4)
  ))
  .map(function(feat){
    return country.intersection(feat);
  });
// Add Uganda outline to the Map as a layer.
Map.centerObject(protectedAreas, 5);
Map.addLayer(protectedAreas);
//print('mask', mask)

// Define additional bounds.
var region = ee.Geometry.Polygon(
        [[[139.81147500927887, -22.584034738321385],
          [139.81147500927887, -39.52984428265981],
          [155.19233438427887, -39.52984428265981],
          [155.19233438427887, -22.584034738321385]]], null, false);

Map.addLayer(region);

// create years and export for every year
var years = Array.apply(null, {length: 21}).map(Number.call, Number) // sequence of 16 numbers
            .map(function(number){
              return exportImagePerYear(col, number + 2000)}); // add 2000 for each year

var scale = col.first().projection().nominalScale().getInfo();
print('Scale', scale);

// function to export
function exportImagePerYear(col, startYear) {

  var colFilt = col.filterDate(String(startYear)+'-01-01', String(startYear+1)+'-01-01');

  // Filter fire with more than 50% confidence and add a new band representing areas where confidence of fire > 50%
  var clipToRegion = function(img) {
    var dateString = ee.Date(img.get('system:time_start')).format('yyyy-MM-dd');
    var clipped = img.clip(protectedAreas);
  //  var burned = clipped.gt(0);
    return clipped.rename(dateString);
  };
  var burned_area = colFilt.map(clipToRegion);

  //var check = ee.Image(burned_area.first());
  //Map.addLayer(check, {palette: ['000000', '00FFFF'], max: 366}, 'check');
  //print('Scale', scale )

  // Stack One layer by year
  //
  var stackCollection = function(collection) {
    // Create an initial image.
    var first = ee.Image(collection.first()).select([]);

    // Write a function that appends a band to an image.
    var appendBands = function(image, previous) {
        var dateString = ee.Date(image.get('system:time_start')).format('yyyy-MM-dd');
        return ee.Image(previous).addBands(image);
    };
    return ee.Image(collection.iterate(appendBands, first));
  };
  var scale = col.first().projection().nominalScale().getInfo();

  var evi_img = stackCollection(burned_area);
  //print("EVI image stack",evi_img);

  Export.image.toDrive({
    image: evi_img,
    description: 'NDVIAustralia'+String(startYear),
    scale: scale,
    region: region,
    fileFormat: 'GeoTIFF',
    formatOptions: {
      cloudOptimized: true
    }
  });
}


texts: Changing shapefiles from a geographic coordinate system to a projected coordinate system in QGIS; converting from WGS84 to EPSG:2263

texts:
I am trying to calculate service area polygons using the QNEAT3 plug in tool for a project I am doing. However, I got an error before these polygons could be created.

2021-01-04T08:03:38     CRITICAL    The TIN-Interpolation algorithm in QGIS is designed to work with projected coordinate systems. Please use a projected coordinate system (eg. UTM zones) instead of geographic coordinate systems (eg. WGS84)!

I thought I had changed the Coordinate Reference System for all my layers to NAD83 / New York Long Island or EPSG:2263 (the project is about New York City), but after doing it again, I got the same error.
I realized that changing the CRS from WGS84 to NAD83/NY Long Island moved some layers from the accurate spot in NYC to some random spot in Pennsylvania, and on a completely different scale. However, other layers that were already in NAD83/NY Long Island were still in the right spot.
The yellow dot is one of the layers that got moved to the wrong spot after being converted from WGS84 to NAD83/NY Long Island.

I added an OSM base map to illustrate the problem.

If you zoom in, the layer is ridiculously small for some reason. A layer that should encompass an entire city is smaller than a single house.
This is what the layer is supposed to look like, and where it should be located.

Something is very wrong, but what is it?


texts: I'm looking for a feature id without going through a for or while cycle

texts:
I'm trying to find some land from a previous vector layer loaded with geojson from the geoserver.
I have a way to search:
This method goes through each field in my vector layer and compares the ids between them and where it finds results, the function operates:
imobile.getSource().forEachFeature(function(feature){

      console.log(feature.getProperties());
});

Console returns:

> e_terra: "73915" geometry: e {disposed_: false, pendingRemovals_: {…},
> dispatching_: {…}, listeners_: {…}, revision_: 2, …} id_primarie: 26
> layer: "TEREN" parcela: "3082" proprietari: "COMUNA SILIȘTEA--"
> suprafata: "429252" tarla: "145" text: 6099
> __proto__: Object

That "text" is the id of the field polygon now I would like to find that id without going through that search cycle above.


texts: Updating Field to Calculated Variable using Arcpy

texts:
I'm trying to use the updateCursor to insert a calculated value in the correct field. However, I keep on receiving an execution error. I think it is an easy fix, but I am new to arcpy and cannot figure out where I'm making a mistake.

Here's the code:

cursor = arcpy.UpdateCursor(outAgri)
for row in cursor:
    arcpy.AddMessage("Field " + str(count) + " of " + str(loop))
    agriFeat = arcpy.SelectLayerByAttribute_management(agri_lyr,       "NEW_SELECTION", '"OBJECTID_1" = ' + "%s" %count)
    count += 1

    marginbound = arcpy.SelectLayerByLocation_management (margin_lyr, "BOUNDARY_TOUCHES", agriFeat, "", "NEW_SELECTION")
    mCursor = arcpy.SearchCursor(marginbound, fields= "Margin_Type; Health_Index; Shape_Area")
    ##mFields = arcpy.ListFields(marginbound)
    mcount = 0
    gArea = 0
    pArea = 0
    rArea = 0
    gHealth = 0
    pHealth = 0
    rHealth = 0

    for row in mCursor:

        if row.getValue("Margin_Type") == "G":
            arcpy.AddMessage("Grass")
            gAreaTest = row.getValue("Shape_Area")
            gArea = (row.getValue("Shape_Area") + gArea)
            gHealth = row.getValue("Health_Index")

        if row.getValue("Margin_Type") == "P":
            arcpy.AddMessage("Planted Hedges")
            rArea = (row.getValue("Shape_Area") + rArea)
            rHealth = row.getValue("Health_Index")

        if row.getValue("Margin_Type") == "R":
            arcpy.AddMessage("Remnant Hedges")
            pArea = (row.getValue("Shape_Area") + pArea)
            pHealth = row.getValue("Health_Index")

        mcount += 1

    tArea = (pArea + gArea + rArea)

    fieldHealth = (((rArea/tArea) * rHealth) + ((gArea/tArea) * gHealth) + ((pArea/tArea) * pHealth))

    arcpy.AddMessage("The Field health is " + str(fieldHealth))
    row.setValue("Field_Health", fieldHealth)
    cursor.updateRow(row)


texts: Raster can't be clipped inside the polygon's boundaries

texts:
Here is the link what it looks like when its clipped: https://prnt.sc/r7wjfi

I just want a raster which is inside of the polygon's boundaries. I work in QGIS. Any suggestions?


texts: Why does "save" overwrite the "capture point" function?

texts:
I am very new at using QGIS. I have the standard installation (no additional plug-ins). I have several polygons and each polygon is formed by a number of segments. What I am trying to do is to connect the polygons with each other with new lines that I create. So for example I have two polygons (left and right are the two polygons) and I want to create a new segment from one point of the first polygon to another point on the second polygon.

I am using the "catch" function so it catches the points automatically. I create the segment (zoom in BEFORE saving, everything ok) BUT when I click  on "save changes" (for the segments layer) both ends of the new segment jump from their position (zoom in AFTER saving, not ok) and are not staying on the polygon-points they were captured on. Why? Also the identification chart of the new segment is empty and doesn't adopt the start and endnode. http://postimage.org/gallery/5fh42dqk/

I than tried to create a segment and type in manually the StartNode and EndNode but QGIS doesn't connect the two nodes at all. How does this work? Do I need a plugin for this?

Hope you can help me on this. Would appreciate it. Thank you, Tud

EDIT: So i narrowed it down to a round up. When i save the changes i made to the layer (so NOT the project save) QGIS rounds up the coordinates to an accuracy of six decimals. Why is that and how can i disable it? Thank you


texts: How to insert features in a MS SQL based SDE without giving ObjectID s?

texts:
I'm trying to load features into ArcSDE without using ArcObjects, and using SQL from my .NET Application instead.

I have followed the answers given in this question: Is it possible to load features into ArcSDE without using ArcObjects?

The Only problem, is that in the insert SQL Statement, I have to give the objectID, like this:

INSERT INTO [mydatabase].[dbo].[DATAPTS]
           ([OBJECTID],[tor]
           ,[value]
           ,[SHAPE])
     VALUES
           (1,dateadd(hour, 11, dateadd(dayofyear, 62 - 1, dateadd(year, 2013 - 1900, 0)))
           ,0.23654
           ,geometry::STGeomFromText('POINT (30000 198572)', 32643))
GO

Is this the best way to give the object Id?
Or is there any way to get SQL Server to auto increment this field?

(The second looks impossible, because the Featureclass created by ArcGIS in the SDE, has int as the OBJECTID, and you can't Add identity to an existing column)


texts: How to run multi level segmentation for LULC in process tree, eCognition software

texts: Can we run the eCognition  multi level (level 1 ,level 2,...) segmentation and classification in the same process tree? I want different outputs to show for each levels in a process tree. I'm in the intial learning stage of eCognition.


texts:
I'm having some trouble using the selectFeatures method with a feature layer. The code below builds a query using input text. If anyone has any pointers on what I could be doing wrong here, it'd be a huge help! I added a function to get some feedback after the selection is complete. It appears to be working fine when I test the search, but the results are not being displayed correctly.

//Query
var query = new Query();
searchButton.onclick=function(){
    var value = dom.byId('searchText').value;
    query.returnGeometry = true
    query.where = '"Type"' + " LIKE '%" + value + "%'"
    console.log('Query Registered');
    roundabouts.selectFeatures(query,FeatureLayer.SELECTION_NEW, function(selection){
            console.log('Selection Query: ' + query.where);
            console.log(selection.length + ' Results Returned');
    });
};


texts:
Is it possible to draw a line or to draw a triangle or another figure by digit the coordinates of the vertex?

I use QGIS 1.8 on windows 7 system.


texts: How do you write classify that you are using a time varaible in R?

texts:
I am using R. I am including time as a co-variate in my model.

10:30, 11:56, 12:00 etc

Is there a way to tell R that my values are time measurements?


texts:
I'm trying to inspect town assets, like pavement and signs. I want to display that information on a web-based map (CartoDB). I also want to be able to pull the most up to date information to generate reports using Crystal Reports or other similar reporting software.

My first attempt was to upload the shapefiles to CartoDB. Then I created a table with all the attributes that I would be filling out in the field that would stay synced in my DropBox. This would allow me to access the table to perform reports or update any information using excel, etc on my desktop. But I needed to merged that table with the shapefile, and if I did that, it would turn off syncing.

Does anyone know how I could get this to work? And what the best method of field collecting the data would be in order to keep everything synced with that one table in my DropBox?


texts:
We have a set of polygons which we would like to sort and number according to their latitude. This will make it easier for local authorities to relate to the numbers of the objects.

I have looked around and not found any "out-of-the-box" solutions.

Any suggestions for how to solve this using Esri tools or open source tools are interesting.


texts: Building footprint extraction in automated way

texts:
I am trying to extract building footprints automatically (even semi automated way will do) from 0.5mts optical imagery.

Currently my study area is Poland, however I would love to have a way that will give me an optimized result across the entire globe.

I have come across two potential solutions as listed below:

Using BREC4GEM software as a plugin for QGIS. However, I am not able to install the plugin in QGIS as the plugin seems to broken. Can you help with the software or the installation part?
Using Orfeo toolbox- I have the software working.  Can you help with any work flow on how to use Orfeo toolbox in extracting building footprints automatically?

Is there any other way to extract the building footprints automatically?


texts:
I have not touched Mapinfo for a good 2-3 years. (Think it was 7.5 back then, around the time it got aquired by Pitney Bowes)
I found it a great package, particularly with Vertical Mapper and utilising MapBasic to automate tedious tasks.

I would like to get an update on peoples recent experiences with the product.  I saw some nice screens of v10 where the interface looked a lot shinier.  Has it changed much over the past few years?

What about the web side of MapInfo, was it Map Extreme or something? I am pretty sure there are no similar web platforms like what ESRI offer (Silverlight/JS/etc) - Are there any plans to catch up on the web side of things?


texts:
I downloaded a shapefile from the Brazilian Institute of Geography and Statistics (IBGE). I can open it on QGIS and work with it perfectly there, but I need to work with it using Python.

I'm trying to open it with GeoPandas using the code line bellow:

shp_pr = gpd.read_file("data/PR/PR.shp")

But I get the following errors:

proj_create: Cannot find proj.db
proj_create: no database context specified
ERROR 4: Unable to open EPSG support file gcs.csv.  Try setting the GDAL_DATA environment variable to point to the directory containing EPSG csv files.
proj_create: Cannot find proj.db
proj_create: no database context specified
Unable to open EPSG support file gcs.csv.  Try setting the GDAL_DATA environment variable to point to the directory containing EPSG csv files.
proj_create: Cannot find proj.db
proj_create: no database context specified

I've also tried with GDAL, using these code lines:

from osgeo import ogr

driver = ogr.GetDriverByName("ESRI Shapefile")
ds = driver.Open("data/PR/PR.shp", 0) #open your shapefile

But I get these errors:

proj_create: Cannot find proj.db
proj_create: no database context specified
ERROR 4: Unable to open EPSG support file gcs.csv.  Try setting the GDAL_DATA environment variable to point to the directory containing EPSG csv files.
proj_create: Cannot find proj.db
proj_create: no database context specified
Unable to open EPSG support file gcs.csv.  Try setting the GDAL_DATA environment variable to point to the directory containing EPSG csv files.

I have all the ".cpg", ".dbf", ".prj", ".qpj", ".shp" and ".shx" files.

What is going wrong?


texts: Landsat footprint artifacts in statewide CCDC map

texts:
I have produced a map for New York State using CCDC as implemented in google earth engine. The map counts the number of disturbances for each pixel occurring between 1990 and 2019 (black = 0 disturbance, white = more disturbances). This map shows artifacts of the landsat scene footprints – specifically, the areas of scene overlap show more disturbances than those without scene overlap. See the striping pattern in the western part of the state. How can I mitigate this?
Disturbance count map:

Landsat scene footprints:


texts: Weird behavior using org.geotools Filter

texts:
I am using filter to catch all the features who fall into a specific polygon. My collection is in EPSG:32632 coordinate reference system and it is located somewhere in this bounding box:

xMin, yMin 470701.73,5048820.54 : xMax,yMax 641679.97,5165204.84

Filter uses polygons expressed in EPSG:4326, but usign org.geotools ReprojectingFeatureCollection this is not a problem. Here are the two polygons:

the_geom bbox POLYGON ((-180 -90, -180 90, 0 90, 0 -90, -180 -90)),
the_geom bbox POLYGON ((0 -90, 0 90, 180 90, 180 -90, 0 -90))

who are respectively reprojected into:

the_geom bbox POLYGON ((-108523640.02189268 -187176133.93183675, -108523640.02189268 187176133.93183675, 274812009.447262 187176133.93183675, 274812009.447262 -187176133.93183675, -108523640.02189268 -187176133.93183675)),
the_geom bbox POLYGON ((-69745878.65576684 -187176133.93183675, -69745878.65576684 187176133.93183675, 274812009.447262 187176133.93183675, 274812009.447262 -187176133.93183675, -69745878.65576684 -187176133.93183675))

What I expect is that applying filter with the first polygon I should get no features and when I apply filter with the second polygon I should get all the features. It's strange that filter gives back the whole features when is applied with both polygons! Here is a test to reproduce

public class FilterTesting {

@BeforeClass
public static void setUp() {
  System.setProperty("org.geotools.referencing.forceXY", "true");
}

@Test
public void filterTesting() throws Exception {
  String file = "whereisyour.shp";
  FileDataStore ds = FileDataStoreFinder.getDataStore(new File(file));
  SimpleFeatureCollection features = ds.getFeatureSource().getFeatures();
  ReprojectingFeatureCollection rfc = new ReprojectingFeatureCollection(features, org.geotools.referencing.CRS.decode("EPSG:4326"));
  Filter filter = buildFilter(-180D, -90D, 0D, 90D);
  SimpleFeatureCollection fc = rfc.subCollection(filter);
  Assert.assertTrue(fc.isEmpty());
  filter = buildFilter(0D, -90D, 180D, 90D);
  fc = rfc.subCollection(filter);
  Assert.assertFalse(fc.isEmpty());
}

private BBOX buildFilter(double x1, double y1, double x2, double y2) throws Exception {
  FilterFactory2 ff = CommonFactoryFinder.getFilterFactory2(null);
  ReferencedEnvelope bbox = new ReferencedEnvelope(x1, x2, y1, y2,
  CRS.decode("EPSG:4326"));
  return ff.bbox(ff.property("the_geom"), bbox);
}
}

I am using org.geotools version 17.2.

UPDATE

You can download a shapefile to reproduce the behavior here


texts: Options for styling Leaflet Draw line segements

texts:
I am playing around with the options for Leaflet-Draw and have came up wit a little annoyance. Using a colorpicker I can set the color for every new polyline I draw by invoking options.draw.polyline.shapeOptions.color with draw:drawstart like this:
map.on('draw:drawstart', function() {
    drawControl.options.draw.polyline.shapeOptions.color = theInput.value;
});

The problem emerges when when changing color and drawing my first polyline - the line segment still have the old color (better explained with a picture):

After finishing the firts polyline it turns to the right color and all new polylines segments have the right color but that first drawn polyline annoys me. Any idea how to solve this? There are no line segment style options in the documentation. And how do I style the dotted line?


texts: How to test spatial correlation between two point layers

texts: I have two point layers, which visually have spatial correlation. I'd like to test it but wonder if there is such a method or index？


texts: Caching tile servers for variable/high dpi - trying to understand what's what

texts:
We need to generate cacheable tiles from some map server - currently we're using geoserver using internally managed road network and other data.

It needs to allow us to get higher resolution tiles for higher resolution devices.

Clients will be openlayers and our own java desktop application(s).

What I've been told by one of my colleagues is:
WMS-C does do variable dpi but it's not really a tile caching standard and  is difficult to use.
WMS-T and a similar standard use fixed (90dpi) resolution so will look nasty on the target high resolution devices.
Proprietary systems (google maps etc) use variable dpi (or multiple pre-defined dpi?) in their apis and so look good.

Assuming the above is correct, can anyone suggest a good route forward? I can switch to something other than geoserver if necessary - though we can't afford an expensive proprietary system.

Besides serving tiles we may also want the server to serve vector data though this is "nice to have" not a requirement, as we can already serve the vector data in wfs through a simple server app we already have.


texts: Defining landuse from remote sensing data?

texts:
I am trying to research ways into obtaining landuse data for a project which calculates the visual impacts of wind parks.

I have a DEM and DOPs for a relatively large region (5000km²), but I have no landuse data.

Does anyone know how I could obtain or create such data to create a database of landuse types .i.e residential, agriculture, industry, forest, woodland etc?

Is is possible to create such a dataset from satelite imagery or would I have to think about digitizing polygons from aerial photos?

I would be very interested in hearing view or reading methodologies for similar projects.


texts: Adding .tif files with various bit depths to a raster mosaic

texts: So I have a few georeferenced .tif files that I need to add to a raster mosaic, however the problem is that some of these are 8 bit unsigned TIFFs while others are 1 bit unsigned.  After adding all the TIFFs into the mosaic, the 1 bit TIFFs display as a black square, while the 8 bit TIFFs display normally.  Is there a way to have all TIFFs display properly regardless of bit depth?  I sort of got it working by changing the image properties to display as a stretched values rather than RGB composite, but this makes all the colour TIFFs black and white which is not ideal.  Any ideas?


texts: Writing shapefile with projection defined crashes Fiona?

texts:
I have the same problem as in Writing shapefile with projection defined crashes fiona

It says that was an issue with GDAL.

I cleaned my site-packages of anything GDAL, fiona, six related. Then installed the following binaries from http://www.lfd.uci.edu/~gohlke/pythonlibs/

I using the Python27 ArcGIS10.2 interpreter

GDAL-1.11.1.win32-py2.7
Fiona-1.4.1.win32-py2.7
six-1.8.0.win32-py2.7

Tried running the code again:

with fiona.open(output_path, 'w', 'ESRI Shapefile', schema, crs=crs.from_epsg(4326)) as layer:
        for lat_long in lat_long_data:
            latitude = lat_long[0]
            longitude = lat_long[1]
            image_name = lat_long[2]
            point = self._create_shapely_point_shapefile(latitude, longitude)
            layer.write({'properties':  {'Name': image_name, 'Lat': latitude, 'Long': longitude},  'geometry': mapping(point)})
    return output_location

Python crashes and I get this error message:

Process finished with exit code -1073741795 (0xC000001D)

If I remove the crs param it works but I need it to work with the crs param.


texts: How do I get all the labels to show when they have exactly the same grid reference using QGIS?

texts: I map species records using QGIS 2.14.2. Many of the records have been submitted with exactly the same grid reference for different species. This can up to twenty records (labels).  But QGIS will only show one of the records (labels) for that grid reference. I have checked "Show all labels".  How do I get QGIS to show all the labels (records) for the one grid reference? Ideally this would be in a column or columns?


texts:
I want to simulate a pizza delivery service for San Francisco.

Because delivering Pizza to homes is taking some time (giving the pizza to clients, waiting for the money, etc.) I thought that it may be interesting to add AdditionalTime to the delivery time.

In the Network Analyst Toolbox there is the Find Routes Tool, which should be able to do so.

It says, that AdditionalTime is "The amount of time spent at the stop, which is added to the total time of the route. The units for this attribute value are specified by the Measurement Units parameter. The attribute value is included in the analysis only when the measurement units are time based. The default value is 0."

But it doesn't say how or where to enter the additional time.

Any Ideas?


texts: Google Earth Engine: Exporting many processed Images

texts:
I am new to Google Earth Engine and with the help of GIS community, I was able to learn really quick. So, kudos to the GIS community.

But I am having a bit of an issue. I am currently exporting entire MODIS Lai band image collection for US continent. My script is working, but the problem is that it's taking forever to export all the images. It's been more than a day and all my tasks are still running, I barely got 10% of my Image Collection.

Originally, I wanted to get all images (from 2002 to present) in one go but if I do that then I am getting an error from GEE. it seems GEE cannot generate that many tasks at once. The number images in modis dataset are around 1630. So, instead all the images in one go I am using filterDate to export the images in parts of 2 years.

// import the modis dataset as well as DEM or shape file
// to get the exact outline of US continent using Image.mask()
var modis = ee.ImageCollection("MODIS/006/MCD15A3H");
var dem = ee.Image("users/username/dem_us_continent");

// select the "Lai" band
var modis_lai = modis.select('Lai')
                .filterDate('2018-01-01','2020-01-01');

// Get the List of all Images in ImageCollection
var size = modis_lai.size().getInfo(); // Total no. of images (around 1630)
var modis_lai_list = modis_lai.toList(size);

// Get the geometry coordinates of dem
// Properties of dem(system:footprint)
//                      type: LinearRing
//                      coordinates: List (17 elements)
var dem_coordinates = dem.get('system:footprint');

// Iterate over all the Images to Export them in the Drive
for (var n=0; n<size; n++) {
  // Get the nth index image from ImageCollection
  var image = ee.Image(modis_lai_list.get(n));

  // Set the file name based on the Dates
  var fileName = ee.Date(image.get('system:time_start'))
              .format('YYYY-MM-DD'); // For continus date: 2002-07-185
  fileName = fileName.getInfo();

  // Masked the image with the DEM to get the US Continent Area
  var masked_image = image.mask(dem);

  // Export the image to the Drive
  Export.image.toDrive({
    image: masked_image,
    description: fileName,
    crs: 'EPSG:4326',       //Projection system of download
    scale: 1000,     //Change the scale of download to matching resolution in meters.
    folder:"GEE_Data",
    region: dem_coordinates,
    dimensions: '8930x3816',
    maxPixels: 2000000000
  });
}

// Display the last MODIS masked_image in the for loop
Map.addLayer(masked_image);

To run all the tasks at once, I have used below script in web-browser's developer mode(console).

/**
 * Batch execute GEE Export task
 *
 * Once the script ran/executed,
 * export tasks will be generated and run button will show-up in Tasks section.
 *
 * Afterwards press F12 get into console and paste the below script in it,
 * and press enter or run the below script in console.
 * All the tasks will start automatically.
 * (Firefox and Chrome are supported.)
 */

function runTaskList(){
    var tasklist = document.getElementsByClassName('task local type-EXPORT_IMAGE awaiting-user-config');
    for (var i = 0; i < tasklist.length; i++)
            tasklist[i].getElementsByClassName('run-button')[0].click();
}

function confirmAll() {
    var ok = document.getElementsByClassName('goog-buttonset-default goog-buttonset-action');
    for (var i = 0; i < ok.length; i++)
        ok[i].click();
}

runTaskList();
confirmAll();

I wanted to know, if I am doing anything wrong here or not. Because it's taking excruciating amount of time to just export all images in my drive. Thanks in advance.


texts: Create a single connected polygon from multiple polygons

texts: I am trying to create a single continuous (neighborhood) polygon from many tax block parcels. The Multipart to Singlepart GP tool does not accomplish the task since I need the polygon to be complete; i.e., without the spatial gaps between the parcels (see image). Any ideas of a quick and dirty way to accomplish this?


texts: How can I create a custom form to display the attributes of the selected object?

texts: I have already appealed to the selected objects, and found a way how to put them in a MsgBox. But I need to output was in the form


texts:
I have uploaded PostGIS/PostgreSQL tables to GeoServer and now need to create SLD style to draw points with a name of locations.

The problem is that point coordinates are in one table, and the name of the location is in another.

There is a common column in both tables, in form of PK>FK. I couldn't find a solution for taking, let say, an external label in SLD. I can't merge this two tables into one and would like to avoid using views in the database for that.


texts:
GeoHash is an excellent method to code the latitude and longitude of a point into a string. It's concise, due to the use of numbers + letters, and the precision can be tuned adding characters at the end.

I'm looking for an equivalent method to code polygons, not points. I saw two ways to do it:

geohash all polygon points, then concatenate them, and then compress the resulting string using classical text compression algorithms.
reference all geohashes included into the polygon. It's the RecursivePrefixTree method - see http://www.opensourceconnections.com/2014/04/11/indexing-polygons-in-lucene-with-accuracy/.

I would like the code could be manipulated by humans, ie. manually copied, not only by machines as an index.

Edit: see also the discussion (in French) at http://seenthis.net/messages/269838.


texts: ImportError with GeoDataFrame.from_file(..) shapefile load - Fiona issue

texts:
I've installed Geopandas and its dependencies. An ImportError results in Fiona's ogrext.so lookup of "_CPLSetThreadLocalConfigOption":

Code executed:

boros = GeoDataFrame.from_file("./data/nybb_14a_av/nybb.shp")

Error result:

Traceback (most recent call last):
  File "/Users/ronaldbjork/Documents/Boundary Solutions/geopandaone.py", line 3, in <module>
    import fiona
  File "/Applications/scipy/anaconda/lib/python2.7/site-packages/fiona/__init__.py", line 72, in <module>
    from fiona.collection import Collection, BytesCollection, vsi_path
  File "/Applications/scipy/anaconda/lib/python2.7/site- packages/fiona/collection.py", line 7, in <module>
    from fiona.ogrext import Iterator, ItemsIterator, KeysIterator
ImportError: dlopen(/Applications/scipy/anaconda/lib/python2.7/site-packages/fiona/ogrext.so, 2): **Symbol not found: _CPLSetThreadLocalConfigOption**
  **Referenced from: /Applications/scipy/anaconda/lib/python2.7/site-packages/fiona/ogrext.so
  Expected in: dynamic lookup**


texts: Natural earth tiles are not detailed enough

texts:
I am new to GIS, but I have been using the following set of tiles in some of my D3 applications:

.tiles.mapbox.com/v3/mapbox.natural-earth-2/

The issue with these tiles is that they are not too 'deep', by which I mean no tiles are displayed at higher zoom levels. Is there another natural earth set of tiles that allows for higher zoom levels?

Purely as an FYI, if needed, this is how I fetch my tiles:

image.enter().append("image")
    .attr("xlink:href", function(d) { return "http://" + "abc"[d[1] % 3] + ".tiles.mapbox.com/v3/mapbox.natural-earth-2/" + d[2] + "/" + d[0] + "/" + d[1] + ".png"; })
    .attr("x", function(d) { return d[0] * 256; })
    .attr("y", function(d) { return d[1] * 256; })
    .attr("width", 256)
    .attr("height", 256);


texts: I have thousands of records to go through and I need to pull out a specific date. Heres the set up of the fields one field is full of Null values and the second field is an ID field that looks like this XX_01052001_092045. The first part is a set of letters, the next is the date being MMDDYYYY. The third part is the time, all I want is to pull out the date and put that into the field full of null values. However the format that is requested is that the new date look like this MM/DD/YYYY. Is there a way to use the field calculator and parse this out?


texts: Cell Statistics failing to read rasters in list

texts:
Heyo. I'm trying to get the rasters from a list (about 32 rasters) and then sum them up with CellStatistics in arcpy. I'm using Python 2.7.8 and ArcGIS 10.3.1

I got it to work the first time around but then I added a new process and then the CellStatistics failed for some reason.

The code I have is:

import arcpy
import arcpy.sa
from arcpy.sa import *
from arcpy import env
import os

#Set Environment
arcpy.CheckOutExtension("Spatial")
arcpy.OverwriteOutput = True
env.extent = "MINOF"
env.cellSize = "MINOF"
Coordsystem = "GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]"

#Enter Input Data
parentdir = r'Z:\Analysis\Dave\WindRiver\P37R30\1999'
fmaskdir = parentdir + "\\" + "SavedOutputs" + "\\" + "Fmask" + "\\" + "Reclassified"
output = r'G:\trash\cloudsum\1999' #output folder
scratch = r'G:\trash\cloudsum\1999\scratch' #scratchfolder
year = 1999

fmaskfolder = fmaskdir
env.workspace = fmaskdir
fmaskrasters = arcpy.ListRasters()
fmaskrasters.sort()
print fmaskrasters

number = len(fmaskrasters)
print str(number)

env.workspace = scratch
outras = scratch + "\\" + "sum" + str(year) + ".tif"
outSUM = CellStatistics([fmaskrasters],"SUM") #the script fails here!
arcpy.CopyRaster_management(outSUM,outras,"","","","","","32_BIT_FLOAT")

factor = float(number)
finras = output + "\\" + "SUM" + str(year) + ".tif"
percentage = outSUM / factor
arcpy.CopyRaster_management(percentage,finras,"","","","","","32_BIT_FLOAT")
print finras

The error I am getting says:

      [u'1999003_Fmask_reclass.tif', u'1999019_Fmask_reclass.tif', u'1999035_Fmask_reclass.tif', u'1999051_Fmask_reclass.tif', u'1999083_Fmask_reclass.tif', u'1999131_Fmask_reclass.tif', u'1999147_Fmask_reclass.tif', u'1999163_Fmask_reclass.tif', u'1999179_Fmask_reclass.tif', u'1999187_Fmask_reclass.tif', u'1999195_Fmask_reclass.tif', u'1999203_Fmask_reclass.tif', u'1999211_Fmask_reclass.tif', u'1999219_Fmask_reclass.tif', u'1999227_Fmask_reclass.tif', u'1999235_Fmask_reclass.tif', u'1999243_Fmask_reclass.tif', u'1999251_Fmask_reclass.tif', u'1999259_Fmask_reclass.tif', u'1999267_Fmask_reclass.tif', u'1999275_Fmask_reclass.tif', u'1999283_Fmask_reclass.tif', u'1999291_Fmask_reclass.tif', u'1999299_Fmask_reclass.tif', u'1999307_Fmask_reclass.tif', u'1999315_Fmask_reclass.tif', u'1999323_Fmask_reclass.tif', u'1999331_Fmask_reclass.tif', u'1999339_Fmask_reclass.tif', u'1999347_Fmask_reclass.tif', u'1999355_Fmask_reclass.tif', u'1999363_Fmask_reclass.tif']
      32

Traceback (most recent call last):

File "C:\Users\mschauer\Desktop\ETScripts\bscripts\clouds_sum.py", line 35, in <module>

outSUM = CellStatistics([fmaskrasters],"SUM")

File "C:\ArcGIS\Desktop10.3\ArcPy\arcpy\sa\Functions.py", line 2967, in CellStatistics
ignore_nodata)

File "C:\ArcGIS\Desktop10.3\ArcPy\arcpy\sa\Utils.py", line 53, in swapper
result = wrapper(*args, **kwargs)

File "C:\ArcGIS\Desktop10.3\ArcPy\arcpy\sa\Functions.py", line 2963, in Wrapper
[function] + Utils.flattenLists(in_rasters_or_constants))

RuntimeError: ERROR 000732: Input Raster: Dataset 1999003_Fmask_reclass.tif does not exist or is not supported

What am I doing wrong?


texts:
I know this is a very simple question, however, I find it difficult to find an answer. How do I add an existing layer to my map using ArcGIS Pro SDK? I saw several examples where the layer is added while being created (for example

LayerFactory.Instance.CreateFeatureLayer(new Uri(strUri), map);

), however, how this is being done for, let's say, FeatureLayer from geodatabase?

Thanks


texts: PostGIS: Speed up intersect by filtering first

texts:
I have a table with 20 million polygons where queries take some time (with an gist index on the geometries). I thought I could speed up the queries by first filtering on an additional id (selectivity: 1/600).
But the queries take basically the same time despite also creating an additional btree_gist index on (id, geometry). The problem is that PostgreSQL runs my filter query and the intersect parallel instead of using the filter query first (Maybe it doesn't make sense? Is a selectivity of 1/600 too small? Does my (id, geometry)-index not work?)
Original query to calculate vector tiles:
WITH bounds as (...),
mvtgeom AS (
    SELECT ST_AsMVTGeom(ST_Transform(t.geom, 3857), bounds.b2d) AS geom,
    FROM data t, bounds
    WHERE ST_Intersects(t.geom, ST_Transform(bounds.geom, 32632)))
SELECT ST_AsMVT(mvtgeom.*) FROM mvtgeom;

Query with filter:
WITH bounds as (...),
WITH filter as (SELECT id FROM ... (WHERE x < n OR x > m) AND (y > n OR y < m)),
mvtgeom AS (
    SELECT ST_AsMVTGeom(ST_Transform(table.geom, 3857), bounds.b2d) AS geom,
    FROM table, bounds, filter
    WHERE ST_Intersects(table.geom, ST_Transform(bounds.geom, 32632)))
    AND table.id = filter.id
SELECT ST_AsMVT(mvtgeom.*) FROM mvtgeom;


texts:
I have been trying to use a PicLayer plugin in the JOSM editor, however after installing the plugin it would never appear in the top menu next to Windows, Audio, Help dropdowns (at least that's where I think it should be). I must say I am a total newbie to JOSM, and I do have no idea on how to even uninstall that plugin, so i could give it another try installing it back (so it hopefully shows now). I'd appreciate any hints on how to uninstall that plugin, thanks.

P.S. the reason I came here as opposed to JOSM forums was because they, for some reason have shut down the registration.


texts: Re-Opening PostGIS Tab in Log Messages Panel

texts:
I accidentally clicked the "X" next to PostGIS in the Log Messages Panel, and I cannot seem to find any way to open it again. I've checked most of the obvious settings dialogs and menus thinking this should be an easy fix, but I haven't been able to track down this particular item. I'm running QGIS v2.14.3.

Any ideas?


texts: Detection of geographical trend of a variable

texts:
What is the best method to detect and graphically display the spatial trend of a variable?
I'm looking for something that tells me the probability that the trend is from "north-west to south-east" or something like that.
Someone wrongly (my opinion) suggested me to use kriging interpolation but I think that it is not correct to infer spatial trend by using estimated variables.
I can use ArcGis or R.
Thanks


texts:
I would like to geotag my images based on a gpx file:
Preconditions:
- Images with time codes
- Gpx track points (with time codes)

Does anyone has a idea how add the images to the track points?


texts: SDE Feature Class Spatial Reference 'does not match' other layers (but it does)

texts:
I've a strange issue where I have a versioned SDE geodatabase (SQL) that contains 7 or 8 feature datasets and roughly 20 or so feature classes within them.  All of the layers have the following spatial reference:

  GCS_WGS_1984 WKID 4326 Authority: EPSG Angular Unit: Degree
  (0.0174532925199433) Prime Meridian: Greenwich (0.0) Datum: D_WGS_1984
  Spheroid: WGS_1984 Semimajor Axis: 6378137.0 Semiminor Axis:
  6356752.314245179 Inverse Flattening: 298.257223563

However, one layer in one feature dataset throws warnings when I attempt to edit the gdb, claiming that it or the other layers do not match the data frame's spatial reference (depending on what was added first.)  Meaning, this one layer just doesn't 'match' the other layers.  More so, features in this layer appear to be shifted to a certain extent.

When I look at the feature extent for this 'odd' layer IN ARC, I see:

  Max X: -67.1451770310964 Max Y: 48.7549151660984 Min X:
  -215.348767986975 Min Y: 13.3773756624525

I've applied the 'recalculate' button already.

However, I also have a map service running on this dataset and when I look at the details page for this one particular layer, I see something odd:

  Extent:
  Xmin: -180 Ymin: 13.77290001454298
  Xmax: 180 Ymax: 48.7548295050966 Spatial Reference: 4326 (4326)

The one feature class doesn't even seem to match up with the other feature classes in the same feature dataset, even though the spatial reference properties appear to be identical.

I'm not sure where to begin, especially as this is a large dataset with many relationship classes and services involved.


texts: PointPlacemark label is not hidden by the terrain

texts:
I am facing this problem where the label for the pointplacemark using the NASA worldwind API is not hidden behind terrain. The symbol hides but the label doesn't. It uses javax.media.opengl.GL and TextRenderer. If anyone has any idea about this problem, please share your views and do let me know if you need any other input.

You can see the problem in a similar thread at:
http://forum.worldwindcentral.com/showthread.php?42600-PointPlacemark-label-is-not-hidden-by-the-terrain


texts: Method for finding circle center points that cover an extent

texts:
If I have a geographic extent, is there any way to find all of the center points of all of the 5km circles that would be needed to cover the entire extent?

I'm working with an API that will return results that occur within a 5km radius of a given latitude and longitude.

I would like to come up with a way to get all of the items occurring within the extent.


texts:
I recently added Mapbox's Geocoder Control and it works wonderfully. I was wondering if anyone could tell me how I could have a marker appear at the location of the address that was searched for? Below is the code I'm currently using for the control.

.addControl(L.mapbox.geocoderControl('mapbox.places', {
    position: 'topright',
    keepOpen: false,
    autocomplete: false
}))


texts: osm2pgsql: Failed to fork helper process

texts:
I am trying to import planet.osm on am Amazon EC2 machine. It is a a machine with 16 processors and 122GB RAM.

The parameters I use are: osm2pgsql -s -H aws.rds.host -d osm -U user -C 92000 --number-processes 16 --cache-strategy dense --hstore --hstore-add-index -K -G -W -p osm --extra-attributes -v planet.osm.pbf

However, after the processing phase I get the following error:

WARNING: Failed to fork helper process 1: Trying to recover.
WARNING: Failed to fork helper process 2: Trying to recover.
...
WARNING: Failed to fork helper process 15: Trying to recover.

There is another question here with a memory issue (https://gis.stackexchange.com/questions/78879/osm-problem-with-memory), but this is not the error I get. There is no mention of "Cannot allocate memory".

Does anyone know of these issues or have any experience importing on EC2 machines?


texts: Inaccurate georeferencing when plotting raster files using rasterio

texts:
I am trying to develop a plot of a raster image with a background map. I used rasterio for reading the raster file and contextily for obtaining the background map following this tutorial.
The plotted map is wrong in terms of geocoordinates. To elaborate, the raster is not aligning with the feature it is representing (river) webmap. How can I fix this? What is going wrong?
Google colab file with a reproducible example.
pth_raster = '/content/raster_data_resampled.tif'
data, out_path, affine = reproject_raster(in_path=pth_raster, band_to_read=1, out_folder="./raster_temp", destination_crs='epsg:3857') #<-- for webmap epsg:3857
src = rasterio.open(out_path)
src.crs

>CRS.from_epsg(3857)

basemap, basemap_extent = ctx.bounds2img(*src.bounds)

import matplotlib.cm as cm
my_cmap = cm.jet
my_cmap.set_under('k', alpha=0)

plt.figure(figsize=(10,10))
plt.imshow(basemap, extent=basemap_extent)
plt.imshow(src.read(1), extent=basemap_extent, cmap=my_cmap, interpolation='none', clim=[-900, 1e5]) #<-- clim for removing no data (-9999)
plt.show()


texts: ArcGIS Conflation Tool Recommendation or Custom Script for Comparing Points

texts:
Can anyone recommend conflation tools for ArcGIS Desktop (10 or 10.1, paid or free) or a scripting outline / procedure (or just a script) to have a technician efficiently process the following rules?

We need to compare and update 2 points of interest (POI) point layers, let's call them MasterPOI and ComparePOI. Aprox 30k points in each. The rules for a technician would be the following:

Find any ComparePOI points within a ~1000m of MasterPOI.
If a ComparePOI point has a better location or a better description (using their human judgement), update MasterPOI with either location or its description field.
Separate step maybe -- then any ComparePOI point not used to update MasterPOI, simply add it to the MasterPOI. Most points would fall into this category. (I'm thinking just put a flag on a feature in ComparePOI when #2 occurs.)

I could script it myself eventually, but wondering if someone has some experience they could share


texts: What is the solution to validate if a 3D linestring is always going down? All linestring nodes must have a lower value than the previous one, to check if Z values in a line are going down continuously ... It is to validate the continuity of water lines.


texts: osm2pgsql --append super slow (node 2k/s vs 1500k/s for --create)

texts:
I'm working with multiple pbf files (one per administrative region). I thought I would import them one by one in postgis, with:

osm2pgsql --proj 4326 --hstore-all --database gis first.pbf --style custom.style --slim --flat-nodes osm2pgsql.flatnodes --create
for file in other_pbfs; do
    osm2pgsql --proj 4326 --hstore-all --database gis $file.pbf --style custom.style --slim --flat-nodes osm2pgsql.flatnodes --append
done

However, the first import (--create) is going quite fast:

  Processing: Node(262182k 1691.5k/s) Way(38130k 125.43k/s) Relation(0 0.00/s)

While the second one (--append) is going… super slow:

  Processing: Node(13k 2.5k/s) Way(0 0.00/s) Relation(0 0.00/s)

If I merge these two PBFs before import (osmium merge *.pbf -o output.pbf), it goes fast again but I'd like to remove that unnecessary step if possible?

I tried with different formats (.pbf, .osm, …) but there is no changes in speed results.


texts:
I have a set of field engineers, each with a specific expertise and I have to plan a schedule and assign customers based on their nearest reachable address based on customer requirement.

Is there an algorithm or any application where i can integrate my existing database and start the route optimization with google, bing maps or any other source ?


texts: Check if PYQT5 Widget is Visible in QGIS Plugin

texts:
I'm trying to check if my widget is alredy open/visible before showing it (or before run this function). I need to check this because in my QGIS plugin it can happen that I call this function many times. I tried isVisible() but this code returns always false. How to do this check?

def run(self):
  self.dlg = Dialog()
  if self.dlg.isVisible():
    print("is visible")
    pass
  else:
    print("is not visible")
    self.dlg.show()

  # Run the dialog event loop
  result = self.dlg.exec_()

  # See if OK was pressed
  if result:
    #Do something


texts: Exporting attributes table to csv from shp using PyQGIS 2.18?

texts:
I'm approaching now for the first time to Python console and I wanted to know step by step how to export in .csv two fields from the attribute table of a vector shapefile(.shp) using just PyQgis.
Let's say that up to now I have edited as follows:

fn = "file path"
vlayer = iface.addVectorLayer(fn, "file path", "ogr")
features = vlayer.getFeatures()
for feature in features:
    print (feature [field], feature [2])

layer = QgsVectorLayer(fn, "ogr")
QgsVectorFileWriter.writeAsVectorFormat (layer, r "C: / path of file / filename_csv", "utf-8", None, "CSV", layerOptions = ['GEOMETRY = AS_XYZ']))

From this I get a .csv file but empty, how should I proceed?


texts: ArcGIS Online Webmap: Summarize values of overlapping polygons in quantitative symbology

texts:
I am attempting to create a dynamic webmap to display varying values across a grid of polygons. Within each grid cell there are multiple overlapping polygons (grid cells) each with a unique value. I would like aggregate these values to give a total value for each grid cell which relates to the sum of each of the overlapping values. These would then be displayed using a quantitative color gradient.

I would dissolve this layer and sum each overlapping polygon's value and display it as a flat single feature per cell, however I have also added the ability for users to filter specific values out relating to a descriptive attribute within the table.

The approach I have tried was to

•Create a flat empty grid with single cells.

•Took the centroid of the each overlapping polygon.

This left me with a point layer of overlapping features (value) and a polygon layer of a grid with only an ID code.

Within the WebMap builder I have set the symbology of the grid using the Arcade expression:

var centroid = FeatureSetByName($map,"centroid_values")
IIf(Contains($feature, centroid), Sum(value), 0)

This aims to sum the "value" of each point contained within the single empty grid cell and display it. This would mean if the user was to filter the points based off an attribute the symbology would be dynamic.

However, I am running into the error:

Parse Error:featuresetbyname is not available

Is it possible to incorporate a separate featureset from the map into the symbology expression for a layer?

The featureset functions appear to be greyed out an labelled 'Not Available" within the function list.

Are there any other approaches to overcome the problem outlined above?

I have experimented with joins however the number of records means it is extremely slow to process and is not dynamic when the user applies filters.


texts: Getting information from one shapefile to another using QGIS

texts:
I am very new to GIS and learning by myself.

I am using QGIS.

I have a shapefile, which categorizes areas (land cover) by numbers (corresponding to land cover classes).

I have created a wkt file with twelve polygons, with the following structure: id|polygon

e.g.:

150-180|POLYGON((559053 4367372,559053 4364372,557553 4366970.07621135,559053 4367372))

120-150|POLYGON((557553 4366970.07621135,559053 4364372,556454.923788647 4365872,557553 4366970.07621135))

90-120|POLYGON((556454.923788647 4365872,559053 4364372,556053 4364372,556454.923788647 4365872))
....

The polygons from the wkt file superimpose categorized areas from the shapefile.

I wish to know the following:
how much of each category (from the shapefile) is included in each polygon from the wkt file?

e.g.:

150-180: 0.25 of the area is category 232 + 0.70 of the area is category 123 + 0.05 of the area is category 657

How can I do this?


texts: Using reclassified raster to make contours, contours contain the class value instead of height value

texts: I reclassified my raster to 12 classes and used the contour tool to make contours from there (they are much smoother). The issue is, the contour lines no longer contain the raster height, they contains the reclass class. I fully understand why this is happening but I was wondering if anyone has a solution or a work around so that my contour lines can still display height.


texts: Selecting point feature on Mouse click using ArcObjects?

texts:
How do I select point feature on mouse click?

I can select polygon feature but I cannot select apoint feature.

ESRI.ArcGIS.Carto.IIdentify identifyLayer = (IIdentify)iLayer;
IArray array = identifyLayer.Identify(identifyPoint);


texts: Does ArcGIS API for Python use ArcPy internally?

texts:
Today I read a GIS SE user opinion in Using Geometry.distance_to() in ArcGIS API for Python that:

  it appears that ArcGIS API for Python uses ArcPy internally

This surprised me because I do not think that this is the case.

Does the ArcGIS API for Python use ArcPy internally or are they two independent Python libraries that nevertheless may call some common underlying libraries?


texts:
I am using hibernate spatial 4.0 with hibernate core 4.3, and we are using the oracle 11g with spatial extension.

This is the entity and dao:

@Entity
public class Pts{
    @Type(type = "org.hibernate.spatial.GeometryType")
    private Geometry Shape;
    ......
}

@Repository
public class PtsDao {
    private GeometryFactory gf = new GeometryFactory(new PrecisionModel(), 8307);

    @PersistenceContext
    private EntityManager em;

    @Transactional
    public void save(Pts Pts) {
        em.persist(Pts);
    }

    public List<Pts> query(double x, double y, double radius) {
        Geometry point = gf.createPoint(new Coordinate(x, y));
        Query query = em.createQuery("select b from Pts b where distance(b.Shape,:filter) < :distance", Pts.class);
        query.setParameter("filter", point);
        query.setParameter("distance", radius);
        return query.getResultList();
    }
}

The test:

@RunWith(SpringJUnit4ClassRunner.class)
@ContextConfiguration(locations = {"/spring/application-config.xml"})
public class PtsTest {
    private GeometryFactory gf = new GeometryFactory(new PrecisionModel(), 8307);

    @Autowired
    private PtsDao dao;

    @Before
    public void setup() {
        Pts b1 = new Pts();
        b1.setShape(gf.createPoint(new Coordinate(21.401829098, 29.0819022380)));
        dao.save(b1);
        ...
        //save only five points to the db
    }
    @Test
    public void testQuery() {
        long st = System.currentTimeMillis();
        List<Pts> re = dao.query(20.82123, 30.12342, 10000);
        System.out.println("-------test distance get result " + re.size() + " " + (System.currentTimeMillis() - st));
    }
}

I found that the query will take almost 400 million seconds, I think this is unacceptable for only five points.

I think this maybe caused by the index, so I create index in oracle manually:

CREATE INDEX pts_sidx ON pts(shape) INDEXTYPE IS MDSYS.SPATIAL_INDEX

However once I run the test again, it still cost 400+ ms.

What's going on?

Since we are using the ORM framework, so the table are generated by hibernate, and this is what we got in the pl/sql developer:

-- Create table

create table PTS
(
  objectid   NUMBER(10) not null,
  address    VARCHAR2(255 CHAR),
  shape      SDO_GEOMETRY,
)
alter table PTS
  add primary key (OBJECTID)
  using index
  tablespace SPRING
  pctfree 10
  initrans 2
  maxtrans 255
  storage
  (
    initial 64K
    next 1M
    minextents 1
    maxextents unlimited
  );

And when we insert data to the database, we save nothing except the objectid(generated) and the geometry, all other column are empty.

And we do not create any query plan in the database.

When I tried to query use the native sql:

 select objectid from pts p where mdsys.sdo_geom.sdo_distance(
       p.shape,
       sdo_geometry(2001, 8307, null, sdo_elem_info_array(1,1,1),
          sdo_ordinate_array(20.82123, 30.12342)),
       0.5
) <10000 ;

It cost me 300+ ms.


texts: Identify Results form for virtual layers in QGIS keeps displaying wrong feature

texts:
When I load up a virtual layer and I use the Identify Results tool, it will display the correct information in the Identify Results panel, but in the feature form popup it will always display only the one feature (I haven't figured out which one, it isn't always the one at the top of the feature list).

In the example below, no matter which polygon I pick, the feature form will always show T1.

This is happening regardless of the complexity of the query, even if it's a simple SELECT id FROM table. And it happens whether I use the autogenerated feature form or drag and drop my own.

Can this be fixed? Is it a limitation of how virtual layers handles queries? Using QGIS 2.18.12.


texts: Calculate field coordinate in DMS (Degrees Minutes Seconds)

texts:
I have a point feature class, I want to calculate its latitude and longitude in Degrees Minutes Seconds. The feature class's spatial reference is in PCS (UTM Zone). Is there any way to calculate it using python? I've tried using code
!SHAPE.CENTROID@DECIMALDEGRESS!.split()[0] for Longitude, and
!SHAPE.CENTROID@DECIMALDEGRESS!.split()[1] for Latitude
It was worked in decimals but it was return to meters unit not in decimal degrees. Now I want to calculate it in DMS.

Any suggest?

My Code is just like this

# Add XY Coordinates
arcpy.AddXY_management(LCC_Change_Point)

# Convert Coordinate Notation
arcpy.ConvertCoordinateNotation_management(LCC_Change_Point, LCC_Change_Point_ConvertCoor, "POINT_X", "POINT_X", "UTM", "DMS", "", "GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]];-400 -400 1000000000;-100000 10000;-100000 10000;8.98315284119521E-09;0.001;0.001;IsHighPrecision")


texts:
I am new to rasterio. I'm trying to open a tif file using this code:

import rasterio
dataset= rasterio.open('/home/zaianir/Documents/DB/serie_temporelle_s2/serie_temporelle_s2/Sentinel2_T31TCJ_MASKS.tif')

But I get this error:

AttributeError: module 'rasterio' has no attribute 'open'

How do I avoid this error?


texts: How to check a polyline is self intersecting

texts:
How to check weather a single polyline is self intersecting to itself.
using arcobject C#


texts: Merging adjacent polygons to one which sum of attribute values of adjacent ones is between 1000 and 2000?

texts: I have a huge number of polygons that contain attribute field that represent number of offices inside each polygon. Now I want to merge all adjacent polygons which count (sum) of offices values is between 1000 and 2000. I want to do that using ArcGIS. I have idea to use Dissolve tool, Integrate tool and Merge tool, but I am not sure about steps in order to solve this.


texts: Keeping a lower value in some cases and a higher value in some other

texts:
Geography student here. We recently started to use Raster Calculator and I have a problem in this one case.
So we have three types of roads with a buffer zone. The road itself has value 0, a 1000-2000 unit buffer around it has value of 10 and the rest has value of 5. How can I add these 3 maps by keeping the lowest value except of the cases where value of 5 meets 10?


texts: How to return Print Layout to its default state

texts:
When I first open Print Layout it has a set of panels on the right hand side - "Item Properties", "Atlas Generation", "Composition". For better or worse, these are easily configurable, so you can delete them or add them back in different places. For me this was too easy, because now I've lost panels I need and others are put in places where it's very hard to access them.
If I just quit Print Layout and open up a new one, it seems to have inherited the layout of panels from my previous Print Layout, in which I had screwed it up completely.
Is there some way I can return the configuration of the Print Layout panels to its default state, the one I started with and which most tutorials have?


texts: BigInt data type in SQL Server mapped as Text in ArcMap

texts:
I have a table in a SQL Server 2008 database containing a BIGINT data column. When I tried to import it into ArcMap 10, it automatically converted as Text of 8 length. This caused the records that have a big integer longer than 8 fail to display. The ArcMap reported "The row contains bad value."

Could anyone help with any solution to this issue? Thanks


texts: Definite way to confirm a geocode as correct

texts:
So you tested few geocoding services for the same street address. You got results varying by 0.7 to 3 kilometers. How can you tell which geocoding service was more right, other than letting your user mark their real location, compared to where you show they are, on a map?

feels a bit like a chicken and egg problem, unless you let the user mark their location on the map.


texts:
The CORS error seen is:

  Uncaught SecurityError: Failed to execute 'toDataURL' on
  'HTMLCanvasElement': Tainted canvases may not be exported

All this doesn't work:

OpenLayers 3: Cross-Origin Request Blocked: The Same Origin Policy disallows

Enabling CORS in OpenLayers

https://stackoverflow.com/questions/20260620/openlayers-3-and-xyz-layer ( the exemple jsfiddle don't even work)

Can't do this (I don't own the server):

http://enable-cors.org/

I'm tring to save my map as PNG (http://openlayers.org/en/v3.0.0/examples/export-map.html).

The example works nice in my server with only the OSM layer, but when I add one more layer from anywhere, I got the CORS error.

My application allow the user to add arbitrary layers from any getMap link, so I have no control over the source Geoserver.

I have an addLayer function to do the job (add the user's layer to the map).

function addLayer( serverUrl, serverLayers, layerName ) {
    var newLayer = new ol.layer.Tile({
        source: new ol.source.TileWMS({
            url: serverUrl,
            isBaseLayer : false,
            crossOrigin: 'anonymous',  <<<<--- PROBLLEM
            params: {
                'layers': serverLayers,
                'format': 'image/png'
            },
            projection: ol.proj.get('EPSG:4326')
        })
    });
    newLayer.set('name', layerName);
    map.addLayer( newLayer );
    return newLayer;
}

When I add the crossOrigin: 'anonymous' nothing works anymore. Without crossOrigin it work if the map have only the OSM layer. By adding any other layer I got the CORS error.

Is there any other way to get the Map as image than using the canvas?

I've searched over all this site and found many different answers but none worked.

I cannot configure server side.

This is my export function:

function exportMap() {
    map.once('postcompose', function(event) {
        var canvas = event.context.canvas;
        var xx = canvas.toDataURL('image/png'); <<------- CORS Error !!
        console.log( xx );
    });
    map.renderSync();
}


texts: LINE PATTERN FILL for plygon symbology in qgis 2.18 by pyqgis

texts:
i would like to create this symbology from this vector layer.

by taking advatage of pyqgis and programming. i try to develope a code but i got stock , can somebody help me?
layer = QgsVectorLayer('C:\Users\ost\Desktop\ASPERG_TEST/eigentum_der_gemeinde.shp','test', 'ogr')
QgsMapLayerRegistry.instance().addMapLayers([layer])

symbol = QgsSymbolV2.defaultSymbol(layer.geometryType())
renderer = QgsSingleSymbolRendererV2(symbol)
symbol_layer = QgsSimpleFillSymbolLayerV2.create('color':'red', 'style': 'f_diagonal')
renderer.symbols()[0].changeSymbolLayer(0, symbol_layer)
layer.setRendererV2(renderer)

QgsMapLayerRegistry.instance().addMapLayers([layer])


texts: Why is vector data not saved permanently in QGIS?

texts:
I'm using QGIS 2.2.0.  I add a vector to a vector layer and save the file.  After I reopen the .qgs file, the new vector is not shown.

(1) Click 'Toggle Editing'

(2) Click 'Add Feature' and draw a polyline.

(3) Click 'Save Layer Edits'

(4) Click 'Save'


texts: How to find number of synchronized records in ArcGIS10.2.2 (Data Management)

texts:
I want to find number of records synchronized while using arcpy.SynchronizeChanges_management tool. Right now, I'm getting only messages like process is Executing, Synchronize successfully completed... and Succeeded at some time.

One way I thought of get count of all feature classes in dataset before and after the tool run and then compare the both counts. In this process I can only get the newly added but not the modified record or one record is added and other existing record is deleted from feature class then too the count will be same.

# sync parameters
#Path to source sde connection file
geodatabase_1 = source_sde
in_replica = "replica_name"
#Path to target sde connection file
geodatabase_2 = target_sde
in_direction = "FROM_GEODATABASE1_TO_2"
conflict_policy = "IN_FAVOR_OF_GDB1"
conflict_defin = "BY_OBJECT"
reconcile = "DO_NOT_RECONCILE"

# call sync tool
sync_result = arcpy.SynchronizeChanges_management(geodatabase_1, in_replica,
                                           geodatabase_2, in_direction,
                                           conflict_policy, conflict_defin,
                                           reconcile)
# want to print synchronized records
print sync_result

How can I find the number of records(newly added and updated) synchronized?


texts: Exporting webmap into PDF without layout using ArcPy with ArcGIS Pro

texts:
I try to export a webmap in using arcpy into a pdf. This works fine, if there is a Layout available using:

aprx = arcpy.mp.ConvertWebMapToArcGISProject(webmap, None, None, notes_gdb)
aprx.importDocument("C:\\path\\to\\the\\layout.pagx")
layout = aprx.listLayouts()[0]
layout.exportToPDF("C:\\temp\\export.pdf")

Is there a posibility to export the webmap without layout?

There are resolution, scale and extent given in the webmap, so it should be possible. But I have nothing found to export a map directly like what is possible in arcpy.mapping.


texts: How do you edit a stack of polygons simultaneously?

texts:
Using ArcMap 10.1.
Lets say you have 100 polygons that all are in the exact same shape and in the exact same place.  They appear to be one polygon.  How would you edit the shape of all 100 at the same time, can you share vertices among the 100 polygons?  Can you join all 100 and keep their separate attributes and Object ID's?


texts:
I want to insert a condition in a geojson filter. This condition is based in a variable called here varState:

<button id='button1'>button1</button>
<button id='button2'>button2</button>

var varState = "Initial value";

var geojsonLayer = L.geoJson(data, {
    pointToLayer: function (feature, latlng) {
      return L.circleMarker(latlng, {
        radius: 5.0,
        fillColor: '#e1118e',
        color: '#000000',
        weight: 1,
        opacity: 1.0,
        fillOpacity: 1.0
        })
      },
      filter: function(feature, layer) {
   if (varState == "01") {
     return (feature.properties.field1 == "yes" && feature.properties.field2 == "yes");
   }
   else if (varState == "02") {
     return (feature.properties.field1 == "yes" && feature.properties.field2 == "no");
    }
    else {
      return feature.properties.field1 == "yes";
    }
}

var feature_group = new L.featureGroup([]);
feature_group.addLayer(geojsonLayer);
feature_group.addTo(map);

$("#button1").click(function() {
    map.removeLayer(geojsonLayer);
    varState = "01";
    map.addLayer(geojsonLayer);
});

$("#button2").click(function() {
    map.removeLayer(geojsonLayer);
    varState = "02";
    map.addLayer(geojsonLayer);
});

It always return me the last condition despite my varState is 01,02 or null.

Does someone has already succeed this ?


texts: I need to know the distance in feet between the coordinates N34-00.000 and N34-00.001 (Degrees, minutes and thousanths of minutes) I computed approx 6 feet but am not confidant in my math skills.


texts: Error adding buffer to polyline in ArcMap 10

texts:
As the next step to what I did in this previous question, I would like to add a buffer to a polyline, with 5 meters on either side.

I tried adding buffer to two versions of the line, one in geographic WGS 1984, and the other in projected UTM 16N, but in either case I get the error "spatial reference does not match data frame".

Does ArcGIS have specific requirements as to whether I need to have the project/data frame in geographical vs. projected coordinate system? And how do I alleviate this problem? Thanks.

P.S. I am running ArcMap 10 in Windows.


texts:
(UPDATE QUESTION)

I have a problem with WFS request, mine is as follows:

var wfs = new OpenLayers.Layer.Vector("WFS",{
    strategies: [new OpenLayers.Strategy.Fixed()],
    projection: new OpenLayers.Projection("EPSG:4326"),
    protocol: new OpenLayers.Protocol.WFS({
        version: "1.1.0",
        url: "http://localhost:8080/geoserver/wfs",
        featureNS : "http://www.opengeospatial.net/cite",
        // field Namespace URI in Edit Workspace
        geometryName: "geom",
        // property of Geometry type in Edit Layer
        featureType: "InfoGeometric",
        // filed Name in Edit LayerLayer Name
        featurePrefix: 'cite'
        // field name in Edit Workspace
    }),
});

In GeoServer I can see geometric features in preview. If I go on

http://localhost:8080/geoserver/wfs/DescribeFeatureType?
    version=1.1.0&
    typename=cite:InfoGeometric

and

http://localhost:8080/geoserver/wfs/GetFeature?
    version=1.1.0&
    typename=cite:InfoGeometric

I see the columns table and data table.

After request on my map there is "WFS" layer but no features. I think that the request is submitted correctly (from my HTML page, through the proxy to localhost:8080/geoserver/wfs, code.  200 OK), but the response is not direct to my map in fact if I see it in firebug I see the proxy.cgi file while I aspect my features!!!!! Can anybody give me an advise?


texts: When I run the Join function on two vector layers the distance field keeps returning "-1"

texts: When joining two point layers by geographic location with the "Each point will be given attribute" selection, I keep getting a NEARDIST of -1.


texts:
I would like to perform the equivalent of the below ArcGIS raster calculator equation in PostGIS using the ST_MapAlgebra function and return it as a tif using ST_AsTIFF.

Con(("Sol_gain_clipped.tif" >= 1491)  &  ("can_cov_01.tif" <= 70)  & ("imp_sur_10_French_Creek.tif" <= 10),1,0)

The equation uses three single band rasters with pixel value specifications for each one, and returns a single band raster with values of 1 for pixels meeting all three requirements, and 0 for pixels that do not. Any help in how to write the SQL, including PL/pgSQL, would be greatly appreciated.


texts:
I´m trying to recover a previous version of a Qgis file, because I started playing around with CRS, and saved the file and closed it, thinking next time it woud be easy to get everything back into place.. BUT now all my layers are flying around all over the place and I can´t get them back to matching one on top of the other even when they apparently ahve the same CRS and projection (I´ve been trying to for hours, apparently there´s something seriously wrong with native projections.. but I did have them neatly lined up, so I don´t know what happened) So, I think the easiest would be to go back to last night´s version of the file. Except, when I right click the file and hit restore previous version, there are no previous versions to choose from... any ideas?
I´m running Qgis Essen on Windows Parallels (I´m on Mac Yosemite OS).


texts: Determining total number of lines connected to a line segment in QGIS?

texts: I have a line shapefile showing the street network for a major city in Africa. For each line segment, I want to calculate the total number of roads/streets connected to it. For example, the road highlighted in red is a line segment with 4 roads connected to it. How do perform this for all line segments in QGIS?


texts: Finding UTM coordinates of intersection points using PostGIS

texts:
I have two line layers, and I would like to generate a another column that shows the UTM coordinates of each intersection point.
How do I do it?
My code:
Select st_intersects(b.geom, a.geom), st_intersection(b.geom, a.geom)
From afile as a, bfile as b


texts: Two versions of ArcGIS Engine SDK on same computer?

texts:
I wish to be able to develop against both ArcGIS Engine 9.3 and 10.2, so that for example, an application could depend on either at runtime.

Are there any issues with having the SDKs for both on a developer PC?


texts:
I am trying to stack my text in an attribute table using the field calculator, but it does not seem to be working.

My field name is: Facilities.
My attributes are in this format: Public Access\Campground

I want the field to read like this:

Public Access
Campground

So I am using this expression to stack the text:

Replace( [Facilities] ,"\",vbnewline)

But when I run it, it doesn't do anything. Can anyone steer me in the right direction?


texts: How to add CADRG RPF A.TOC Raster Data Source in Geoserver 2.6

texts:
I'm currently trying to add CADRG as a datasource in Geoserver 2.6.  I've installed the necessary GDAL plugins (1.9.2) for Geoserver 2.6.

I placed my RPF folder in

C:\Program Files\Apache Software Foundation\Tomcat 7.0\webapps\geoserver\data\data\CADRG\RPF which contains the a.toc file and the folders that contain my CADRG tiles.

I'm getting the error

Could not list layers for this store, an error occurred retrieving them: null

I know the a.toc file is properly created as I've tested it in other software, and the tiles are good as I've viewed them in QGIS 2.4.

As far as I know I have all of the tiles in the correct folders, so I don't think the paths are wrong.  Does anyone have any experience loading CADRG and could help me guide me in the right direction?

Edit:

I know my GDAL plugin installation worked correctly as I have been able to add a geoTIFF store successfully, and I know my CADRG tiles are valid as I've been able to add an individual tile as an NITF store.