texts: Requesting remote data using openlayers for multiple URLS and storing respose value according to URL

texts:
I have problem creating feature count window. I have urls stored in as json object. I iterate it with jquery and send request to remote data. The I have featureStore which is initially set to 0. When the request is done i try to update it with new feature count but it only updates the last url of array here "S". The console of handler has always last value of the array here "S". Is there any way I can synchronize this? I think there are mulitple request and some request returns value faster than other without considering which url request is send. How do I synchronize so that I can update featureStore eaxctly as the elements in URL.

                var current ;
                var featuresStore =
                                    {
                                     "F" : 0,
                                     "S" : 0
                                     };

                          var urls =
                                        {
                                         "F" : "http://192.168.69.88:8080/geoserver/wfs?request=GetFeature&typename=Data:sites_2&version=1.1.1&resultType=hits&filter=(<Filter><PropertyIsEqualTo><PropertyName>Vendor</PropertyName><Literal>TE</Literal></PropertyIsEqualTo></Filter>)",
                                         "S" : "http://192.168.69.88:8080/geoserver/wfs?request=GetFeature&typename=Data:sites_3&version=1.1.1&resultType=hits&filter=(<Filter><PropertyIsEqualTo><PropertyName>Vendor</PropertyName><Literal>TE</Literal></PropertyIsEqualTo></Filter>)"
                                        };

                  //Feature Store before
                  console.log(featuresStore);

                   //reads nof of features
                  var LegendHandler = function (request) {
                          var noOfFeature;
                             try {
                                    var xmlFormat  = new OpenLayers.Format.XML();
                                    var data  = xmlFormat.read(request.responseText);
                                    var format = new OpenLayers.Format.WFST.v1_1_0({});
                                    var result = format.read(data, {output: "object"});
                                    noOfFeature = result.numberOfFeatures;

                                }catch(e){
                                    $("#messageText").val("Error "+e);
                                    console.log("Error "+e);
                                    noOfFeature = 0;
                                }
                                   featuresStore.current = noOfFeature;
                                   console.log("The current value is "+current);
                            }
            //Iterating over each URLs
             $.each(urls,function(i){
                                current = i;
                                console.log("The current value is "+current);
                        var request = OpenLayers.Request.GET({
                                    url: urls[i],
                                    callback: LegendHandler
                                });

                          });

                          //Feature Store after
                          console.log(featuresStore);


texts: Dealing with spatio-temporal data: sp::over function in R

texts:
I have a dataset called carcass, which contains occurrences of livestock carcasses in Norway in the last 5 years (each raw corresponds to a single carcass, which has an ID number, along with longitude, latitude and year). I have also a shapefile called regions, a vector file over the different regions in Norway (every region has an ID number as well, specified in the column number_reg).

I managed to plot the map of the regions in R, and to plot the points corresponding to each carcass over the map.

regions<-readShapeSpatial("\\\\homer.uit.no/fma023/Desktop/Filippo Marolla/Datasets/rbd2012_siida.shp")
plot(regions)

carcass<-readShapePoints("\\\\homer.uit.no/fma023/Desktop/Filippo Marolla/Datasets/carcass_data_edit.shp")
plot(carcass, add=T, pch=20, cex=.5, col="red")

Then, I managed to count how many points fall within each district using the over function:

res<-table(over(carcass, regions)$number_reg

I also managed to plot the temporal trend of carcass occurrences for the entire country:

plot(table(carcass$YEAR), type="o", ylab="Carcass occurrences", main="Temporal trend of carcass occurrences in Norway")

Now I need to plot carcass occurrences over time for every single region, i.e. I want a graph with number of carcasses on y-axis and time (year) on x-axis for region 1, same graph for region 2, same graph for region 3, and so forth.

However, I am in trouble since the over function in R "(...) at the spatial location of object x retrieves the indexes or attributes from spatial object y", therefore in my res object I don't have the column year included (because year is in the carcass dataset, not in the regions dataset).


texts: QGIS CLOUD expression to display multiple fields?

texts: I have a point layer with a few fields that I pushed through QGIS CLOUD. When I select the point in the cloud wep map the only field that appears is the ID. Im aware that going to the properties and selecting display allows to edit the fields displayed but I dont know what the expression would look like to display multiple fields. What would the syntax to display multiple fields look like?


texts: Splitting (stacking) labels in ArcMap in multiple lines using vbScript?

texts:
I want to show labels in ArcMap on multiple lines. The label I have is written as: Caritas Lebanon Migrant Center, DRC, Handicap International, NRC, SCI.
I want to show it in separate lines. e.g.:

Caritas Lebanon Migrant Center, DRC,
Handicap International, NRC, SCI

Or:

Caritas Lebanon Migrant Center, DRC,
Handicap International, NRC,
SCI

Is that possible to do using an expression?


texts: Values change when converting NumPyArray to raster?

texts:
I am trying to convert a numpy array to a raster that I can open in ArcGIS. When queried through numpy, the values of the array seem reasonable and correct for my purposes. An image of the array plotted through Python looks exactly like expected.

>>> numpy.min(three)
0.0
>>> numpy.max(three)
36.222499999999997
>>> numpy.mean(three)
0.017112844231034222

However, when I convert the numpy array to a raster and open it in ArcGIS, the raster contains erroneous values and will not display.

newdsm = arcpy.NumPyArrayToRaster(three, "", 0.1, 0.1)
newdsm.save("newdsm.tif")

What do I need to do in order to convert the numpy array correctly? I am using Python v2.7 and ArcGIS v10.4. It is a single-band square raster.


texts: How do can I calculate biomass for some agriculture by using VV and VH polarisation SAR images? I know there are some correlation for HV RADARSAT-2 to get biomass for some agriculture like corn, wheat, soy bean. So, what about sentinel-1?


texts:
I'm having some difficulty when using a Logical OR in a WFS query to get multiple postcodes. This code works fine when it's only 1 PropertyIsEqualTo filter, but when I wrap the OR around more than one PropertyIsEqualTo filter I get an exception.

POST http://intranet.nottscc.gov.uk/maps/wfs HTTP/1.1
User-Agent: Fiddler
Content-Type: text/xml
Host: intranet.nottscc.gov.uk
Content-Length: 1875
Expect: 100-continue

<wfs:GetFeature service="WFS" version="1.1.0" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns:wfs="http://www.opengis.net/wfs" xmlns:ogc="http://www.opengis.net/ogc" xmlns:gml="http://www.opengis.net/gml" xmlns:ncc="http://nottinghamshire.gov.uk" xsi:schemaLocation="http://www.opengis.net/wfs http://schemas.opengis.net/wfs/1.1.0/wfs.xsd">
<wfs:Query typeName="ncc:OSCodepoint_15KM">
<wfs:PropertyName>ncc:FORMAT1</wfs:PropertyName>
<wfs:PropertyName>ncc:EASTING</wfs:PropertyName>
<wfs:PropertyName>ncc:NORTHING</wfs:PropertyName>
<ogc:Filter>
  <ogc:OR>
    <ogc:PropertyIsEqualTo>
      <ogc:PropertyName>ncc:FORMAT3</ogc:PropertyName>
      <ogc:Literal>NG156LR</ogc:Literal>
    </ogc:PropertyIsEqualTo>
    <ogc:PropertyIsEqualTo>
      <ogc:PropertyName>ncc:FORMAT3</ogc:PropertyName>
      <ogc:Literal>NG163GP</ogc:Literal>
    </ogc:PropertyIsEqualTo>
    <ogc:PropertyIsEqualTo>
      <ogc:PropertyName>ncc:FORMAT3</ogc:PropertyName>
      <ogc:Literal>NG158AY</ogc:Literal>
    </ogc:PropertyIsEqualTo>
    <ogc:PropertyIsEqualTo>
      <ogc:PropertyName>ncc:FORMAT3</ogc:PropertyName>
      <ogc:Literal>NG171ER</ogc:Literal>
    </ogc:PropertyIsEqualTo>
      </ogc:OR>
    </ogc:Filter>
  </wfs:Query>
</wfs:GetFeature>

This is the Error message I get.

<?xml version="1.0" encoding="UTF-8"?>
<ows:ExceptionReport version="1.0.0" xsi:schemaLocation="http://www.opengis.net/ows http://intranet.nottscc.gov.uk/maps/schemas/ows/1.0.0/owsExceptionReport.xsd"  xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns:ows="http://www.opengis.net/ows">
  <ows:Exception exceptionCode="NoApplicableCode">
    <ows:ExceptionText>Error occurred getting features Invalid fid filter provides, has no fids inside Invalid fid filter provides, has no fids inside</ows:ExceptionText>
  </ows:Exception>
</ows:ExceptionReport>

I'm using GeoServer to process the query.

Any help is much appreciated.


texts:
I'm wondering if it's possible to move the position of legend when using CartoDB.js.

I can see in the layer api I can turn it on or off, but ideally I'd like to be able to specify a div for it to be nested within. Is this possible?


texts: ArcGIS Javascript Spatial Query with Lat Lng or X Y?

texts:
I am using Google Maps on the front end as all I need on this view is to show the points. When user clicks on the map in the view, it selects a point hence giving out lat and lng. I am trying to use the lat and lng to query my feature layers and to find the features containing the point but it is giving the same error again and again.

TypeError: undefined is not a function {stack: (...), message: "undefined is not a function"} "TypeError: undefined is not a function
at p.toJson (http://js.arcgis.com/3.13/init.js:589:112)
at p.execute (http://js.arcgis.com/3.13/init.js:267:406)
at A._createWrappers.p.forEach.c.(anonymous function).g.addCallbacks.d (http://js.arcgis.com/3.13/init.js:723:454)
at c (http://js.arcgis.com/3.13/init.js:76:221)
at b.Deferred.A.then.then [as then] (http://js.arcgis.com/3.13/init.js:78:301)
at b.Deferred.addCallbacks (http://js.arcgis.com/3.13/init.js:78:144)
at A._createWrappers.p.forEach.c.(anonymous function) [as execute] (http://js.arcgis.com/3.13/init.js:723:382)
at http://localhost:26196/Projects/Add:4678:27
at ha (http://js.arcgis.com/3.13/init.js:32:473)
at http://js.arcgis.com/3.13/init.js:25:61"

The code that is throwing the error is pasted below.

require([
          "esri/map", "esri/tasks/query", "esri/tasks/QueryTask", "esri/SpatialReference", "esri/geometry/Point", "dojo/dom"
        ], function (Map, Query, QueryTask, SpatialReference, Point, dom) {

            var query = new Query();
            query.outSpatialReference = { wkid: 4326 };
            query.spatialRelationship = Query.SPATIAL_REL_INTERSECTS;
            query.returnGeometry = true;
            query.outFields = "*";
            var point = new Point();
            point.setX(lng);
            point.setY(lat);
            query.geometry = point;
            query.where = "1=1";
            var queryTask = new QueryTask("http://myServer/arcgis/rest/services/ApplicationServices/DataService/MapServer/6/");
            queryTask.execute(query, function(a){ debugger }, function(e){ debugger });
        });

I found the snippet somewhere on line, the only difference is that query.geometry was set as event.Point in the example while I am setting a point myself.


texts:
Please can anyone have experience with netcdf data with irregular grid.
my data file is taken from NCEP ADP Global Upper Air and Surface Weather Observations (PREPBUFR format)

So i need to extract data in a ASCII file like this :

STATION_ID   LAT  LON     OBS_VAL ...
  70585     45.5  60.5    1005.5    .....

the header information is :

[hachelaf@skynet1 test]$ ncdump -h gdas.44512.2010040100.nc
netcdf gdas.44512.2010040100 {
dimensions:
        mxstr = 16 ;
        hdr_arr_len = 3 ;
        obs_arr_len = 5 ;
        nobs = UNLIMITED ; // (98012 currently)
        nhdr = 30766 ;
variables:
        char obs_qty(nobs, mxstr) ;
                obs_qty:long_name = "quality flag" ;
        float obs_arr(nobs, obs_arr_len) ;
                obs_arr:long_name = "array of observation values" ;
                obs_arr:_fill_value = -9999.f ;
                obs_arr:columns = "hdr_id gc lvl hgt ob" ;
                obs_arr:hdr_id_long_name = "index of matching header data" ;
                obs_arr:gc_long_name = "grib code corresponding to the observation type" ;
                obs_arr:lvl_long_name = "pressure level (hPa) or accumulation interval (sec)" ;
                obs_arr:hgt_long_name = "height in meters above sea level (msl)" ;
                obs_arr:ob_long_name = "observation value" ;
        char hdr_typ(nhdr, mxstr) ;
                hdr_typ:long_name = "message type" ;
        char hdr_sid(nhdr, mxstr) ;
                hdr_sid:long_name = "station identification" ;
        char hdr_vld(nhdr, mxstr) ;
                hdr_vld:long_name = "valid time" ;
                hdr_vld:units = "YYYYMMDD_HHMMSS" ;
        float hdr_arr(nhdr, hdr_arr_len) ;
                hdr_arr:long_name = "array of observation station header values" ;
                hdr_arr:_fill_value = -9999.f ;
                hdr_arr:columns = "lat lon elv" ;
                hdr_arr:lat_long_name = "latitude" ;
                hdr_arr:lat_units = "degrees_north" ;
                hdr_arr:lon_long_name = "longitude" ;
                hdr_arr:lon_units = "degrees_east" ;
                hdr_arr:elv_long_name = "elevation" ;
                hdr_arr:elv_units = "meters above sea level (msl)" ;

// global attributes:
                :FileOrigins = "File /glade/p/rda/transfer/dsrqst/RABAH44512/gdas.44512.2010040100.nc generated 20130912_180331 UTC on host geyser01 by the MET pb2nc tool" ;
                :MET_version = "V4.1" ;
                :MET_tool = "pb2nc"


texts: Error: Unable to open image while ingesting into earth engine

texts:
I am trying to upload and ingest a raster file of TIF format into the earth engine. The data is stored in GCS and I used CLI to make the upload. Upload is complete, but ingestion is interrupted throwing an error:

Error: Unable to open sample_image.tif

I am using jupyter notebook to implement the code and here is the code to upload data
!earthengine upload image --asset_id=users/ksvsp/sample_image pyramiding_policy=sample gs://cloud-ksvsp-ee/sample_image.tif`

What am I doing wrong?


texts: How to get Z attribute of point using Python and ArcGIS Desktop 9.3?

texts:
I am making a python code for arcgis 9.3. In this code i create a new point with x,y and z attributes as following :

rowsLine = gp.updatecursor(line_Feature_Class)
rowLine = rowsLine.next()

rowsPoint = gp.InsertCursor(point_Feature_Class)
feat = rowsPoint.NewRow()

while rowLine:

    #...

    ptObj = gp.CreateObject('Point')
    ptObj.x = xPremierPoint
    ptObj.y = yPremierPoint
    ptObj.z = zPremierPoint
    feat.Shape = ptObj
    rowsPoint.InsertRow(feat)

I want to add the z value in a field "zvalue" in my "point_Feature_Class". The problem is  i can't use an insert cursor (rowsPoint) and an update cursor at the same time.

So, at the end of the code, I want to delete the insert cursor, create an upsdate cursor, to loop in my "point_Feature_Class" and get the z attribute of each point and add it in the fiels "zvalue".

The code should be something like that :

rows = gp.updatecursor(point_Feature_Class)
row = rows.next()

While row :

    Zvalue = #???? row.shape.getz ????
    row.SetValue("zvalue", Zvalue)

rows.updateRow(row)
row = rows.Next()

Thanks for your help ;)


texts: ArcGIS Online Symbology - Match Custom Used Symbol

texts:
Created a WebMap and published it through ArcPro about a year ago. It is now hosted on AGOL, I am now updating the layers of the map and adding new ones.
When clicking on the "Change Style" for the layer, I can see my symbol (which is not a standard symbol) at the top, but cannot select it to use anywhere else. I don't remember if I used an image (I've searched my HD and my online content, not found), and I've searched all the standard symbols. Is there a way to use the same symbology from an existing layer to a newly uploaded layer?
Edit: Upon review, the answer below from Pfalbaum points out that these can only be referenced online (there is no upload from HD option), and looking at my symbol it does not show the web address I used for the symbol.
I'm afraid there may not be an answer.


texts:
Statistics (mean, min, max) change quite a lot after importing a .TIF grid file (DEM) into GRASS. According to GDAL queries in QGIS, the TIF has the following properties:

STATISTICS_MAXIMUM=2.2400000095367
STATISTICS_MEAN=0.05137734831008
STATISTICS_MINIMUM=-3.5300002098083
STATISTICS_STDDEV=0.43984890430363

After import in GRASS, with region set according to this grid (extent & cellsize), GRASS lists the properties as follows:

...
minimum: -4.21
maximum: 2.99
...
mean: 0.0451878
standard deviation: 0.435221

Especially the difference in the 'mean' worries me a bit.  I checked extent and cell size, they are identical.

Does anybody have any idea why the results (especially the 'mean') are different?

I'm new to this forum, and pretty much a beginner with GRASS (7.0.0beta3 installed atm).


texts:
My software produce WKT and I'm trying to import it into an ODBC connected Geopackage:

SELECT load_extension('mod_spatialite');

SELECT IsValidGPB(
    AsGPB (
        ST_GeomFromText('POLYGON((  318830.663 301484.952 18.120,
                                    318821.415 301490.594 18.405,
                                    318814.969 301480.037 18.489,
                                    318824.223 301474.375 18.132,
                                    318830.663 301484.952 18.120))'
                        , 3163
        )
    )
) AS "IsValid?";

Result: 0

Is there a trick to make a valid GPB geom?


texts: Postgres DB updates visible to Mapinfo Application

texts:
I design a MapInfo 10.5 app that connects to a Remote Linux Server hosting a PostgreSQL DB. When starting the app I'm reading some tables from DB and create a local version of them that I work on it (query info, display thematic maps etc.). I'm using commit/refresh procedure to update info on DB. On the server hosting the DB I'm using another app (develop with Java) that updates the very same tables I'm using locally on MapInfo. Is there a method to be inform in MapInfo app when these changes take place on the server?

ex. I read in MapInfo table GSMCells from PostreSQL and working on it locally query or display some info.

"Server hdbc Link Table "SELECT  *  FROM  ""public"".""GSMCells""" Into GSMCells    File "C:\Temp\GSMCells"

Meanwhile a new version of GSMCells is available on the server. Me, as a user, I'm not aware of these change unless I'm doing a refresh:

"Server Refresh GSMCells"

Since I don't know when the new version is available and not willing to refresh periodically I will like to be informed when this is happening. Any idea how?


texts: What output colour format is expected by the QGIS data defined properties widget?

texts:
I'm trying to write a python function that takes a BGSREF attribute from a BGS geology shapefile and uses it to colour up the geology polygons of a shapefile.

I posted my non-code working as a part answer to Data-defined Styles in QGIS as it shows the maths behind taking a BGSREF code and converting it to RGB values, but I would like to know what needs to be corrected to have the full solution.

Using Nathan's QGIS Blog post on User defined expression functions for QGIS as a basis I created my code like:

from qgis.utils import qgsfunction
from qgis.core import QGis

aColourAdjust = [0,7,14,21,31,42,54,67,80,100]

@qgsfunction(1, "Python")
def BGSrefToQGIScolorRGB(values, feature, parent):
    """
       Converts BGSref colour (attribution) to RGB values.
       Output as integers in the range 0 - 255, in the style:
       color_rgb(rrr,ggg,bbb) or color_rgba(rrr,ggg,bbb,aaa)
    """
    # set some variables
    #rgbValue = "color_rgb(255,255,255)"
    rgbValue = "255,255,255"
    #rgbValue = "<255>,<255>,<255>"
    isTransparent = 'False'

    bgsref2Convert = values[0]

    # Test the length of the input value and pad with zeros if applicable.
    if len(str(bgsref2Convert)) == 1:
        bgsref2Convert = "00" + str(bgsref2Convert)
    elif len(str(bgsref2Convert)) == 2:
        bgsref2Convert = "0" + str(bgsref2Convert)
    elif len(str(bgsref2Convert)) == 3:
        pass
    else:
        isTransparent = 'True'

    #Now we can calculate the RGB values...
    if isTransparent:
        #rgbValue = "color_rgba(255,255,255,255)"
        rgbValue = "255,255,255,255"
        #rgbValue = ""<255>,<255>,<255>,<255>"
    else:
        yellow = aColourAdjust[int(bgsref2Convert[:1])]
        cyan = aColourAdjust[int(bgsref2Convert[-1:])]
        magenta = aColourAdjust[int(bgsref2Convert[1:2])]

        # 100.0 is deliberate to get a floating point result, otherwise the result of division will always be zero.
        red = 255-((cyan/100.0)*255)
        green = 255-((magenta/100.0)*255)
        blue = 255-((yellow/100.0)*255)

        #Now just need to find the expected QGIS format for the colour output...

        #rgbValue = "color_rgb(" + str(int(red)) + "," + str(int(green)) + "," + str(int(blue)) + ")"
        rgbValue = str(int(red)) + "," + str(int(green)) + "," + str(int(blue))
        #rgbValue = "<" + str(int(red)) + ">,<" + str(int(green)) + ">,<" + str(int(blue)) + ">"
    return rgbValue

I start QGIS (2.2.0), open the Python Console, import my function, and it shows in the expression string builder.

For my chosen layer style properties, I choose Single Symbol, data defined properties, Expression and create an expression like:  BGSrefToQGIScolorRGB(  "BGSREF" ) and apply this, but I never seem to get any colours in the map.

The help tip tells me the colour format is: '<red>,<green>,<blue>,<alpha>' but what does this actually mean?

I've tried outputting the results as the following strings (255 here indicates number format not the literal output):

* <255><255><255><255>
* <255><255><255>
* 255,255,255
* color_rgb(255,255,255)

Nothing seems to work for me.

What format is expected?


texts: How do I create a 'heatmap' in QGIS based on the values of the different points?  I did do this once before with rainfall, but can't remember how!


texts:
I am using Using QGIS 18.16 to try and create a grid measured in feet.
I am having some issues with the best way to create a accurate grid in a project projected in WGS 84.  I know that it uses degrees but i need to come up with a way to create it in feet.  I have tried the Vector Grid creation and a few plugins but i feel like i am missing something and making this a little to hard.


texts: Conversion of array to tif_Expected array of dim 2

texts:
I'm trying to perform a simple operation on an image, which is applying a median filter to it (for now). I'm trying to do it with QGIS and the console Python. My code is :

import numpy as np
from scipy import ndimage
import gdal
from osgeo import gdal
from osgeo.gdalconst import *

#Image reading
rasterpath = r"//sites.fr/DVR/RHONE-ALPES/RCA/Traitement_PTHR3D/PFE/Caroline GRUSSI - 2020/Données Test/Tunnel Croix Rousse - Lyon/2018_gris.tif"
ds = gdal.Open(rasterpath)
band1 = ds.GetRasterBand(1)

#computation
im_med = ndimage.median_filter(band1, 3)

#ouput image creation
cols = ds.RasterXSize
rows = ds.RasterYSize
in_band = ds.GetRasterBand(1)

driver = gdal.GetDriverByName('GTiff')
outRaster = driver.Create("//sites.fr/DVR/RHONE-ALPES/RCA/Traitement_PTHR3D/PFE/Caroline GRUSSI - 2020/Données Test/Tunnel Croix Rousse - Lyon/2018_mediane.tif", cols, rows, 1, in_band.DataType)
outRaster.SetGeoTransform(ds.GetGeoTransform())
outRaster.SetProjection(ds.GetProjection())
outband = outRaster.GetRasterBand(1)
outband.WriteArray(im_med)

The error message i'm getting is "ValueError : expected array of dim 2".
So I checked, and wrote

im_med.shape

To which the response is " () ". So dimension zero ?
Does anyone has an idea about what I am doing wrong ?


texts: How to Import the geometry/feature collection which is in GeoJSON format into Google Earth Engine?

texts:
I have a GeoJSON file. I have opened it in Sublime Text Editor

How can I import it in Google Earth Engine?


texts:
I am teaching myself GIS in baby steps, I recently acquired aerial imagery (employing RGB and NIR-GB sensors) using a small UAV with the intent of generating a map in a small plot of land. After stitching the images, the resulting georeferenced orthomosaics came out misaligned due to the differing GSD (the sensors aren't identical) or slightly different image geotags between flights.

This picture shows the misalignment in question after applying some transparency.

Is it possible to align and resample both of the layers in QGIS? If so, what procedure can I use to accomplish this task?


texts: Is there some custom project existing of geoexplorer free downloadable? I'm student and I'm doing a stage. I need to have a query mask for my map but I'm not able to do it with sdk-suite.


texts: Inserting LiDAR points (from laspy) in GeoDataFrame without using a numpy array?

texts:
Working with Python 3.6 (Anaconda) under Windows using a laspy fork that works with Python 3 (https://github.com/sethrh/laspy).

I am loading a LAS file to a GeoPandas GeoDataFrame to eventually perform a spatial join (using sjoin) with a polygon layer. I am sure there is a more elegant way to do what I am doing using the kludge below:

import numpy as np
from laspy.file import File
from pandas import DataFrame
from geopandas import GeoDataFrame
from shapely.geometry import Point

#Read LAS file
inFile = File("s428_7568.las", mode = "r")

#Import LAS into numpy array (X=raw integer value x=scaled float value)
lidar_points = np.array((inFile.x,inFile.y,inFile.z,inFile.intensity,
               inFile.raw_classification,inFile.scan_angle_rank)).transpose()

#Transform to pandas DataFrame
lidar_df=DataFrame(lidar_points)

#Transform to geopandas GeoDataFrame
crs = None
geometry = [Point(xyz) for xyz in zip(inFile.x,inFile.y,inFile.z)]
lidar_geodf = GeoDataFrame(lidar_df, crs=crs, geometry=geometry)
lidar_geodf.crs = {'init' :'epsg:2959'} # set correct spatial reference

So, is there a better way, more efficient, more pythonesque way of injecting the .las file points (from laspy) to the GeoPandas dataframe without passing through a numpy array?


texts: ArcGIS Server Manager requesting log on for every page of web interface, how can I prevent this from happening?

texts:
The first issue that came up was that the Certificate was not recognized. I saved this to my machine and imported it to all tabs of the Certificate Authority. I also realized that the link being used didn't have the FQDN. When I entered this it seemed to correct the Initial Certificate Error which was thought to be the issue.

As it turns out this wasn't the case. When I log in with the Server Administrator account and I try to go to, for example, the Security Link. It will send me immediately to a log in page and I would have to enter the credentials again to proceed. I was also informed by the Senior Admin who passed on this troubleshooting issue to me that the server was changed and the new license was placed on the machine and this was not the issue before.

How can I rectify this constant having to log in every time a feature is accessed from the Web Interface?


texts: Combining DEM, GeoTiff and CityGML in QGIS

texts:
I'm trying to combine DEM, GeoTiff, and CityGML in QGIS to create a 3D representation like in this example.
The problem I have is that:-

The CityGML buildings are hovering above the 3D view and

The GeoTiff doesn't take the shape of the CityGML buildings like it does with the DEM. Here is what it looks like:

I'm working with open data. They offer the LoD2 data in the CityGML, DXF, and Shape format. In QGIS the CityGML buildings are just displayed as 2D contour lines, whereas the DXF and Shapefiles are properly displayed as 3D. Is the data corrupted? For reference here is the data I used:
DEM, CityGML, DXF, Shape, GeoTiff


texts:
I've created a polygon from lines in a couple of cases where I'm unable to select the polygon using the standard Select Features button. I am able to select the feature using the Attributes Table, but in neither case am I able to actually cut and paste the polygon from that layer to another layer.

Any ideas?

(Using QGIS 2.18 on OSX)


texts: Loop exits after 100 count while creating domains using python

texts:
I am trying to create domains on feature class or table class by reading a table which holds all look up information. Following is example code for what I am doing:

tableIndex = tableClasses.index("LookupTable")
tableGnLookup = tableClasses[tableIndex]

with arcpy.da.SearchCursor(tableGnLookup,
                           ["DATA_TABLE", "DATA_REF", "LOOKUP_TABLE", "LOOKUP_REF", "LOOKUP_DESCRIPTION"],
                           where_clause="DATA_TABLE IS NOT NULL AND LOOKUP_TABLE IS NOT NULL",
                           sql_clause=["", "ORDER BY LOOKUP_TABLE, DATA_TABLE"])  as cursor:

    # Itrate each row
    for row in cursor:

        # Get row values
        dataTableName = str(row[0])
        dataRefField = str(row[1])
        lookupTableName = str(row[2])
        lookupRefField = str(row[3])
        lookupDescriptionField = str(row[4])

        # Find table in all tables
        try:
            tableIndex = tableClasses.index(lookupTableName)
        except:
            tableIndex = -1

        # Get table if table index is found
        if tableIndex >= 0:
            lookupTable = tableClasses[tableIndex]
            arcpy.TableToDomain_management(lookupTable, lookupRefField, lookupDescriptionField, database, lookupTableName, lookupTableName)

        # Find feature class in feature class list
        try:
            tableIndex = featureClasses.index(dataTableName)
        except:
            tableIndex = -1

        # Get feature class at found index
        if tableIndex >= 0:
            fc = featureClasses[tableIndex]
            try:
                arcpy.RemoveDomainFromField_management(fc, dataRefField)
            except Exception as e:
                print e.message

            try:
                # Assign domain to field
                arcpy.AssignDomainToField_management(fc, dataRefField, domainName)

            except Exception as e:
                print e.message

This function works fine for 100 loops but exits after that without any error or warnings. I tried to iterate loop by just printing output which give me ~500 rows. Could anyone please help me to identify the problem here?

Edit: Just tried creating more than 100 domains and assigning them to field which worked well so there must be something wrong in this code


texts: Using QGIS Slider for filtering based on attributes?

texts:
For each of my layers, I have various attributes that I might be filtering on, like attribute A < 1000 AND B > 300 etc..

Just like we can add an opacity slider to the legend, I want to be able to add sliders for the various attributes I am filtering on, without having to go into the filter screen and manually change the filter text.

I have seen this answer:
Slider for non time attributes in QGIS

And while the hacky solution does work (use Time Manager), it doesn't really work if you want to have multiple sliders to explore your data.

I'm thinking of making something along these lines, but before I do that I figured it's such a useful feature that perhaps there already is a plugin that does it, or some other mechanism. I couldn't find it.


texts: Why does loading speed of GeoServer services vary so much, how to optimized settings?

texts: We are running web services in GeoServer 2.14. Loading tiled WMS layers in an openlayers application. Polygons, lines and points. One layer published from a table holding several million records and growing. Backend is PostgreSQL/PostGIS database. Spatial indexes added to all data. Drawing scale limits imposed. The services will load fairly quickly at some points and very slowly at other times. We have not really modified the settings in terms of the connections pool or anything else. The slowness is not my local network, it is definitely specific to the application/web services as I have no issue with other sites. I'm unsure of how to troubleshoot but think that there must be some optimization of the configuration that could take place. I will add that aside from the drawing time, when you click on the feature to access the popup, at times it responds right away, at other times, when the features are slow to draw, it can take awhile for the popup to appear after clicking on the feature.


texts: How to merge multiple GeoPackage files into one file with a single layer using ogr2ogr?

texts:
I'm using ogr2ogr in a bash script to merge multiple gpkg files. What I'm doing basically looks like this:

ogr2ogr -f 'gpkg' -append merged.gpkg input_file_1.gpkg
ogr2ogr -f 'gpkg' -append merged.gpkg input_file_2.gpkg
ogr2ogr -f 'gpkg' -append merged.gpkg input_file_3.gpkg

So all the input_file_x.gpkg end up in merged.gpkg. This works fine, except a new vector layer is created for every input file, which I want to avoid. As far as I understand it, the -append flag should take care of this.

From man ogr2ogr:

  -append: Append to existing layer instead of creating new

All of the input files look identical in terms of attributes and geometry type.


texts:
I recently discovered the Parallel Processing environment setting in ArcGIS (this must be new to 10.1), however, there is no mention as to what tools are designed to utilize this environment setting.  ESRI provides one example of how to set the parallel processing environment:

import arcpy

# Use half of the cores on the machine.
arcpy.env.parallelProcessingFactor = "50%"

My questions are very closely related, so I am including them in the same thread:

Which geoprocessing tools honor the parallel processing environment?
Are these local or global settings (i.e. can you set the environment
at the beginning of the Arcpy script and all respective tools will
honor the environment setting thereafter?
Are most geoprocessing tools already set to:
arcpy.env.parallelProcessingFactor = "100%"by default?


texts: Yandex Navigator / Maps - Determining user location

texts:
Kindly redirect me to a more appropriate forum in case there exists one.

I am an active user of Yandex.Navigator App (available for all kind of mobile platforms). Whenever I open the App, even though my GPS is turned off, it shows me my current location upto a great extent of precision on the map. How they are obtaining my location without asking me to turn GPS on? Are they using my mobile network connection, 3G internet connection, some new technology or what?


texts: vba code for algorithm to assign track file gps UTM coordinates (East Africa) to an overlay 5x5km grid cell

texts:
The GPS coordinates are temporarily imported in to an Access database from DNR Garmin files (downloaded as dbf file(s)), and then exported for mapping in AEJEE (ArcExplorer Education Edition).

I need to be able to run the vba code on importing the dbf aerial track files and populate another database table with the grid cell reference associated with the track file coordinates, so that I get grid cell 'hits' per flight -          ie. the number of times each 5x5km grid cell has been 'visited' by the aircraft.

Any help with this would be fantastic!


texts: Comparing ESRI Demographic Data from different years

texts: I have in my possession ESRI Demographic update data for multiple years: 2013/2018, 2014/2019, 2015/2020, 2016/2021, 2017/2022.  I am currently trying to see if it is possible to compare each of these data sets against each other to ascertain potential customers for highspeed internet service in my county.  For example, a time series analysis looking at population change using the total population estimates for each year.  I have tried looking up this information on ESRI's website but I have not found anything that could help me answer my question, so I am hoping someone on the stack exchange with knowledge of ESRI Demographic Data could provide with one.  If it is possible to do this, I am also considering using other variables too such as median household income, educational attainment, disposable income, age and other socioeconomic measures for marketing purposes.


texts: Initialise a database with SpatiaLite 4.3

texts:
I am starting to work with SpatiaLite 4.3 in an empty database. SpatiaLite requires a number of meta-data tables, where it stores information about geometry fields. For previous versions there are SQL scripts that create and initialise these tables; for instance, the release page for version 2.3 lists a script named init_spatialite-2.3.sql.

I have not been able to find a similar script for version 4.3. Is the database initialised in different in this version?


texts: How to emulate the modulo operator in a style rule?

texts:
I have SHP layer loaded in QGIS. It consists of lines with attribute column z, containing integer value associated with line elevation in meters.
Lines are drawn at 10 meter increment, and I want to apply style to lines on 50 meter increment.

For that purpose in "Style" tab in "Layer Properties" dialog, I use "rule based" style and want to define above described rule.
From available operators I can see only one possibility with this filter expression:

z/50 - toint(z/50) = 0

However it returns true for all lines, while I expected just those with 50 meter increment.

There is also "Expression builder" dialog, but it does not return results from typed expression, for some reason, and I'm not able to determine what is the problem?

Any ideas?


texts: Converted GeoJSON have 3 elements in each coordinate entry


texts: Why OpenScales show tiles in cache in zoom in/out WMS layer?

texts:
I use OpenScales to show a WMS layer and it is strange that when zoom in/out sometimes OpenScales show tiles in cache, so Multi-resolution tiles showed in the layer.
When paning the layer, incorrect tiles will disappear. How can I fix this? Thanks!!!


texts: How to parse an ArcGIS json in google code?

texts:
I'd like to be able to programmatically query a feature service api and pull the attribute data (I don't really need the geometry for this) into a google sheet. It seems like it ought to be simple but maybe I'm just missing something. I've been playing around with the code at this link here gist and below but so far haven't been able to get it to work, and it seems like it might be out of date.

I'm currently stuck at the code below. When I run the debugger, it looks like while dataAll is populating correctly, the feature variable is not. It looks like 'dataAll[i]' isn't working. Does anyone know why or how to fix it?

  function pullJSON() {   var ss =
  SpreadsheetApp.getActiveSpreadsheet();   var sheets = ss.getSheets();
  var sheet = ss.getActiveSheet();

var url="http://your url here"; // Paste your JSON URL here
 var response = UrlFetchApp.fetch(url); // get feed

   var dataAll = JSON.parse(response); //   Logger.log(dataAll)
for (i = 0; i<30; i++) {
  feature = dataAll[i];
  name = feature.NAME;
  Logger.log (name)
}


texts: Generate node / vertex at starting point of line connections

texts:
I have line and point layer. I connected line layer with points using connect nodes to line in processing toolbox. After creating a connection line (green color), I used line intersection from processing toolbox to extract the nodes/vertex of each starting point of the connection line as shown in red color. The problem here is:

The line intersection generates nodes for every intersection of line as circled with black color in image. I only want to extract nodes of lines which connect point and a line layer.
I also tried Extract specific vertices tool as suggested in previous questions (here), but it does the same thing as line intersection function.


texts: Changing index order of points within the same polygon

texts:
I am currently trying to get rid of a waterway that divides my polygon. The problem is that the points are part of the same polygon but are around 190 rows away for each other.

Is there a way to rearrange the point so that the points in the black circle are connected to each, same for the red circled points?


texts: Google Earth Engine: "Parameter 'index' is required" error when mapping to an Image Collection

texts:
I am trying to apply a distance function to a collection of masked images, but am receiving a "Parameter 'index' is required" error. I've found answers for other "Parameter 'xyz' is required/missing" errors but nothing for 'index'.
The function works when I input independent years (e.g. '7' or '12') but not when applied to over the list of years. It seems like there is some trouble indexing each image after applying the function but I'm not sure why. I've tried setting the system index (set("system:index", y)) within the function but still no luck. I've even tried nesting functions instead of having two arguments within one function, but the same error shows up. I'm by no means a GGE expert so please let me know if I'm missing something easy!
Here is the code for the non-nested function modified to an 'ImageCollection' for generalizability:
// List of last two digits of years 2000-2016 to call from a list of 17 tree cover images (1 from each
// year)

var years = ee.List.sequence(0,16);

// function to calculate distance (for each image in an Image Collection) from each masked pixel to
// nearest unmasked pixel

function distance(coll, y) {
  var list = coll.toList(years.size()); // converting to list to call by indexed number
  var non_tc = ee.Image(list.get(Math.round(y))).lte(0);   // using 'round' to convert floats to integers
  var tc = ee.Image(list.get(Math.round(y))).gt(0);
  var dist = non_tc.cumulativeCost({
                      source: tc,
                      maxDistance: 10000
    });
  return dist
    .set('year', y);
}

// When applied to one year, it works fine:
var test12 = distance(ImageCollection, 12);

// When I try to apply it to all years, I get the error:
var dist_coll = distance(ImageCollection, years);
print('dist_coll', dist_coll);

Here's another version of the function using "map", which returns the same error:
function distance1(coll) {
  var list = coll.toList(years.size());
  var non_tc = ee.Image(list.get(Math.round(yeears))).lte(0);
  var tc = ee.Image(list.get(Math.round(years))).gt(0);
  var dist = non_tc.cumulativeCost({
                      source: tc,
                      maxDistance: 10000
    });
  return dist
    .set('year', y);
}

var dist_coll1 = years.map(distance1(ImageCollection));
print('dist_coll1', dist_coll1);

And trying to nest functions, with the same result:
var distance2 = function(coll) {
  var dist = function(y) {
       var list = coll.toList(years.size());
       var non_tc = ee.Image(list.get(Math.round(y))).lte(0);
       var tc = ee.Image(list.get(Math.round(y))).gt(0);
       var cost = non_tc.cumulativeCost({
                        source: tc,
                        maxDistance: 10000});
          return cost
            .set('year', y);
      };
   return coll.map(dist);
};

var dist_coll2 = years.map(distance2(ImageCollection));
print('dist coll2', dist_coll2);


texts: Comparison between using VB.net or C# in ArcObjects or ArcGIS Server

texts: Is there any differences between using VB.net or C# in ArcObjects or ArcGIS Server based projects, in functionalities, samples, Q&A community? Is there any comparison which covers these aspects to help choose between them ?


texts: R: jsonlite/geojsonio: Why are these geoJSON polygons being read in different formats? And how to deal with them for shp exportation?

texts:
EDIT: Turns out fromJSON() has a variable called simplifyVector that by default goes TRUE. By setting it to FALSE I was able to turn the string into Json, from Json into geoJson and from geoJson into spatialpolygonsdataframe. Still don't know why it happens, though.

I have these two polygons in geoJSON format I got from an API and I want to export them to shp. When I read these two polygons using fromJSON() I get two different formats. Apparently, the reason would be because the second polygon includes a third coordinate in some of its vertex. Why is it happening?

It's worth noting that I've been successfully able to export the polygons to shp in the second format of the reproducible example below, but I haven't had the same success with the first format. How could I "standardize" them so I could  read multiple polygons in these two possible formats, and export them to shp files?

Here's a reproducible example:

Polygon1 <- '{"type":"Polygon","coordinates":[[[-70.3248,-27.3708],[-70.3247,-27.3713],[-70.3246,-27.3719],[-70.3245,-27.3723],[-70.3239,-27.3743],[-70.3198,-27.374],[-70.3194,-27.3724],[-70.3193,-27.3719],[-70.3196,-27.3704],[-70.3201,-27.3696],[-70.3204,-27.3695],[-70.3222,-27.3693],[-70.3248,-27.3708]]]}'

Polygon2 <- '{"type":"Polygon","coordinates":[[[-70.597,-33.4169,0.0],[-70.5952,-33.4144],[-70.5947,-33.4139],[-70.5926,-33.4133],[-70.5921,-33.4133],[-70.5907,-33.4138],[-70.59,-33.4142],[-70.5885,-33.4152],[-70.5894,-33.4163],[-70.5916,-33.4195],[-70.5918,-33.4197],[-70.5929,-33.4191],[-70.5943,-33.4184],[-70.5955,-33.4177],[-70.597,-33.4169,0.0]]]}'

Geojson1 <- fromJSON(Polygon1)
print(Geojson1)
Geojson2 <- fromJSON(Polygon2)
print(Geojson2)

And the output of the above code:

print(Geojson1)
$type
[1] "Polygon"
$coordinates
, , 1
         [,1]     [,2]     [,3]     [,4]     [,5]     [,6]     [,7]     [,8]     [,9]    [,10]
[1,] -70.3248 -70.3247 -70.3246 -70.3245 -70.3239 -70.3198 -70.3194 -70.3193 -70.3196 -70.3201
        [,11]    [,12]    [,13]
[1,] -70.3204 -70.3222 -70.3248
, , 2
         [,1]     [,2]     [,3]     [,4]     [,5]    [,6]     [,7]     [,8]     [,9]    [,10]
[1,] -27.3708 -27.3713 -27.3719 -27.3723 -27.3743 -27.374 -27.3724 -27.3719 -27.3704 -27.3696
        [,11]    [,12]    [,13]
[1,] -27.3695 -27.3693 -27.3708

print(Geojson2)
$type
[1] "Polygon"
$coordinates
$coordinates[[1]]
$coordinates[[1]][[1]]
[1] -70.5970 -33.4169   0.0000
$coordinates[[1]][[2]]
[1] -70.5952 -33.4144
$coordinates[[1]][[3]]
[1] -70.5947 -33.4139
$coordinates[[1]][[4]]
[1] -70.5926 -33.4133
$coordinates[[1]][[5]]
[1] -70.5921 -33.4133
$coordinates[[1]][[6]]
[1] -70.5907 -33.4138
$coordinates[[1]][[7]]
[1] -70.5900 -33.4142
$coordinates[[1]][[8]]
[1] -70.5885 -33.4152
$coordinates[[1]][[9]]
[1] -70.5894 -33.4163
$coordinates[[1]][[10]]
[1] -70.5916 -33.4195
$coordinates[[1]][[11]]
[1] -70.5918 -33.4197
$coordinates[[1]][[12]]
[1] -70.5929 -33.4191
$coordinates[[1]][[13]]
[1] -70.5943 -33.4184
$coordinates[[1]][[14]]
[1] -70.5955 -33.4177
$coordinates[[1]][[15]]
[1] -70.5970 -33.4169   0.0000


texts:
Spatial graph databases are capable to find a shortest path from the node A to the node B. Are there any DB which allows complex conditions on nodes and edges between them?

For example:

The node is closed between 01:00 - 06:00 UTC
The edge is available only on weekends
The edge is available if we passed the node X and going to the node Y
The edge is not available in winter
The node is temporarily disabled
The polygon closes all edges intersecting it
etc

Are there any DBs which allows such geo-restrictions?

Say we can call some function when considering the edge as a candidate?


texts:
To say it upfront: This is a question that has already come up before in Stack Overflow but did not provide a manageable and legal solution for my problem (see below).

I have a database (in STATA format) of around 1,3 million addresses in Germany with street name, house number, zip code and city stored separately in different variables.
For each of those I need the longitude and latitude which I then want to store in a (local) database and use later on in a statistical analysis (in the context of a research project). Note that this is a one-time request for each of those addresses!

Essentially, this poses four requirements on the service to be used:

It needs to cover Germany (in particular its rural areas)
It needs to be sufficiently exact, i.e. I am willing to take the imprecision of around 100m relative to the object placed at the corresponding address, preferably less of course
My usage needs to comply with the terms of services
It should be for free or the price should be reasonable / It should take a reasonable time to obtain the data (<= 4 weeks)

Before I look into licenses for commercial usage, I want to get an overview of free / reasonable priced services. Seen against this background, the third condition seems to rule out the big players like Google and Yahoo since they do not allow to store the data (see also discussion in GIS SE at Geocode address once in google maps, save in database, read geocoded result from there on).

I checked with an employee of Google Germany who claimed they would charge 10,000 EUR (which appears to be way too high, though).

Do any services come to your mind that suit those requirements?
If so, do program codes exists to collect and store the data?

I have also raised this question in Stack Overflow where it was put on hold. Since I did not find restrictions for this kind of "preference question" for the GIS Q&A section, I am posting it here.


texts: Is it possible to build up a system based only on individual markers without postGIS?

texts:
We are a team of 2 developers and we've been thinking this month on how to build up a system for showing locations on a map. After reading some articles and taking some courses seems like the best aproach would be PostGIS for saving the data.

First some background. We have 2 months to build a demo and the stack looks as follows:

ReactJS
Leaflet
PostgreSQL
Node

Since we only need to show certain locations on maps, seems like geometries are not neccesary. So the way that we would save the location would be just LatLng and the object itself (user, place, etc...) in standard tables. Distances would be calculated by leaflet if neccesary. Tiles would be provided by CartoDB.We would build a simple function to translate the data to geoJSON before passing it to leaflet.

My question is:

Should I consider PostGIS or just standard SQL table structure?

We seek any suggestion or corrections since this is our first app with gis data.


texts: Display the gml:id of WFS layer as an attribute in QGIS

texts:
I have loaded a WFS layer in QGIS and want to display the gml:id of the features in the attribute table. Is this somehow possible?
I can see the gml:ids in other applications (ex FME).


texts: Moving Table of Contents, Search window, Toolbox, etc to second screen?

texts:
Table of Contents, Search, Toolboxes etc. all compete for horizontal space on the computer screen, often obscuring the map display itself.

If I use two screens can I put these windows on one screen and keep my map display on the other?


texts: Is there a way to reference a selection in model builder or python?

texts:
Is there a way to reference a selection in model builder or python so that I can edit my python script to only run if there is a selection. And even better if there is only one row selected?

I made a model to make a selection with python script to ask for variables, make a selection based on those variables, and then zoom to the selection. When the tool does not make a selection for whatever reason it takes forever to finish and zooms to the entire extent of the map or crashes ARCmap entirely.


texts:
I'm trying to sort a dataset by field StartTime (datatype date, indicates when a technician is supposed to start a visit). I am including a SQL where_clause to avoid an earlier problem when it complained it couldn't compare a date to a NULL value.

for Technician in TechList:
    insertFields = ["StartTime", "AssignedTech", "VisitOrder"]
    sql = """"StartTime" IS NOT NULL"""
    print sql
    counter = 1
    with sorted(arcpy.da.UpdateCursor(fcStops, insertFields, sql)) as cursor:
        for row in cursor:
            if row[1] == Technician:
                row[2] = counter
                cursor.updateRow(row)
                counter += 1

The problem is in the sorted (an idea I got from the ArcMap SearchCursor help page, Example 5). If I take that out, it runs but of course is all out of order.

How can I sort the UpdateCursor so that I'm assigning VisitOrder the correct ordered value -- or is there another way ensure that the first visit gets 1, second 2, etc.?


texts:
I want to set the rotation of an ol.style.Icon to a given feature property. See line 12 of this:

var styleCache_airports={}
var style_airports = function(feature, resolution){
    var value = ""
    var size = 0;
    var style = [ new ol.style.Style({
        image: new ol.style.Icon({
                  imgSize: [32, 32],
                  scale: 0.59375,
                  anchor: [9, 9],
                  anchorXUnits: "pixels",
                  anchorYUnits: "pixels",
                  rotation: function(feature){return feature.get("ELEV") * 0.0174533;},
                  src: "styles/a18_sluis.svg"
            })
    })];
    if ("" !== null && resolution > 28 && resolution < 1400) {
        var labelText = String("");
    } else {
        var labelText = ""
    }
    var key = value + "_" + labelText

    if (!styleCache_airports[key]){
        var text = new ol.style.Text({
              font: '10.725px \'MS Shell Dlg 2\', sans-serif',
              text: labelText,
              textBaseline: "center",
              textAlign: "left",
              offsetX: 5,
              offsetY: 3,
              fill: new ol.style.Fill({
                color: "rgba(0, 0, 0, 255)"
              }),
            });
        styleCache_airports[key] = new ol.style.Style({"text": text})
    }
    var allStyles = [styleCache_airports[key]];
    allStyles.push.apply(allStyles, style);
    return allStyles;
};

This doesn't work, and doesn't throw any errors. If I replace it with a static value, it works. What am I doing wrong?


texts: Require a 3rd Party GIS service to embed in my application

texts:
I am looking for a service (Free or Premium) that allows offline access to maps and also gives the functionality to search map using an address or Latitude / Longitude. I want to embed it into a Web Application. I am particularly interested in Maps. Please tell me if any one know of any such service. Thanks.

Note: Offline means the map's data will not reside on the Internet but on a private (intranet) server.


texts: ST_GEOMETRY: spatial column must appear immediately after spatial relationship function for spatial index to be used

texts:
The docs for SDE.ST_GEOMETRY have a blurb that says the following:

When are spatial indexes used with ST_Geometry?
...the spatial column must appear immediately after the
spatial relationship function in the query for the spatial index to be
used.

I'm not clear on what that means or why it's necessary.

That blurb seems to imply that spatial relationship functions would only be used in the SELECT clause. I would have thought that spatial relationship functions could be used in the FROM and WHERE clauses too.

Question:
What does the documentation mean by "must appear immediately after the spatial relationship function" and why is it necessary?

Oracle 18c; ArcGIS Enterprise 10.7.1 (the docs seem to apply to PostgreSQL too)


texts: GCS_WGS_1984 and 900913. What is the difference?

texts:
I use WMS Service with mymapproxy/mapproxy.yaml :

services:
  demo:
  kml:
  tms:
  wmts:
  wms:
    srs: ['EPSG:4326', 'EPSG:900913']
    md:
      title: MapProxy OSM WMS Proxy
      abstract: This is the fantastic MapProxy.
layers:
  - name: osm
    title: Local OSM WMS
    sources: [osm_cache]

caches:
  osm_cache:
    grids: [global_mercator_osm, global_geodetic_sqrt2]
    sources: [osm_mapnik2]

sources:
  osm_mapnik:
    type: mapnik
    mapfile: /home/oleg/src/mapnik-style/osm.xml
  osm_mapnik2:
    type: tile
    grid: global_mercator_osm
    url: http://localhost/osm_tiles/%(z)s/%(x)s/%(y)s.png

grids:
  global_mercator_osm:
    base: GLOBAL_MERCATOR
    srs: 'EPSG:900913'
    origin: nw
  global_geodetic_sqrt2:
    base: GLOBAL_GEODETIC
    res_factor: 'sqrt2'
    srs: 'EPSG:900913'
    origin: nw

Then I display this WMS service in Geoserver. I have not any problems (if i set Native и Declared SRS as 900913)

But I'd like display this service in ArcMap 10.1

This is the coordinate system of my ArcGis service.

So, I don't see the 900913 coordinate system. And when I add this service to TOC (left panel with layers in ArcMap) I have a lot of errrors:

Sometimes I do not see a map or this map is not good:


texts: Getting ExecuteError: ERROR 001558: Error parsing JSON file in ArcPy?

texts:
I'm trying to convert some open-source data in JSON format to a shapefile for use in ArcGIS Desktop.  I'm attempting to do this in ArcPy with ArcMap 10.7.  The data I'm trying to access is here, in JSON format:

data from data.SF.gov

And my script:

#access the JSON from web link
bikeLanesRaw = 'https://data.sfgov.org/resource/ygmz-vaxd.json'
bikeLanesRawResponse = urllib2.urlopen(bikeLanesRaw)
bikeLanesRaw_JSON = bikeLanesRawResponse.read()
#Write each response to JSON
foo = open(bikeFacilitiesFolder + "\\" + "bikeLanes_SF.json", "wb")
foo.write(bikeLanesRaw_JSON);
foo.close()

#Convert JSON to shapefile - THIS FAILS
arcpy.JSONToFeatures_conversion(bikeFacilitiesFolder + "\\" + "bikeLanes_SF.json", bikeFacilitiesFolder + "\\" + "bikeLanes_SF.shp")

Everything works as expected until the final line, where I try and convert the JSON to features.  It gives me the error:

 ExecuteError: ERROR 001558: Error parsing json file 'Bike Facilities Data\bikeLanes_SF.json'.
Failed to execute (JSONToFeatures).

I've read more of ESRI documentation, and they say:

Review the JSON structure in the input file. The JSON must have at least the geometryType, spatialReference, fields, and features (with geometry and attributes) property.

What do I do if the JSON does not have this information to begin with?


texts: How to batch process spatial analyst tools?

texts:
I have written the following code to convert my point files of daily max temps to rasters:

# Set local variables for Point to Raster Conversion
inFeatures = Day12015Max
valField = "Max_2015"
outRaster = ("Rstr12015Max")
assignmentType = "MOST_FREQUENT"
priorityField = ""
cellSize = 0.05

# Point to Raster conversion for Max values
arcpy.PointToRaster_conversion(inFeatures, valField, outRaster, assignmentType, priorityField, cellSize)

I have years worth of data so this will be 365x files to process. If all my point files are in a folder, how can I write a script that will go into the folder and go down the list and convert everything and output to a new folder?

I am working with Python 2.7 and ArcGIS 10.3

I have tried the following method:

env.workspace = "C:\Users\uma.bhandaram\Desktop\DailyMax2015\Shapefiles"

#folder to list
shpFileList = arcpy.ListFiles("*.shp")

for shpFile in shpFileList:
    # geoprocessing steps
    shpFileName = os.path.splitext (shpFile)[0]
    rastFile = shpFileName + "Rast.tif"

arcpy.PointToRaster_conversion(shpFile, "Max_2015", rastFile, "MOST_FREQUENT", "", 0.05)

But, this only processes one file and does not seem to loop through.


texts: GeoServer style multiple label inside shield

texts:
I want to add multiple labels inside a circle shield.

I also want to set different locations for all labels inside the circle shield.

I made this for one label.

How can I do this for multiple labels?

If multiple labels are not possible, how can I change the label location inside the shield?

* {
    stroke: black,lightgray;
    stroke-width: 3,2;
    label: [name];

    font-family: 'Ariel';
    font-size: 5m;
    font-fill: black;
    shield: symbol(circle);
}
:shield {
    fill: white;
    stroke: black;
    size: 10m;
}


texts: Can a shapefile, or other file, be made from Excel without coordinates

texts:
I have two Excel spreadsheets I'd like to make into shapefiles, or some equivalent. I have access to ArcGIS at work and want to try using GRASS as well. I have seen the help on doing so, but only if I have Lat and Long.

What I'm trying to do is create a shapefile from an old, 1870s, map that has property owner names on them. I've warped the map and I have an Excel file that has all the names already in it. I'd like to just be able to select a name and put a point where it would be on the map.

I also have a list of cemeteries, some with coordinates, that I'd like to map out. But since not all cemeteries have been located I don't have good GPS.

My question is, can I do this? Can I just make empty Lat/Lon columns? Should I put a 0 value in each cell? Or is this not going to work at all?

edit:
The maps are in no way georeferenced, they are not that well done either, at least cartographically. They are useful if you know the road name and one can get an idea of where the property used to be. There are some 6000+ properties shown, I basically wanted to remake the map to make it easier to find the properties.

https://jscholarship.library.jhu.edu/bitstream/handle/1774.2/34119/p.07%20Buckeystown%2c%20District%201.jpg?sequence=44

is what one of the maps look like. They are all for one county in Maryland. I already have the county outline, roads and other features. I also have a spreadsheet with all the names typed out, I just didn't want to have to retype 6000 names.

2nd edit:
I guess I'm confused as to why it seems so difficult, what I was hoping to be able to do is just select a name and say, place a point here. I understand that the program wants a GP, but I would have thought it was possible to add them later.

As for the maps, they do not warp well at all. They are not the type of map that was surveyed at all. When I went to warp them it didn't work out at all. It works locally, but on the whole it doesn't. The exact placement isn't the point, just a general idea of the location. So I planned on warping as best as possible and then shifting the raster map around. The roads from then pretty much match the roads now so getting a general location is easy enough. It's just I'd rather not retyp3 6000+ names if at all possible.


texts: ArcGIS Calling Wrong fme.exe When Python Script Converted To An ArcGIS Script Tool

texts:
I'm using this this bit of code to execute a FME Workbench file through a Python script ...

import os
os.system("fme.exe P:\\Mapping\\Scripts\\FME\\Easements_ROW_Extract_From_CAD.fmw")

However, when I try using it as an ArcGIS script tool, ArcGIS tries to execute it using the Data Interop's fme.exe. I only have the stand alone version of FME.  I tried using the full path to the fme.exe i'm trying to use, but that's not working ( the script runs without error, but nothing is produced).  Any ideas on how to get this script to work, using the stand alone version of FME,  as an ArcGIS script tool?


texts: What's the fastest way to convert large csvs to point shapefile

texts: I'm working in QGIS, and my current method to change csvs to shapefiles is to import them to QGIS, right click and save as shapefile. For large CSVs (over 100 million rows) it can take quite some time to save. Is there a quicker, more efficient way to do this? Perhaps using geopandas?


texts: Is there any tool to execute Postgis-like queries without a postgres server?

texts:
I would like to know if there is any way/tool of using geoprocesing queries like the ones you use with POSTGIS but without the need of having an actual database server available.

The particular aplication in this case is the linear referencing of arbitrary points over a set of linestrings representing a road network. Basically I'm using a "static" postgis database (static in the sense that there is just an initial load and no further updates/deletes/inserts).

I've heard of Spatialite but I'm not sure what kind of infrastructure do you need to do this kind of operations or if it has the same capabilities than the POSTGIS extension.

The actual aplication is being written in c#, although more than native libraries, I would be interest in running queries with little modifications, as I already have previous knowledge on POSTGIS queries.


texts:
I need a solar radiation data in central Asia. Is there any global dataset of solar radiation?

If someone know about it, please send to me the link.


texts: Autoregressive moving average (ARMA) model in GEE

texts:
Is there a way to implement an autoregressive moving average (ARMA) model in GEE? For example, for MODIS NDVI dataset spanning 1 year:
var modisNDVIcoll = ee.ImageCollection("MODIS/006/MOD13A2");
var modisNDVIcoll = modisNDVIcoll.filter(ee.Filter.calendarRange(2009, null, 'year'));


texts: Which framework is best to implement an application robustly, using ArcGIS Server JS API?

texts:
Over on the Esri forum, user Matt Moyles suggested that the approach used in the Esri JS samples isn't suitable for robust development of a web mapping application using HTML, JavaScript and CSS:

ESRI's suggested approach to javascript application development is ancient and outdated. I would not recommend following the samples. Dojo 1.7 supports AMD with Asynchronous dependency loading. I would start with the dojo boilerplate template and "work" the arcgis api into that. The samples are not suitable for serious application developers. They are mostly just proof of concept snippets.
A serious application should be developed with the aid of some sort of framework to help structure things. I have been using dojox.mvc with great success! But other options include things like backbone.js, spine.js, or even javascript MVC.
Dojo Boilerplate - https://github.com/csnover/dojo-boilerplate

does anyone agree/disagree with this statement?
are there any online examples of better approaches to web application development using the ArcGIS Server JS API?


texts: Draw line to next feature with distance as label (dynamic)

texts:
I would like to dynamically draw the shortest line to the next object (Geometry Generator?).

Example: I have a network-layer (streets) and a point-layer. I know the distance from every point to the closes street (static, by using v.distance or one of the plugins NNjoin or refFunctions) but I would like to have a dynamically drawn line to the next street with the distance as label. It's important to have it dynamic, so if I move a point, the line and distance should update.
Since I use GeoPackage a SQL-Solution would also be nice.

In the screenshot below you can see my static solution:


texts: Computing distance of a georeferenced .jpg to each polygon (town) in a shapefile

texts:
I have a .jpg that is a highway that spans multiple towns, and a shapefile with boundaries of these towns.  I have georeferenced the .jpg so that it is aligned with the shapefile of the towns.  My ultimate goal is to compute the distance of each town--ideally from its centroid--to the nearest point along this highway.

I am able to convert the georeferenced highway into a shapefile.  I have converted the shpaefile of the towns to a set of points to ease the distance calculation.  When I do the Near calculation under Toolbox, however, I do not get any distances computed.

When I try to do a Near calculation, I get Null for the distance.  When I try a Spatial Join, I get a -1.  I cannot figure out what issue I can be having.

I am using Arcmap in ArcGis 10.2.


texts: Missing libraries when upgrading to PostGIS 2.1 and PostgreSQL 9.3.1 using homebrew

texts:
In the process of upgrading my PostgreSQL from version 9.2.4 to 9.3.1 (via homebrew on OS X) I came across an odd problem. These are the steps I took so far

PostgreSQL, PostGIS and required libraries installed (no errors)
run initdb on the new database
stopped both servers
running pg_upgrade

pg_upgrade performs the necessary checks, creates dumps of the old cluster, but when importing into the new cluster I get the following error:

> ./pg_upgrade -b /usr/local/Cellar/postgresql/9.2.4/bin/ -B /usr/local/Cellar/postgresql/9.3.1/bin -d /usr/local/var/postgres/ -D /usr/local/var/postgres9.3.1 -u postgres
Performing Consistency Checks
-----------------------------
Checking cluster versions                                   ok
Checking database user is a superuser                       ok
Checking for prepared transactions                          ok
Checking for reg* system OID user data types                ok
Checking for contrib/isn with bigint-passing mismatch       ok
Creating dump of global objects                             ok
Creating dump of database schemas
                                                            ok
Checking for presence of required libraries                 fatal

Your installation references loadable libraries that are missing from the
new installation.  You can add these libraries to the new installation,
or remove the functions using them from the old installation.  A list of
problem libraries is in the file:
    loadable_libraries.txt

Failure, exiting

It appears as though PostgreSQL 9.3.1 tries to use PostGIS 2.0, which is not compatible

Could not load library "$libdir/postgis-2.0"
ERROR:  could not access file "$libdir/postgis-2.0": No such file or directory

Could not load library "$libdir/rtpostgis-2.0"
ERROR:  could not access file "$libdir/rtpostgis-2.0": No such file or directory

Has anyone run into the same problem?


texts: Can an ArcGIS Desktop Add-In Tool have two or more input events?

texts:
I'm trying to create an add-in for ArcMap in Visual Studio that adds points on mouse clicks and then executes some geoprocessing on those points when the enter key is pressed.

So far sifting through the ArcObjects help has yielded a few things that I've tried unsuccessfully to put together:

http://resources.arcgis.com/en/help/arcobjects-net/componenthelp/index.html#//001v000000pq000000
http://resources.arcgis.com/en/help/arcobjects-net/componenthelp/index.html#/OnMouseDown_Method/001v000000w1000000/

I'm wondering if it's even possible to have more than one of these "events" activated by a tool or if I'm barking up the wrong tree and should be attempting a different method?

Here is a skeleton of some code I have so far:

namespace MultiEventTest
{
    public class Tool1 : ESRI.ArcGIS.Desktop.AddIns.Tool
    {
        public Tool1()
        {
        }

        protected override void OnUpdate()
        {
            Enabled = ArcMap.Application != null;
        }

        protected override void OnMouseDown(MouseEventArgs arg)
        {
        //code to add points to feature class
        }

        protected override void OnKeyDown(KeyEventArgs arg)
        {
            if (arg.ModifierKeys == (System.Windows.Forms.Keys.Enter))
            {
                //code to execute geoprocessing
            }
        }

        }
    }


texts: How to calculate the extent of a projection?

texts:
I am trying to change the projection of my Open Layers map from 'EPSG:3857' to 'EPSG:26915'. I tried to use

var projection=new ol.proj.Projection({
    code: 'EPSG: 26915',
    extent:extent_projection,
    units: 'm'
  });

  var  view =new ol.View({
         projection: projection,
         center: [-57685.58, 4615858.09],
         zoom: 2
        });

The map looks extremely blurry. I believe that I need correct extent_projection to remove the blurriness of my openLayers map. I tried going to www.epsg.io and used projected bounds of EPSG:26915 for my extent array. However, it didn't change anything. The following is the image of my blurry map. I need to correct the center and extent of my map. How can I get the extent for my projection?


texts: Why doesn't wx graphical modeler Grass gis 7 work properly in Spanish?

texts:
I am working with the wx graphical modeler in grass gis 7.

When I try to use a variable call "archivo" in spanish, wich is "file" in english, it does not work. The variable dont let me choose a file.

But if I change the language to english it work fine. even if having the interface in spanish I use the name type "file" in place of "archivo" it works.

Anyone know why this is happening?

This are some screenshot
This is in spanish

This in English

This in Spanglish


texts: Remote sensing data for different scales

texts:
I'm looking for different satellite data concerning different scales. The first one I'm looking for is remote data for soil sealing and water bodies at a scale of 1:30,000, could be generally everywhere. The other one I'm looking for is for data with Mapping of parks in urban areas and the sub-classes: lawn, bushes, trees, tree groups, water bodies, lane (sealed) and lane (unsealed) at a scale of 1:5000.
I was checking the page of USGS but it doesnt really provide a lot of satelitte data for the scales I am looking for. Could anybody recommend some descent pages where I can get this images without have to pay too high amounts to get the info?


texts: How display outlines in a Map, only from a scale, but the colored polygon in all scales?

texts: I want to display a shape with UMN-Mapserver. It has different coloured polygons and also outlines. But I want to display the outlines only, from a specific scale. I could add MINSCALE 150000, but in that case the polygon is completely not rendered, if that doesn't match. How I display the polygon without outline in that case?


texts: Find mean of a group of every 10 layers in raster stack in R

texts: I have a raster stack (s1) of 2000 layers and I need to get mew stack of 20 layers from this (s2). Such that first layer in s2 is mean of first 10 layers of s1 and second s2 layer is mean of 11-20 layers of s1 and so on. I need to do this in R. I can use calc but I am confused what function would be like?


texts: delete the markers in routing mode of MapQuest

texts:
I am using MapQuest plugin for LeafletJS.

After it draw the map, it sets a marker on each stop point on the route, e.g: A, B, C, D, etc.

Is there a way to remove those markers, I just need blue routing lines.

In Google Maps, there is an option, such as:

suppressMarkers: true

Is there a way in MapQuest?


texts: Connection Error on a PostgreSQL x64 9.5 database from ArcGIS Pro 1.3.0

texts:
I have the following PostgreSQL database named "cfa_fem_gc":

CREATE DATABASE cfa_fem_gc
  WITH OWNER = gc
       ENCODING = 'UTF8'
       TABLESPACE = pg_default
       LC_COLLATE = 'English_Australia.1252'
       LC_CTYPE = 'English_Australia.1252'
       CONNECTION LIMIT = -1;

ALTER DATABASE cfa_fem_gc
  SET search_path = "$user", public, tiger, sde;

It was enabled with PostGIS 2.3 and all GIS data has been stored on the "gc" schema with the owner "gc".

I do not have any difficulty connecting to the database and viewing the shapefile on QGIS.

However, a connection error was thrown from ArcGIS Pro 1.3 if I tried to make a new Database Connection and connect to the cfa_fem_gc database.


texts: Raster Parameter Error from Python Script Tool?

texts:
I borrowed a simple script I found online to clip a shapefile from the extent of a raster. It works perfectly from the ArcMap built in python command prompt. I modified is slightly to use in a tool. It's having issues accepting the raster (inraster) in:

extent = arcpy.Raster(inRaster).extent

I get the following error:

  extent = arcpy.Raster(inRaster).extent TypeError: expected a raster or
  layer name

I'm using "Raster Layer" as my inRaster parameter type. The raster type I am using is an ERDAS IMAGINE .img file. I'm thinking it must be conflicting with the parameter type and the .img input raster, but I can't figure out what it is. As I mentioned, it worked from the command prompt, just not in the tool. Any suggestions are welcome.

Code is as follows:

import arcpy

inLines = arcpy.GetParameter(0)
inRaster = arcpy.GetParameter(1)
outRaster = arcpy.GetParameter(2)

pnt_array = arcpy.Array()
extent = arcpy.Raster(inRaster).extent
pnt_array.add(extent.lowerLeft)
pnt_array.add(extent.lowerRight)
pnt_array.add(extent.upperRight)
pnt_array.add(extent.upperLeft)

poly = arcpy.Polygon(pnt_array)

arcpy.Clip_analysis(inLines, poly, outRaster)


texts: Exporting Data Driven Pages with a display extent depndant feature layer

texts:
I am creating a topography map of many sites (200) around the UK in ArcMap 10.6. The map just has to show the site polygon, LiDAR and a legend.
The LiDAR is a stretched colour ramp which scales to the extent of the display map. As I flick through the sites in layout view, the LiDAR refreshes and scales appropriate for each site map. However, when I export to individual PDFs (using the standard DDP function), the LiDAR scale is fixed to the scale and range of the particular site I happened to be looking at the time. Meaning all sites have the same scale and range, and some therefore show no LiDAR data at all if they are not within the range of the "chosen one".
I am also running a script to export the DDP enabled mxd to JPEGs. This outputs something different from the PDFs. The JPEGs all share the same scale, but this is the scale of the entire data set (sea level to the highest point) meaning the relatively small changes in elevation at the site are not visible.
From looking around online, there's something called dynamic legends which I don't have access to that tier of license. Is it possible to have a refresh command between each export to refresh the legend and scale? and would that work? or any other method?


texts: Weird artefacts/pixels in QGIS-created shaded relief map

texts:
I'm working with 30m SRTM data on an terrain map of Cuba. I downloaded the necessary DEM files, merged them and went on to create the shaded relief map in QGIS 2.8 using Raster > Analyze > DEM (Terrain models).

I tried a bunch of different settings but I always get this same result:

Not finding what was wrong (I first thought it was a problem of converting degrees to meters), I double-checked the .tif file in Aperçu and it looks just fine there:

To clarify, those are the same file, loaded in QGIS as a raster and simply opened in OS X Aperçu. Any ideas what the source of the problem is in QGIS?

EDIT: Here's what the gdalinfo command returns about the image:

Driver: GTiff/GeoTIFF
Files: Cuba30Terrain.tif
Cuba30Terrain.tif.aux.xml
Size is 39601, 18001
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (-85.000138888888884,24.000138888888888)
Pixel Size = (0.000277777777778,-0.000277777777778)
Metadata:
  AREA_OR_POINT=Area
 Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( -85.0001389,  24.0001389) ( 85d 0' 0.50"W, 24d 0' 0.50"N)
Lower Left  ( -85.0001389,  18.9998611) ( 85d 0' 0.50"W, 18d59'59.50"N)
Upper Right ( -73.9998611,  24.0001389) ( 73d59'59.50"W, 24d 0' 0.50"N)
Lower Right ( -73.9998611,  18.9998611) ( 73d59'59.50"W, 18d59'59.50"N)
Center      ( -79.5000000,  21.5000000) ( 79d30' 0.00"W, 21d30' 0.00"N)
Band 1 Block=39601x1 Type=Byte, ColorInterp=Gray
  Min=1.000 Max=255.000
  Minimum=1.000, Maximum=255.000, Mean=180.682, StdDev=7.098
  NoData Value=0
  Metadata:
    STATISTICS_MAXIMUM=255
    STATISTICS_MEAN=180.68177132469
    STATISTICS_MINIMUM=1
    STATISTICS_STDDEV=7.0979989753955


texts: How to mass-create maps with different features centered?

texts:
Can anyone explain if this is something that can be done? I have map consisting of a little under 40,000 acres. I would like to print a map based on the legal location with each section being centered on a separate map at a scale of 1:1,000, for a total of 130 individual maps. Each map will have the same data, the only thing that will change from map to map is the feature that is centered.

Is it possible to have arc auto-create a map for each feature? I am using ArcMap 10.1 for desktop.


texts: Iterating over years for features in feature collection using Google Earth Engine?

texts:
I am currently working with the Hansen data (Global Forest Change) in Earth Engine. I have also imported a fusion table representing districts in a specific country. My goal is as follows:

I would like to create a table which aggregates deforestation data for each district in the feature table. For example, a given district ('feature') should eventually have the following associated properties: sum of gained forest area for each year, sum of lost forest area for each year, and base year tree cover. Once all of those are gathered, I plan to export the data.

I have figured out how to do this for one year, but the problem is I can't figure out how to iterate over years without doing everything manually. I am also not entirely sure whether my filtering actually gives me the year that I expect, as I copied the filtering from an example in the tutorials.

Here is some of my current code, which aggregates the tree cover for each district, and adds a 'loss' property to that feature collection:

//Load and filter the Hansen data
var gfc2014 = ee.Image('UMD/hansen/global_forest_change_2015').select(['treecover2000','loss','gain','lossyear']);

//add country districts as a feature collection
var distr = ee.FeatureCollection('ft:1U7sXFHXtxQ--g7XMeXlvPhNXPBcDtPg8Yzr2pvsg', 'geometry');

//look at tree cover, find the area
var treeCover = gfc2014.select(['treecover2000']);
var areaCover = treeCover.multiply(ee.Image.pixelArea());

var lossIn2001 = gfc2014.select(['lossyear']).eq(1);
var areaLoss = lossIn2001.multiply(ee.Image.pixelArea());

//need to check that .eq(1) is actually returning data for
//2001
var gainIn2001 = gfc2014.select(['gain']).eq(1);
var areaGain = gainIn2001.multiply(ee.Image.pixelArea());

var districtSums = areaCover.reduceRegions({
  collection: distr,
  reducer: ee.Reducer.sum(),
  scale: 100,
});

//This function computes pixels lost for each district
var addLoss = function(feature) {
  var loss = areaLoss.reduceRegion({
    reducer: ee.Reducer.sum(),
    geometry: feature.geometry(),
    scale: 100
  });
  return feature.set({loss2001: loss.get('lossyear')});
};

// Map the area getting function over the FeatureCollection.
var areaLosses = districtSums.map(addLoss);

I have looked to the following pages for help, but unfortunately neither appears to answer my question:

Using reducer over space and time in Google Earth Engine?
Computing ndvi 10 yr. time-series from fusion table features using Google Earth Engine?

This is my first time working with geographic data so the data structures are not intuitive to me yet.


texts: Using st_read_db in R to read postGIS rasters

texts:
I use the st_read_db function in the sf package to read vector data from PostGIS into R. Can i use it to query rasters? ideally extract a raster value at a point:

require(sf)
# create connection to PostGIS
conn <- dbConnect("PostgreSQL",user="postgres",password="password",port=5432,dbname="postgis")

# example using query for polygons - fine
qry <- "SELECT * FROM polygontable WHERE arg1 = 10"
pols <- st_read_db(conn, query=qry, geom="geom")

# the following SQL works directly in PostGIS
qryras <- "SELECT rid, ST_Value(rast, ST_SetSRID(ST_MakePoint(x,y),27700))
FROM rastertable WHERE ST_Intersects(rast, ST_SetSRID(ST_MakePoint(x,y),27700))"

# but does not work when passed to st_read_db
st_read_db(conn, query=qryras, geom="rast")

"Error in `[[<-.data.frame`(`*tmp*`, geom_column, value = list()) :
  replacement has 0 rows, data has 1"

With st_read_db being in the sf package, perhaps it only works for vector data?

*edit: this will work using dbGetQuery inside the RPostgreSQL package, but i specifically would like to know if sf will support what i want to do


texts: QGIS interactive gradient editor - already available?

texts:
i have found several links on the "new interactive gradient editor" - which seems pretty interesting and comfortable !!! the features shown in a youtube video are promising, but i didn't find any reference information there ...
does anyone know whether or not this tool is already available ?


texts: Esri Story Map access to everyone - data "endangered"?

texts:
If a story map is created that is accessible for everyone, the concerning web maps with my layers have to be published/shared with everyone before.

Doesn't this mean that the maps, layers, etc. also can be used by the whole ArcGIS online community for their web maps?

The reason why I ask is, that I would like to create a story map and provide access for everyone, but I do not want my data, which has been a huge effort to capture, to be used by everyone.


texts: mapbox.js 1.6.2 - change tile set on zoom levels

texts:
Alright, this is what I've got. What you've suggested doesn't seem to be working. I need to change the basemap at a certain zoom level...

    var map = L.mapbox.map('map', 'wmarci.i6n42nl5', {
      center: new L.LatLng(38.8929,-100.0252),
      zoom: 4,
      minZoom: 0, maxZoom: 5,
      layers: [ L.mapbox.tileLayer('wmarci.k300ifal', { maxZoom: 13, minZoom: 5, }) ]
    });

new L.Control.Zoom({ position: 'topright' }).addTo(map);

 //Navigate to Places
 document.getElementById('navigation').onclick = function(e) {
   var pos = e.target.getAttribute('data-position');
   if (pos) {
    var loc = pos.split(',');
    map.setView(loc, 10);
}
}

//Navigate to Places
document.getElementById('zoomout').onclick = function(e) {
var pos = e.target.getAttribute('data-position');
if (pos) {
    var loc = pos.split(',');
    map.setView(loc, 4);
}
}
    //Mapbox layers for switcher
    var ui = document.getElementById('map-ui-zoning');

 //  Layers and Grids!!!!!!!!
    var layers = [{
        name: 'basemap',
        minZoom: 0,
        maxZoom: 19,
          layer: L.layerGroup([
           L.mapbox.tileLayer('wmarci.i6n42nl5'),
           L.mapbox.gridLayer('')])
      },
      {
        name: '1775',
        minZoom: 0,
        maxZoom: 19,
        layer: L.layerGroup([
          L.mapbox.tileLayer('wmarci.g2ufpqfr'),
          L.mapbox.gridLayer('wmarci.g2ufpqfr')])
      },
      {
        name: '1805',
        minZoom: 0,
        maxZoom: 9,
        layer: L.layerGroup([
          L.mapbox.tileLayer('wmarci.tlk2csor'),
          L.mapbox.gridLayer('wmarci.tlk2csor')])
      },
      {
        name: '1835',
        layer: L.layerGroup([
          L.mapbox.tileLayer('wmarci.mt4t2o6r'),
          L.mapbox.gridLayer('wmarci.mt4t2o6r')])
      },
      {
        name: '1865',
        layer: L.layerGroup([
          L.mapbox.tileLayer('wmarci.m8peewmi'),
          L.mapbox.gridLayer('wmarci.m8peewmi')])
      },
      {
        name: '1895',
        layer: L.layerGroup([
          L.mapbox.tileLayer('wmarci.hafhto6r'),
          L.mapbox.gridLayer('wmarci.hafhto6r')])
      },
      {
        name: '2014',
        layer: L.layerGroup([
          L.mapbox.tileLayer('wmarci.bqz4obt9'),
          L.mapbox.gridLayer('wmarci.bqz4obt9')])
      },
      {
        name: 'contemp',
        minZoom:7,
        maxZoom:19,
        layer: L.layerGroup([
          L.mapbox.tileLayer('wmarci.k300ifal'),
          L.mapbox.gridLayer('wmarci.k300ifal')])
      }
    ];


texts: How to dynamically measure and display radius of circle in OpenLayers?

texts:
I want to give my user the ability to click on the map and give him a dynamic circle drawing tool which also reads out the radius measurement.  That way the user can see how big the circle he creates is.

For example, something like this :

Any idea how I can add such a thing to OpenLayers?


texts:
I have an OS dataset loaded into Geoserver as an 'ImageMosaic'. There are 10588 individual GeoTIFFs.

Geoserver obviously recognises and reads the GeoTIFFs, because if I try 'Layer Preview', it gives me a preview with the correct extents (both in OpenLayers and KML/Google Earth). However, there is no data displayed.

I have also tried loading a single GeoTIFF. The result is the same.

It should not be a problem with the GeoTIFFs themselves, as they open without any problem in QGIS etc., and they are official OS data, so you would assume they are good...

Does anyone know why my rasters are invisible?

GDAL info output:

Driver: GTiff/GeoTIFF
Files: TV69NW.tif
Size is 5000, 5000
Coordinate System is:
PROJCS["OSGB 1936 / British National Grid",
    GEOGCS["OSGB 1936",
        DATUM["OSGB_1936",
            SPHEROID["Airy 1830",6377563.396,299.3249646000044,
                AUTHORITY["EPSG","7001"]],
            TOWGS84[446.448,-125.157,542.06,0.15,0.247,0.842,-20.489],
            AUTHORITY["EPSG","6277"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4277"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",49],
    PARAMETER["central_meridian",-2],
    PARAMETER["scale_factor",0.9996012717],
    PARAMETER["false_easting",400000],
    PARAMETER["false_northing",-100000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","27700"]]
Origin = (560000.000000000000000,100000.000000000000000)
Pixel Size = (1.000000000000000,-1.000000000000000)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_COPYRIGHT=▒ Crown copyright and/or database right 2016 OS
  TIFFTAG_IMAGEDESCRIPTION=OS Open Map - Local - Full colour raster - TILE TV69NW
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
  TIFFTAG_XRESOLUTION=254
  TIFFTAG_YRESOLUTION=254
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  560000.000,  100000.000) (  0d16'10.78"E, 50d46'36.88"N)
Lower Left  (  560000.000,   95000.000) (  0d16' 2.96"E, 50d43'55.12"N)
Upper Right (  565000.000,  100000.000) (  0d20'25.89"E, 50d46'31.84"N)
Lower Right (  565000.000,   95000.000) (  0d20'17.82"E, 50d43'50.09"N)
Center      (  562500.000,   97500.000) (  0d18'14.36"E, 50d45'13.50"N)
Band 1 Block=5000x32 Type=Byte, ColorInterp=Palette
  Color Table (RGB with 256 entries)
    0: 0,0,0,255
    1: 8,7,7,255
    2: 24,24,24,255
    3: 40,40,40,255
    4: 52,52,52,255
    5: 85,15,45,255
etc...


texts: Why does OpenLayers render two AGSCache layers differently?

texts:
I am using OpenLayers to try to display aerial maps for two cities (McHenry, IL and Cupertino, CA).

By following the tutorial for OpenLayers I have gotten my McHenry example to render correctly, as shown in this jsfiddle: http://jsfiddle.net/je45Y/

I should be able to use the exact same code and change only the URL and the zoomToExtent() arguments at the end (such that they match the Initial Extent reported on the ESRI info page). However, when I do this, all the map tile requests generate 404 errors as shown in this jsfiddle: http://jsfiddle.net/9y90mvLo/

Eventually, if I zoom in enough, the tiles will display. This is pretty confusing to me. Presumably, there are different values for zoomToExtent() that I can provide to make this work properly right from the start, but I have no idea how I'd find out what these values are.

Is their ESRI server misconfigured somehow or does OpenLayers not work properly with this newer version of ESRI Server?

It's confusing that the same code works for the first server but not the second.


texts:
In my PostGIS table, I do not have any column with feature ID or FID.

But when I open it in GeoServer layer, It shows the FID

but I don't want to show it when i click over the polygon,
Is there any way to do that?
And even if I am writing the PostGres table with fid, it is not using this as fid

however I am using this column as Identifier in the Geoserver SQL view, but still getting the postgres generated fid.


texts: Protected general and admin settings in GeoODK app

texts:
I was wondering how to prevent the user to change the general and admin settings of the GeoODK (android-collect)? I mean to be protected with a password, for example.

I've tried to go to settings > admin settings > and unchecked whatever the user shouldn't access in the general settings tab, but still the "admin setting" tab is available and the user can get in there and modify anything!!


texts: Using raster iterator, failed raster to point conversion

texts:
I am trying to resample a gdb of rasters, extract them to my study area extent, and convert them to points.  The purpose of this is to have two raster data sets represented as points so I can export their tables and line the data up side by side for analysis. I am new to using the iterator in ModelBuilder.
Here is the model I am using,

However, I am getting the following error message:

  Possible empty feature class.

  Failed to execute (Raster to Point).

Does anyone know what I am missing?


texts:
From the file osm.pbf I created a shape that contains lines of a street network and attributes for 'lanes'.

Problem: The attribute 'lanes' is combined with many others in a column "other_tags".

After doing some research, I found that PostgreSQL and the hstore extension can filter individual attributes (documentation).

I want all OSM IDs and their corresponding values for "lanes".

Now:

SELECT *
FROM
  table
WHERE
  other_tags :: hstore -> 'lanes'

comes unknown operator.

Do I still have to set something specific to use hstore?


texts:
While working with DSAS 5.0 I could cast transects but while attempting for rate calculations an error displays 'unable to find database table for baseline layer'. The ArcMap 10.6 version is not crack and the entire feature class were created in personal geodatabase.
Can anyone guide me on this error?


texts:
Using ArcMap 10.2: when I'm exporting data-driven pages can I exclude single pages using some existing formatting rules (e.g. (-5) or [-5])? It would be a tad easier than my typical method of using ranges to exclude pages (e.g. 1-4, 6-10).

I realize I can set the definition query of the index layer to exclude the page to gain the same effect. I'm wondering if this type of option exists in the "Export" dialog.


texts: How to concatenate values of two fields with update cursor with ArcPy?

texts:
I have a small block of code where I concatenate two date fields and use field calculator to calculate a the values of a third field.  My date fields are ORIGDTDATE and ORIGDTTIME.

How can I change the block of code that is below to concatenate the fields if the value for the row in emergency_overdue is NULL?

import arcpy

calcExpression =  "!ORIGDTDATE!+\" \"+ !ORIGDTTIME!"

arcpy.CalculateField_management(self.line_ticket_feature_class,"emergency_overdue",calcExpression,"PYTHON_9.3")

arcpy.CalculateField_management(self.poly_ticket_feature_class,"emergency_overdue",calcExpression,"PYTHON_9.3")


texts: Generating random points inside raster boundary using R?

texts:
I have a projected raster shown in the image below (left).  I want to generate random points inside the raster like the middle image.  However, the result I am getting is the same as the image on the right side. How would I get the random points inside the extent of the projected raster?

library(rgdal)
library(raster)

a <- brick("E:/raster.tif")

b <- as(extent(a), "SpatialPolygons")

c <- spsample(b, n = 20, type = "random")

par(mar=c(1,1,1,1))
plot(c)


texts: Creating raster with distance to feature using QGIS?

texts:
I'm new to GIS (and QGIS in particular). I'm trying to create a raster layer that contains the distance to shoreline for a certain region.

I have both vector layers of the shoreline (line) and the region in question (polygon).

Is there a tool to create the raster using the line as the object to measure the distance and the polygon as a mask?


texts: Selecting max value from group using ArcGIS definition query?

texts:
In ArcGIS 10 I am trying to do a Definition Query that will select the records that are the max values from a specific group.

Here is what I have tried in the Query Builder but I am getting an "the SQL statement is invalid." error.

"TL" IN (SELECT max( "TL" ) FROM "Instrument_Locations2_20120227" GROUP BY "Site_ID" )


texts:
I have a set of field measurements for a bird species (pitch of its song) along with coordinates for each measurement.

The default choropleth style in CARTO works perfectly for showing how the pitch varies across its geographic distribution (see image; red means higher pitch). However, many of my data points overlap. How can I cluster them together and display the average pitch for overlapping points in choropleth-style?


texts: Display a layer stored in GeoServer with Leaflet

texts:
I wanted to display a layer stored in GeoServer with Leaflet, I used this code but it's not working.
var map = L.map('map').setView([4.477411, 21.309284], 3);

L.tileLayer('http://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png',
{
attribution: '&copy; OpenStreetMap'
})
addTo(map);

// Instantiation de différentes couches provenant de GeoServer

var exemple = L.tileLayer.wms("http://localhost:8080/geoserver/nurc/wms", {
    layers: 'nurc:Arc_Sample',
    format: 'image/png',
    transparent: true
});

exemple.addTo(map);

The error reported is:
Uncaught ReferenceError: L is not defined


texts: Preventing Project pane of ArcGIS Pro from re-confirming hundreds of Delete requests from multiple select?

texts:
Using ArcGIS Pro 1.2 ModelBuilder I just created a couple of hundred shapefiles in a folder to which I had a Folder Connection in the Project pane, and thought it would be very easy to use the Project pane to multiple select them and then right-click to Delete them.

The multiple-select part is easy, but when I right-click to Delete I am presented with the Delete dialog and asked:

  Are you sure that you want to permanently delete this item?

It is probably good that ArcGIS Pro double-checks before going ahead, but the problem is that it now requires me to re-confirm the deletion of a couple of hundred shapefiles one at a time, and offers me no way to bail out.

Is there an option to avoid this situation like configuring a checkbox on this dialog that says something like "Do this for the rest"?


texts: Multiplying values for each field with unique id

texts:
For one country I have

population densities (highlighted with a red rectangular) per watershed
literacy rate 4.48% (highlighted with a green rectangular)
watershed IDs, e.g. 3222, 3271 etc.

I'd like to multiply the population of each watershed by literacy rate. Is there a way to do so?

Perhaps using the field calculator?