texts: Updating Secure service on Portal for ArcGIS with PortalPy and requests?

texts:
I'm using the below Portalpy code to update my secure Feature Services on ESRI Portal which have a stored password.  The return on print post_url in the script is: https://somedomain.com/portal/sharing/rest/content/users/siteadmin/f43542fd62c74db6997a5d16c9eec16e/items/c3aa864f30ee4d8fbc5241134d35762c/update with a satus code of 200.  The above URL is identical when I do it manually and watch it with fiddler.  However the script is not changing the U/P in portal.  I assume I am passing the parameters incorrectly, but am following the user docs for construction and parameters.  Anyone have any suggestions?

Portal: 10.3v

import portalpy
import requests
import json

#Login Info
username = 'XXXX'
password = 'XXXX'
portal_url = 'https:/somedomain.com/portal'
portal = portalpy.Portal(portal_url, username, password)

#get token
token = portal.generate_token(username, password, expiration=60)

#set new username
new_username = 'newu'
new_password = 'newp'

payload = {'serviceUsername':new_username, 'servicePassword':new_password}

#print  portal.get_folder_id('siteadmin')

#Search for layers
search = portal.search('Feature Service')

for i in search:
    try:
        a = i['url'].split('/')
        folder_id = portal.get_folder_id('siteadmin', a[9])
        post_url = portal.resturl + 'content/users/siteadmin/' + folder_id + '/items/' + i['id'] + '/update'
        print post_url
        post = requests.post(post_url, data=payload)
        print post.status_code
    except Exception:
        pass


texts: Connect two points with fixed longitudinal slope over DEM - QGIS

texts:
Is it possible to connect two points with fixed longitudinal slope over DEM terrain or connect points with maximum longitudinal slope? Make road line between two point, for example 7% or make it with maximum slope, for example 10%?
Maybe somebody have some python code? (Qgis 2.14. essen)


texts: Accumulated cost surface (Tobler's hiking function) ignoring slope using R gdistance?

texts:
Background

I wish to calculate in R the isochrones moving from a location outwards, implementing the Tobler's hiking function via the gdistance package (https://cran.r-project.org/web/packages/gdistance/vignettes/gdistance1.pdf). An example of the implementation is provided at the above link, pp. 14-46; please note that in that reference, the issue is approached from a slighlty different perspective, i.e. least-cost path between two locations A and B.

Issue

I can't understand why I get two different results usign the same code over two  dataset (the one provided by the mentioned reference vs. my toy DTM). One result make intuitively sense, the other I believe is wrong.

Reproducible code

The code below is what I used for the two analyses. The first code uses a toy DTM which comes with the gdistance package itself. The dtm and point shapefile for the second analysis are available here:

https://www.dropbox.com/s/8t0fkuovehrsf0z/mydtm.rda?dl=0

https://www.dropbox.com/s/rlzclkzlv37v2v0/loc.rda?dl=0

Following the abovementioned source (with sligth modifications), I get the first (reasonable) result attached (1st and 2nd image):

dtm <- raster(system.file("external/maungawhau.grd", package="gdistance"))
altDiff <- function(x){x[2] - x[1]}
hd <- transition(dtm, altDiff, 8, symm=FALSE)
slope <- geoCorrection(hd)
adj <- adjacent(dtm, cells=1:ncell(dtm), pairs=TRUE, directions=8)
speed <- slope
speed[adj] <- 1/(6 * exp(-3.5 * abs(slope[adj] + 0.05)))
Conductance <- geoCorrection(speed)
A <- c(2667670, 6479000)
accum <- accCost(Conductance, A)
plot(dtm)
plot(accum)
points(A[1], A[2], pch=20)

At the best of my understanding, the accum raster is (correctly, I think) representing the time it takes to move from the source location. Side note: I used the inverse of the Tobler's function to get time not speed).

Different story with my DTM (3rd and 4th image attached):

dtm <- mydtm
altDiff <- function(x){x[2] - x[1]}
hd <- transition(dtm, altDiff, 8, symm=FALSE)
slope <- geoCorrection(hd)
adj <- adjacent(dtm, cells=1:ncell(dtm), pairs=TRUE, directions=8)
speed <- slope
speed[adj] <- 1/(6 * exp(-3.5 * abs(slope[adj] + 0.05)))
Conductance <- geoCorrection(speed)
A <- cbind(loc@coords[,1], loc@coords[,2])
accum <- accCost(Conductance, A)
plot(dtm)
plot(accum)
plot(loc, add=TRUE, pch=20)

Question

At the best of my understanding, from the "exagonal-ish" outlook of raster, I have the feeling that something is going wrong in the process of accounting for the slope. I have seen something like this in case of isotropic accumulated cost-surfaces, and I do not expect to see this since the above code should implement on anisotropic cost-surface analysis (as correctly shown using the toy dtm which comes with the package). I keep getting the same "strange" result even if I re-run the analysis using a larger extent of the dtm (featuring a large variation in the elevation values). Can anyone please help me in locating what I am doing wrong?


texts: Adding heatmap legend to print composer in QGIS

texts: I have created a heatmap from points layer's style property. However when I add a legend to the print composer I do not get anything in the legend except the layer name, how can I add the color ramp besides the layer name?


texts: QGIS 3.10.1 and 3.6 don't pull to front when clicked on taskbar?

texts:
There's a similar question already (QGIS 3.4 won't come to the foreground), but it's been left unanswered and I don't have enough reputation to add a bounty, so please don't yell at me.

Like others in the linked thread, I think it's possible that it started after running OSGeo4W. I tried installing different versions (same result), through OSGeo4W and by using standalone installers. I got rid of OSGeo4W off my disc (just the main file), but it didn't do anything. There's also this cmake file that's included in every version and OSGeo4W? The only version that works fine on my laptop (3.2, installed long before that) does not have that file, so I removed it from 3.10.1 but it didn't help.

It's extremely annoying to have to minimise windows all the time.


texts: Field Error in ArcPy TableToDbase conversion

texts:
I am trying to convert a .csv file in Python/ArcPy to dBase format.

I am using this code:

import arcpy
inputcsv = "D:/2016_08_21.csv"
outputfolder = "D:/"
arcpy.TableToDBASE_conversion(inputcsv, outputfolder)

That works but the problem is that the first field in the CSV ("DWD_ID") does not always contain integers. At first the IDs start like 10321 and so on but eventually there are letters too. However the Conversion tool recognizes the first value in the field as an integer and sets the field as integer type and later on the other values are kicked out and never appear in the .DBF.

Is there a better way to convert a .csv to .dbf or can I explicitly define the field "DWD_ID" as Text?

I have tried the DBF module of python but I have not been able to successfully convert a csv file to dBase. My output file was always empty. I am only seeking programmatic solutions in Python; I do not want to use any other software.


texts: Identifying user triggering ArcGIS Server invalid username or password error

texts: Can anyone tell me if there is a way to identify the culprit of an "Invalid Username or Password" error in ArcGIS Server logs? For source, it just says 'Admin' for source, and nothing for User. Is there a .txt log stored somewhere that would tell me which username is causing the bad login?


texts: Creating Projection/Datum TransvMerc/OrdSurvGrtBrit in QGIS?

texts:
I have Memory Map files QCT exported to ECW.

The Projection/Datum is Transverse Mercator / Ordnance Survey Great Britain which is not mentioned in the list of Projections/Datums available.

How can I create it from the scratch or by modifying any available?


texts: Creating new QGIS color ramp from r,g,b text file?

texts:
I am using QGIS 2.18.22

I would like to create the Cividis color ramp in QGIS.  Here's a (non-cartographic) display example of Cividis:

Source

Cividis' creators have made public a 256-line r,g,b text file defining the color ramp from dark to light.  That file takes the form of:

0,32,77

0,33,78

0,34,80

...

How do I create the Cividis color ramp in QGIS from the 256-line r,g,b text file?

I have tried Settings > Style Manager > Color Ramps > Import Symbol(s), but it asks for an xml file, which I don't know anything about.


texts: How to replace the legacy function "doQueries " to AMD ESRI Javascript format?

texts:
I'm trying to understand the migration of ESRI API for javascript, from legacy to AMD. In the case of this example, how replace the  fucntion doqueries????:

function init() {
    map = new esri.Map("map", {
      extent: new esri.geometry.Extent(-102.61, 36.2, -93.82, 40.5, new esri.SpatialReference({wkid: 4326})),
      slider: false
    });
    dojo.connect(map, "onLoad", doQueries);
    map.addLayer(new esri.layers.ArcGISTiledMapServiceLayer("http://server.arcgisonline.com/ArcGIS/rest/services/ESRI_StreetMap_World_2D/MapServer"));
  }

  function doQueries(map) {
    //Query all counties in Kansas
    var countyQueryTask = new esri.tasks.QueryTask("http://sampleserver1.arcgisonline.com/ArcGIS/rest/services/Demographics/ESRI_Census_USA/MapServer/3");
    var countyQuery = new esri.tasks.Query();
    countyQuery.outFields = ["*"];
    countyQuery.returnGeometry = true;
    countyQuery.outSpatialReference = map.spatialReference;
    countyQuery.where = "STATE_NAME = 'Kansas'";
    countyQueryTask.execute(countyQuery, addCountyFeatureSetToMap);

    //Query all cities in Kansas
    var cityQueryTask = new esri.tasks.QueryTask("http://sampleserver1.arcgisonline.com/ArcGIS/rest/services/Specialty/ESRI_StatesCitiesRivers_USA/MapServer/0/");
    var cityQuery = new esri.tasks.Query();
    cityQuery.outFields = ["*"];
    cityQuery.returnGeometry = true;
    cityQuery.outSpatialReference = map.spatialReference;
    cityQuery.where = "STATE_NAME = 'Kansas'";
    cityQueryTask.execute(cityQuery, addCityFeatureSetToMap);
  }

  function addCountyFeatureSetToMap(featureSet) {
    var symbol = new esri.symbol.SimpleFillSymbol();
    symbol.setColor(new dojo.Color([150,150,150,0.5]));

    //Create graphics layer for counties
    var countyLayer = new esri.layers.GraphicsLayer();
    map.addLayer(countyLayer);

    var infoTemplate = new esri.InfoTemplate("${NAME}","${*}");

    //Add counties to the graphics layer
    dojo.forEach(featureSet.features, function(feature) {
      countyLayer.add(feature.setSymbol(symbol).setInfoTemplate(infoTemplate));
    });
  }

  function addCityFeatureSetToMap(featureSet) {
    var symbol = new esri.symbol.SimpleMarkerSymbol();
    symbol.setColor(new dojo.Color([0,0,255]));

    //Create graphics layer for cities
    var cityLayer = new esri.layers.GraphicsLayer();
    map.addLayer(cityLayer);
    map.reorderLayer(cityLayer,1);

    var infoTemplate = new esri.InfoTemplate("${CITY_NAME}","${*}");

    //Add cities to the graphics layer
    dojo.forEach(featureSet.features, function(feature) {
      cityLayer.add(feature.setSymbol(symbol).setInfoTemplate(infoTemplate));
    });
  }

  dojo.addOnLoad(init);


texts: GEE Change the number of pixels on chart

texts:
I'm trying to create a time series chart that will tell me that mean precipitation for a given imageCollection, but my polygon is too big and I get the next error message:

  "Error generating chart: Collection.first: Error in map(ID=19991001):
  Image.reduceRegion: Too many pixels in the region. Found 67144509, but
  only 10000000 allowed."

How can I change the number of pixels and allow it to calculate the mean precipitation value for a given polygon? Whenever I run my code with the big polygon I get the error messege.

This is my code:

var geometry=table;
var dataset = ee.ImageCollection('NOAA/PERSIANN-CDR')
                  .filter(ee.Filter.date('1999-10-01', '2018-03-31'))
                  .filterBounds(geometry);

var precipitation = dataset.select('precipitation');

var onlyPercipitation = precipitation.filter(ee.Filter.calendarRange(10, 3,'month'));
print(onlyPercipitation);

print(ui.Chart.image.series(onlyPercipitation, geometry, ee.Reducer.mean(),30));


texts: Is there a way to allow a FeatureCollectionTable to update the geometry of existing features?

texts:
I am using the ArcGIS Runtime SDK for .NET version 100.2. I am new to the SDK in particular and GIS in general. My goal is to take data from a locally available database with associated lat/long and show it on a map. When the database is updated with a new location, I get a push notification and I would like to update that feature on the map.

This is a (contrived) example of what I would like to do.

private async Task<Map> BuildAnUpdateableMap()
{
    var map = new Map(Basemap.CreateTopographic());
    var fields = new[] { new Field(FieldType.Text, "foo", "bar", 10) };
    var table = new FeatureCollectionTable(fields, GeometryType.Point, SpatialReferences.Wgs84);
    var collection = new FeatureCollection(new[] { table });
    var layer = new FeatureCollectionLayer(collection);
    var feature = table.CreateFeature();
    feature.SetAttributeValue(fields[0], "baz");
    feature.Geometry = new MapPoint(-119, 36);

    // Add the feature to the table. This works fine.
    await table.AddFeatureAsync(feature);
    map.OperationalLayers.Add(layer);

    // Update feature geometry.
    feature.Geometry = new MapPoint(-119, 30);
    // Push the updated feature. This throws an exception.
    await table.UpdateFeatureAsync(feature);

    return map;
}

The initial feature and geometry shows up just fine but when trying to update the geometry, an exception is thrown:

Esri.ArcGISRuntime.ArcGISRuntimeException: 'Ownership-based access control does not permit this operation on this feature.: allowGeometryUpdates=false'

I can accomplish the update by first removing the feature from the table and then adding it again after updating the geometry but that results in an unsightly flicker (in the WPF control at least) and I expect that wouldn't scale well either.

Is there a way to update feature geometry with FeatureCollectionTable? Or do I need to use something different?


texts: Displaying rasters with transparent background in tilemill

texts:
I have created a set of annotations for different web zoom levels in ArcMap, and would like to use these labels in Tilemill rather than using the tilemill labeling engine.

I am aware that I can't import annotations to tilemill as if they are features, so what I would like to do is export each zoom level's labels as a raster at scale (as high resolution geotiff with a transparent background), and overlay these on each zoom in my tilemill project.

I can bring the label rasters into tilemill no problem and they show up in the right place, however I can't figure out how to make a transparent background work.

I have explored some of the raster-colorizer methods with no luck, anyone have any other solution to this?


texts: Comparing projections of two datasets when SRID not available

texts:
Similar to this question:
How to test if two GDAL datasets are in the same projection?

However, the above does not cover cases where an SRID is not available.
It can often be the case that an SRID is not specified and the WKT between 2 datasets of the same projection does not precisely match.

What is your strategy for determining whether 2 datasets are of the same projection (using python) in this case?

It is annoying to go through the whole process of reprojecting a dataset just because the WKT strings do not exactly match, even though they are in the same projection.

osr's Spatial Reference class can allow us toe get the projection information in a number of formats, and I often find that comparing the results of ExportToProj4() rather than ExportToWkt() to be more reliable in matching projections.
However, I do not know if it is 'bullet-proof'?

Another method could be to parse the WKT strings and match the 'important' parts. This has 2 questions:
1. Does anyone know an efficient method (in python) of breaking a WKT string into it's constituent parts which can be easily compared?
2. What are the relevant parts of a WKT string, or how can you determine this?


texts:
I get the following error each time I run the raster2pgsql:
[postgres@li1959-236 ~]$ raster2pgsql -I -C -s 4326 -M -F /home/temp/wrf_d02_202005010000.tif grib_test |  psql -d cep
Processing 1/1: /home/temp/wrf_d02_202005010000.tif
BEGIN
CREATE TABLE
free(): invalid pointer
ERROR:  unterminated quoted string at or near "'

I converted a band of a Grib into a GeoTIFF and I did nothing else.
My PG version is the 12 and posting 3.0.
Any suggestion?
# gdalinfo /home/temp/wrf_d02_202005010000.tif |more
Driver: GTiff/GeoTIFF
Files: /home/temp/wrf_d02_202005010000.tif
Size is 393, 348
Coordinate System is:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
Origin = (2762460.472548419144005,3366910.845273251645267)
Pixel Size = (7200.000000000000000,-7200.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( 2762460.473, 3366910.845) (Invalid angle,Invalid angle)
Lower Left  ( 2762460.473,  861310.845) (Invalid angle,Invalid angle)
Upper Right ( 5592060.473, 3366910.845) (Invalid angle,Invalid angle)
Lower Right ( 5592060.473,  861310.845) (Invalid angle,Invalid angle)
Center      ( 4177260.473, 2114110.845) (Invalid angle,Invalid angle)
Band 1 Block=393x2 Type=Float64, ColorInterp=Gray
  Description = 2[m] HTGL (Specified height level above ground)
  Metadata:
    GRIB_COMMENT=Temperature [C]
    GRIB_ELEMENT=TMP
    GRIB_FORECAST_SECONDS=259200 sec
    GRIB_REF_TIME=1588032000 sec UTC
    GRIB_SHORT_NAME=2-HTGL
    GRIB_UNIT=[C]
    GRIB_VALID_TIME=1588291200 sec UTC


texts:
I am using ArcMap 10.7.

I am somewhat new to geospatial data analysis and have an analysis question that I am hoping the community here might be able to assist with. I am working with a large raster file (NLCD 2016 CONUS) covering the contiguous United States (file type == '.img').

My object is (theoretically) simple: I want to use a polygon shapefile of the lower 48 states to clip (my initial guess for the appropriate tool, but open to other suggestions) the raster so that I have pixel information for each state in the attribute table. This would allow me to more easily perform downstream statistical summaries of the data that are state-specific.

My initial approach was to use the Data Management Tools --> Raster --> Raster Processing --> Clip function with a shapefile of the lower 48 as the Output Extent. I did choose ouput raster format as '.tif'; not sure if this is problematic.

However, after nearly 8 hours of processing time, the resulting clipped raster did not have an associated attribute table and did not retain the land cover pixel information from the original file.

I am hoping for some general feedback on this issue and modifications to my approach (maybe clip isn't the best tool...) that might be more successful. I am exploring other tools, but the enormous size of the dataset limits my ability to experiment with several different tools.


texts: Reloading own lib in QGIS python console?

texts:
I´m working on a main script with additional code based in a lib folder like this:

main.py
\mylib
  util.py (contains the class Utility)
  node.py
in main.py I´ve imported mylib like this:

from mylib.util import Utility
x = Utility()

It seems like QGIS python console caches the mylib on startup - changes made to the Utility class made after startup is ignored.

Is there a way to reload mylib without restarting QGIS?


texts:
How to load QGIS Python tools in ArcMap?

For example, if I want to use any GDAL or GRASS tools in ArcMap.


texts: openlayers 2.13.1 vector layer and CQL does not work in IE11

texts:
I have implemented a search functionality for a vector layer in openlayers 2.13.1. I use CQL.

I take the values from an html form and I use javascript to pass them to a CQL filter and execute it.

Works fine in Chrome  31.0.1650.63  and Firefox 25.0.1. But not in Internet Explorer 11. Was working in IE10, not in IE 11.

I search for some features and I know that should show features, but nothing shows on the map.

I hit map.layers[2].features on IE11's console and I get [functions], __proto__, layer, renderIntent and length is 1 but, no fid, no geometry.

Also, if I hit refresh in the page, in IE11 console I get The object does not support the property or method 'destroy' (translated from Greek, see screenshot).

Here is my code

    //fallback html5 websockets with flash -libraries
    <script type="text/javascript" src="swfobject.js"></script>
    <script type="text/javascript" src="web_socket.js"></script>

    <script src="jquery-1.8.2.min.js"></script>

    //polyfiller for html5 elements
    <script src="js-webshim/minified/extras/modernizr-custom.js"></script>
    <script src="js-webshim/minified/polyfiller.js"></script>
        <script>
          $.webshims.polyfill();
        </script>

    //polyfiller for css3 media queries
    <script src="../respond.src.js"></script>

    <script src="http://maps.google.com/maps/api/js?sensor=false"></script>

    <script type='text/javascript' src='OpenLayers.js'></script>

    <script type='text/javascript'>

    function init(){//called with body onlad, inits map, layers etc.
      //init the layer
      pins = new OpenLayers.Layer.Vector("Pins", {renderers: ["Canvas", "SVG", "VML"]})

      //set index
      map.setLayerIndex(pins, 70);

      //set styles for pins layer
      var vector_style_map = new OpenLayers.StyleMap({'default': style_d,'select': style_select});
      pins.styleMap = vector_style_map;

    }

   //called when the search button of the form is clicked
   function zoomSearch(){

      //"clean" layer first
      pins.removeFeatures(pins.features);
      pins.refresh({force:true});

      //set array
      var filter_list = [];

      //take the values
      if (document.getElementById("nm").value != "") {
        var filter_c = new OpenLayers.Filter.Comparison({
            type: OpenLayers.Filter.Comparison.LIKE,
            property: 'p_name',
            value: '%'+document.getElementById("nm").value+'%'
        });
        filter_list.push(filter_c);}

      if (document.getElementById("descr").value != "") {
        var filter_d = new OpenLayers.Filter.Comparison({
            type: OpenLayers.Filter.Comparison.LIKE,
            property: 'p_descr',
            value: '%'+document.getElementById("descr").value+'%'
        });
        filter_list.push(filter_d);}

      //create the final filter,combining the above,  to search
      var filter_final = new OpenLayers.Filter.Logical({
        type: OpenLayers.Filter.Logical.AND,
        filters: filter_list});

      //pass the filter to the layer
      var prot =  new OpenLayers.Protocol.WFS({
        url:  "http://localhost/geoserver/wfs",
        featureType: "pins",
        featureNS: "http://www.mysite.gr",
        defaultFilter: filter_final});

      var _CallBack = function(resp) {
          pins.addFeatures(resp.features)
         var cb = pins.features.length;
          if (cb == 0){alert"nothing. try again";}
          if(cb==1){map.zoomToExtent(pins.getDataExtent());   map.zoomTo(8);}
         else{map.zoomToExtent(pins.getDataExtent());}
        };

      var response = prot.read({callback: _CallBack});

      //refresh to render POIs
      pins.refresh({force:true});

    }

How can I fix this? Is it me  or is it a bug?

Thanks in advance


texts: Inserting Polygon data into PostGIS using Python

texts:
I am able to insert polygon data manually into a postGres DB with one column having polygon datatype however, when I try to do the same dynamically it is not working.
This works.
cur.execute("INSERT into  public.test(timestamp,poly) values (now(),'((2,2),(3,4),(3,6),(1,1))');")

This doesnt.
    coords='((2,2),(3,4),(3,6),(1,1))'
    cur = cnn.cursor()
    cur.execute("INSERT into  public.test(timestamp,poly) values (now(),'coords');")

so only difference is that , this time I am converting the same co ordinates into a tuple of tuples and passing that as a variable into the query.

Error: psycopg2.errors.InternalError_: parse error - invalid geometry


texts:
I need help with some python in selecting correct input and near feature classes in Near Table analysis.

The project aim is to run Near Table Analysis (distance) on point locations collected at time 1[day] and time2[night].

I have managed to group the files I want to input into each iteration of the model based on filename and Year using groupby with the first 0:26 characters in the file name (all files are in the same directory).

My problem is I need to separate the inputs for the GenerateNearTable Tool into in_features (Day) and near_features (Night). So split each group again into Day and Night files.

Here is my code.

import arcpy, itertools, os

# Define the input and output workspaces
arcpy.env.workspace = r"C:\Users\Jenny\roos"
outws = r"C:\Users\Jenny\roos\gmehome.gdb"

# List the fcs in your workspace
my_pts = arcpy.ListFeatureClasses()

# Group the items based on first 26 characters
grouped = [list(g) for _, g in itertools.groupby(sorted(my_pts), lambda x: x[0:26])]

# Loop over the groups and perform the NearTable operation
for group in grouped:
    groupedDay = ?
    groupedNight = ?
    outname = group[0][:26] + "_DayNightDist"
    arcpy.GenerateNearTable_analysis(in_features=groupedDay, near_features=groupedNight, out_table= outname, search_radius="", location="NO_LOCATION", angle="NO_ANGLE", closest="ALL", closest_count="0", method="PLANAR")

Other useful information is the file names

  [[Jenxxx_Annual_HR_Poly_Yr1Day_SCV,
  Jenxxx_Annual_HR_Poly_Yr1Night_SCV]
  [Jenxxx_Annual_HR_Poly_Yr2Day_SCV,
  Jenxxx_Annual_HR_Poly_Yr2Night_SCV], [Barry_Annual_HR_Poly_Yr1Day_SCV,
  Barry_Annual_HR_Poly_Yr1Night_SCV]] (and so on)

UPDATE: SCREENSHOT OF ITERATION ISSUE


texts: Fragmentation in multicoloum geometry data type

texts: I have a table with place-identifier(big int) and geometry column combined together as b-tree unique index in PostgreSQL. Upon insertion of 20k records the dB table is getting fragmented around 45-50%. Not sure how I apply Multicolumn geometry index and avoid fragmentation. Any information?


texts: Exporting classification error matrix Google Earth Engine

texts:
I did a supervised classification of land cover and I need to export the Validation error matrix to .csv file. When using 'Export.table.toDrive' an error message 'Invalid argument: 'collection' must be a FeatureCollection' appears. Is there any way how to do that?

var year99=ee.Image('LANDSAT/LE07/C01/T1_SR/LE07_191025_19990915')
.clip(table)
.updateMask(FinalMask);

// Input image

var input=year99;

// Select and visualize MODIS land cover, IGBP classification, for training.

var modis=ee.Image('MODIS/051/MCD12Q1/2013_01_01')
.select('Land_Cover_Type_1')
.clip(table)
.updateMask(FinalMask);

// Define a palette for the 18 distinct land cover classes.

var igbpPalette = [
  'aec3d4', // water
  '152106', '225129', '369b47', '30eb5b', '387242', // forest
  '6a2325', 'c3aa69', 'b76031', 'd9903d', '91af40',  // shrub, grass
  '111149', // wetlands
  'cdb33b', // croplands
  'cc0013', // urban
  '33280d', // crop mosaic
  'd7cdcc', // snow and ice
  'f7e084', // barren
  '6f6f6f'  // tundra
];

// Specify the min and max labels and the color palette matching the labels.

Map.setCenter(15.293, 50.659, 8);

Map.addLayer(modis,
             {min: 0, max: 17, palette: igbpPalette},
             'IGBP classification');

// Visualize input image

var visParams = {bands: ['B3', 'B2', 'B1'], min: 16, max: 1412};

Map.addLayer(input,visParams, 'input');

// Print input image and Modis Landcover image

print ('input',input);
print ('modis',modis);

// Sample the input imagery to get a FeatureCollection of training data.
var training = input.addBands(modis).sample({
  numPixels: 500,
  seed: 0,
  region: table
});

// Make a Random Forest classifier and train it.
var classifier = ee.Classifier.randomForest(10)
    .train(training, 'Land_Cover_Type_1');

// Classify the input imagery.
var classified = input.classify(classifier);

// Get a confusion matrix representing resubstitution accuracy.
var trainAccuracy = classifier.confusionMatrix();
print('Resubstitution error matrix: ', trainAccuracy);
print('Training overall accuracy: ', trainAccuracy.accuracy());

// Sample the input with a different random seed to get validation data.
var validation = input.addBands(modis).sample({
  numPixels: 500,
  seed: 1,
  region:table
  // Filter the result to get rid of any null pixels.
}).filter(ee.Filter.neq('B1', null));

// Classify the validation data.
var validated = validation.classify(classifier);

// Get a confusion matrix representing expected accuracy.
var testAccuracy = validated.errorMatrix('Land_Cover_Type_1',         'classification');
print('Validation error matrix: ', testAccuracy);
print('Validation overall accuracy: ', testAccuracy.accuracy());

// Display the input and the classification.
Map.centerObject(table, 8);
Map.addLayer(input, {bands: ['B3', 'B2', 'B1'], max: 0.4}, 'landsat');
Map.addLayer(classified, {palette: igbpPalette, min: 0, max: 17},     'classification');

print('classified',classified);


texts: Definition Query / SQL problem in ArcGIS

texts:
I have been working on some Definition Queries for layers in ArcGIS (10.2.0.3348). The layers connect to data in a File Geodatabase. Now I have attempted to implement the following query (1), that returns 3051 records when I open the attribute table. However, when I run queries 2 and 3, I get 3280 and 2 records respectively. I would have expected the difference in the number or records between query 1 and 2 to be 2 records, as per the result of query 3, but obviously, this is not the case. I guess my understanding of the processing of these queries by the database engine is wrong, and especially the "AND NOT" part of it, but can someone give a better insight why these 3 queries don't match up in terms of number of records?

Also note the result of query 4: it seems the database engine ignores the extra statements against osm_level, osm_layer and  osm_location, as the number of records in query 4 is equal to the difference in query 1 and 2?

What am I missing? I have also tried to restructure the query to the one under 5, but that returns the same number of records as query 1...

I realize I could put the required values for "railway" in the IN statement, instead of trying to exclude with the "AND NOT", but I have some specific reason to maintain this structure (catch all for unknown values that get symbolized by the "all other values" item in an ArcGIS layer), and need to live within the limitations of the Definition Query and Layer framework of ArcGIS.

(I would have loved to post this question on the new GeoNet as well, but just can't work with all the clutter...)

QUERY 1:

  (railway NOT IN (Null,'') AND osmSupportingElement = 'no') AND NOT
  (railway IN ('platform') AND (osm_level LIKE ('%-%') OR osm_layer LIKE
  ('%-%') OR osm_location IN ('underground','underwater')))

Number of returned records: 3051

QUERY 2:

  railway NOT IN (Null,'') AND osmSupportingElement = 'no'

Number of returned records: 3280

QUERY 3:

  (railway IN ('platform') AND (osm_level LIKE ('%-%') OR osm_layer LIKE
  ('%-%') OR osm_location IN ('underground','underwater')))

Number of returned records: 2

QUERY 4:

  railway IN ('platform')

Number of returned records: 229

QUERY 5:

  ((railway NOT IN (Null,'') AND osmSupportingElement = 'no') AND NOT
  ((railway IN ('platform') AND osm_level LIKE ('%-%')) OR (railway IN
  ('platform') AND osm_layer LIKE ('%-%')) OR (railway IN ('platform')
  AND osm_location IN ('underground','underwater'))))

Number of returned records: 3051


texts: Creating contours using bare-earth LiDAR .XYZI point data

texts: I have 42 tiles of LiDAR .XYZI point data containing approximately 800,000 points each and I need to create 1 foot interval contours out of the data. I have access to ArcMap 10.3 and the 3D Analyst and Spatial Analyst toolbar's. The final goal of the project would be to display the 1' contours associated with flow lines.


texts: WKID (Well Known ID) for national projection for ArcGIS REST API

texts: I have to use my national projection system which does not have corresponding WKID listed for ArcGIS REST API.   Without it I am unable to use this projection system for my web application.  Any pointers which can help me out to use my data in its native projection format?  Thanks


texts:
I try to show the field names (alias names) of a feature class. I found code snippet but I don't succeed to get it to run :(
Can somebody tell me please what is wrong in the code? (System: Arcgis 10.2, Visual Studio2012)

Imports ESRI.ArcGIS.esriSystem
Imports ESRI.ArcGIS.Carto
Imports ESRI.ArcGIS.ArcMapUI
Imports ESRI.ArcGIS.Geodatabase

Public Class FieldNameButton
    Inherits ESRI.ArcGIS.Desktop.AddIns.Button

    Public Sub New()
    End Sub

    Protected Overrides Sub OnClick()
        ShowDistinctFieldAliasNames(My.ArcMap.Document)
    End Sub

    Shared Sub ShowDistinctFieldAliasNames(ByVal featureClass As IFeatureClass)

        ' Get the Fields collection from the feature class.
        Dim fields As IFields = featureClass.Fields
        Dim field As IField = Nothing

        ' On a zero-based index, iterate through the fields in the collection.
        For i As Integer = 0 To fields.FieldCount - 1
            ' Get the field at the given index.
            field = fields.Field(i)
            If field.Name <> field.AliasName Then
                Console.WriteLine("{0} : {1}", field.Name, field.AliasName)
                MsgBox(field.AliasName)
            End If
        Next i
        My.ArcMap.Application.CurrentTool = Nothing
    End Sub

    Protected Overrides Sub OnUpdate()
        Enabled = My.ArcMap.Application IsNot Nothing
    End Sub

    Private Sub DisplayDistinctFieldAliasNames()
        Throw New NotImplementedException
    End Sub

End Class


texts:
I have a dataset on some villages like this on a grid:

I want to fill the space of the whole grid so that each grid cell has a value that reflects the pressure from the surrounding villages (i.e., the value would be something like average distance to all villages, weighted by the size of each village).

As a first attempt, I tried to smooth it with:

    villageRaster <- focal(villageRaster, w=matrix(1, 9, 9), mean)

to give:

But there are still large white areas with no values.

Can someone suggest a nice way to smooth across the whole space?

I am using R for this preferably.


texts: Wanted to ask about accuracy of areas calculated by GIS compared to actual ones on ground. Because I have data of my country without areas for the districts within Uganda. Are the areas obtained using GIS accurate and do they reflect the actual areas on ground or in the field?


texts:
I am interested in learning how to write standalone scripts using QGIS. An answer to my previous question directed me to the PyGIS CookBook which has the following instructions for importing the library in python:

from qgis.core import *

I am using a Mac with OSX 10.7.4. To set the library and pythonpath, I used the following commands in my terminal, based on advice on this mailing list:

export DYLD_LIBRARY_PATH=/Applications/Qgis.app/Contents/MacOS/lib:/Applications/Qgis.app/Contents/Frameworks
export PYTHONPATH=/Applications/Qgis.app/Contents/Resources/python

Then, when I started python and typed from qgis.core import * I got the following error message:

python: VERSIONER_PYTHON_VERSION environment variable error (ignored)
Python 2.6.7 (r267:88850, Jul 31 2011, 19:30:54)
[GCC 4.2.1 (Based on Apple Inc. build 5658) (LLVM build 2335.15.00)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> from qgis.core import *
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ImportError: dlopen(/Applications/Qgis.app/Contents/Resources/python/qgis/core.so, 2): Symbol not found: _PyCapsule_Type
  Referenced from: /Applications/Qgis.app/Contents/Resources/python/qgis/core.so
  Expected in: flat namespace
 in /Applications/Qgis.app/Contents/Resources/python/qgis/core.so

Am I setting the paths correctly? I visually checked these and they appear to make sense based on what is in the folder, but I'm not familiar enough with this to know if there is a problem or not.

Edit:
In response to @lynxlynxlynx inquiring if a more recent version of python solved the problem, I got a different error message.

When using QGIS, I checked the python version and saw the following message:

>>>sys.version_info
sys.version_info(major=2, minor=7, micro=1, releaselevel='final', serial=0)

In my terminal I get the following message:

>>> sys.version_info
sys.version_info(major=2, minor=7, micro=2, releaselevel='final', serial=0)

If I trying using Python 2.7.2, I get a different error message:

Python 2.7.2 |EPD_free 7.2-2 (32-bit)| (default, Sep  7 2011, 09:16:50)
[GCC 4.0.1 (Apple Inc. build 5493)] on darwin
Type "packages", "demo", "upgrade" or "enthought" for more information.
>>> from qgis.core import *
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ImportError: dlopen(/Applications/Qgis.app/Contents/Resources/python/qgis/core.so, 2): no suitable image found.  Did find:
    /Applications/Qgis.app/Contents/Resources/python/qgis/core.so: mach-o, but wrong architecture


texts: Geoserver netcdf plugin to publish nc file - raster values wrong

texts:
the following image shows the geoserver i usedI used netcdf plugin to publish a nc file. I use wms to access it (see following picture for result), but I found that the value of the raster is wrong.

Another picture of the same nc file opened by the PanoplyWin (a pro nc software):

This is the structure of the nc file:

Visualization of mean Daily relative humidity given a specific time and a specific level. The mean daily relative humidity value is between 0 and 100, but the as shown above an example value is -26891, out of this range.


texts:
I am trying to merge quarter quarter sections of polygons from PLSS into a quarter section.
However, I also want it assigned to correct Township/Range/Section while merging.
I know I can hand merge all of them using the Editing tool, but since I will have to do to the entire state of Kansas.
I was thinking of writing a script to help.  The problem is that the Township/Range/Section field in the attribute table is in one field and I don't know how to do that when is all together.

I will need to merge those four polygon that belong to the same Township/Range/Sec/qtr as following:
01S01E01NWNE
01S01E01NWNW
01S01E01NWSE
01S01E01NWSW

and I need to leave in 01S01E01NWNW.

Thanks!


texts: Cropping a raster by a polygon - cells missing that are partially outside the polygon

texts:
I am trying to crop a raster in RStudio using the crop() function in the Raster package.

I have used the following code:

studyarea<-crop(chlorophyll,mcp2015)

Where 'chlorophyll' is the raster, and 'mcp2015' is a 'SpatialPolygon'

I plotted the output and noticed a problem:

plot(studyarea)
plot(mcp2015, add=TRUE)

The cropped raster is missing cells that are, at least in part, inside the SpatialPolygon. I have used this exact same code, but cropped the raster by a different SpatialPolygon, and I did not have this problem.

I also checked the extent of the cropped raster and the SpatialPolygon, and they are indeed different (ruling out a plotting issue).

>mcp2015

class       : SpatialPolygons
features    : 1
extent      : -130.4485, -129.095, 50.64663, 51.22328  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0

studyarea

class       : RasterLayer
dimensions  : 13, 33, 429  (nrow, ncol, ncell)
resolution  : 0.04166667, 0.04166667  (x, y)
extent      : -130.4583, -129.0833, 50.66666, 51.20833  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
data source : in memory
names       : Chlorophyll.Concentration..OCI.Algorithm
values      : 0.444908, 20.26956  (min, max)

Any ideas?


texts: View in ArcGIS JavaScript (Rest/Services) and API's not working with cache tile services

texts:
In AGS 10.0 I shared the C:/arcgisserver folder, and now my apps and rest endpoints will not display any of my cached services. I can view them fine in ArcCatalog however.  I can build new caches but they are only visible in ArcCatalog.

I am having a similar problem to this post: http://forums.arcgis.com/threads/51216-View-in-ArcGIS-JavaScript-(Rest-Services)-not-working-with-cache-tile-services

AGS Site seems to be battling Windows Server for permissions. Here's how my permissions are set up.

Name:                    Permission Level
Admin                    Owner
AGS SOC                  Read/Write
AGS SOM                  Read/Write
User                     Read/Write


texts: How to programmatically partition a large Polygon in smaller parts within a given Area Threshold?

texts:
I'm looking for a Java library, or even an algorithm to do the following:

I have a large polygon, about the size of a large city. I need to break it down into smaller polygons such that:

The area of each part is less than 1 sq KM. (This threshold should be a parameter to the algorithm or function)
Each Part should have 3 or more sides.
The Part need not be a regular polygon.
Each Polygon can have different area.
These Polygons should cover the entire source Polygon.
The solution should minimize the number of polygons.

How do I Partition the Polygon in such a way?


texts: Domains not connecting to create drop down list?

texts: I am trying to create a file geodatabase with features that have drop down list. I have tried going about this by creating a file geodatabase then adding in feature classes by going to new> Feature Class then right clicking the geodatabase hitting properties going to the domain tab and then adding in the code/description. But when I go into the table it does not have the domain listed in properties. I have made sure the domain is labeled the same as the column. I am running ArcGIS 10.6.


texts: The correct definition of ESRI geodatabase

texts:
Iet's see if i can explain my question correctly.

As the title says I'm looking for an definition of the databasetype that you create in ArcCatalog, file geodatabase and personal geodatabase.

The reason for this is that I'm looking for a list which types of databases that you can work directly with in ArcGIS Desktop Basic without extensions.


texts: Only Interpolating within a polygon feature in ArcMap 10.1

texts: I have sample points within a single polygon feature class that I would like to interpolate, likely using kriging.  The problem is the sample points and polygon boundaries are many miles apart and therefore the kriging interpolates across these many miles.  Is there anyway I could force the kriging to only interpolate within the polygons and not across what looks to be extents of the polygons.  Another way of saying this is let's assume I have sample data points in Missouri and Ohio and I want to compare across these two states so any chance I could merge the 2 states in single feature class and interpolate them together?  It just creates a big rectangle to the extents when I try kriging this data.


texts:
There's no guide to installing OpenLayers on the website and Google is proving fruitless.

Is there an installation file or is it just a bunch of css and js scripts?

If so, where do you paste them?


texts: Converting label expression from VBScript to Python in ArcGIS Desktop?

texts:
I have a feature class with label expression in VB script. I do not have the skill to convert a VB script to a python script.

Here is the VBScript:-

Function FindLabel ( [DESCRIPTION]  )
  if (left( [name] ,4) = "ABCD") then
   FindLabel =PCase( MID( [DESCRIPTION] , 19, 100))
  else
   FindLabel = PCase( [DESCRIPTION] )
  end if
End Function

Function PCase(strInput)
 Dim iPosition
 Dim iSpace
 Dim strOutput
 iPosition = 1
 Do While InStr(iPosition, strInput, " ", 1) <> 0
 iSpace = InStr(iPosition, strInput, " ", 1)
 strOutput = strOutput & UCase(Mid(strInput, iPosition, 1))
 strOutput = strOutput & LCase(Mid(strInput, iPosition + 1, iSpace - iPosition))
 iPosition = iSpace + 1
 Loop
 strOutput = strOutput & UCase(Mid(strInput, iPosition, 1))
 strOutput = strOutput & LCase(Mid(strInput, iPosition + 1))
 PCase = strOutput
 End Function

I came up with following python code. But it did not work:-

def FindLabel([DESCRIPTION]):
    if [name][:4]=='ABCD':
        FindLabel = PCase( [DESCRIPTION][19:100] )
    else:
        FindLabel = PCase( [DESCRIPTION] )

def PCase(strInput):
    iSpace=strInput.find(" ", beg=0, end=len(strInput))
    iPosition=1
    while iSpace!=0:
        iSpace=strInput.find(" ", beg=0, end=len(strInput))
        strOutput=strOutput+strInput[0:1].upper()
        strOutput=strOutput+strInput[iPosition+1:iSpace-iPosition].lower()
        iPosition=iSpace+1
    strOutput=strOutput+strInput[iPosition:1].upper()
    strOutput=strOutput+strInput[iPosition+1:].lower()
    return strOutput


texts: Is there a native JAI implementation for Ubuntu 14.04 64 bit? I'm trying to boost geoserver performance

texts: Is there a native JAI implementation for Ubuntu 14.04  64 bit? I can't find it anywhere.


texts: Two-value raster to generalized polygons (join pixels close together into one polygon) with short computation time in QGIS or ArcMap

texts:
I want to create application map for agriculture based on NDVI map. I have theoretical way how to do it.

My raster has two values (1 - red, 0 - green), the raster is quite big (100 hectares with 6 mm pixel size).

I need to create polygons which contains all the red pixels (as on the picture) - ideally with some buffer. Result should be shapefile (so Raster to Polygon) with maximum number of 50 000 polygons, so I need to generalize polygons and merge to achieve smaller amount and shorter computation time.

Polygonize often use big amount of time. I also tried make just worse resolution, it works quite well, BUT I need that the new pixel (for example 1 m resolution) need to be red even if only one former pixel was red. Any red place can't disappear. Best would be if I don't have to polygonize it first because of waste of time, but generalize and merge it first and then polygonize.

In QGIS or ArcGIS.


texts: Retrieve the WMS Request URL of a Geoserver layer on button click

texts:
I have some Geoserver layers, which are displayed using leaflet or geoexplorer.
What I need is to add a button, which onclick gets the WMS url for this layer.

Using the developers tools of Chrome, I can see the WMS request. I could try to build the URL by concatenating all the different parameters.

I wonder if there is another more elegant way to accomplish this. I was checking the gsconfig library but the get_layer and getresource functions don't think they return the WMS request URL.

Any ideas?


texts:
I have 2 multiband rasters and I want to mask some areas of one of them with values of the other one, so I write this code (I can't figure out why is not working):

import os, rasterio
import numpy as np

feb = r'path\to\multiband\stack_1.img'
may = r'path\to\multiband\stack_2.img'
mask = r'path\to\mask_raster' #this one is just one band
outfile = r'path\to\outfile' #multiband

with rasterio.open(mask) as mascara:

    MASK = mascara.read()

with rasterio.open(feb) as febrero:

    FEB = febrero.read([1, 2, 3, 4, 5, 6, 7, 8, 9, 10])

with rasterio.open(may) as mayo:

    MAY = mayo.read([1, 2, 3, 4, 5, 6, 7, 8, 9, 10])

for i in range(1,11):

    MAY = np.where(MASK == 1, FEB[i], MAY[i])

profile = mayo.meta
profile.update(dtype=rasterio.int16)

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

I get this error when running the code:

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-21-054d59b9e8d1> in <module>()
     22 for i in range(1,11):
     23
---> 24     MAY = np.where(MASK == 1, FEB[i], MAY[i])
     25
     26 profile = mayo.meta

IndexError: index 2 is out of bounds for axis 0 with size 1

I am able to do it one by one band, but then I have to do the stack again and overall I would like to know why the code doesn't works.


texts: How to do a Point to Polygon analysis with multiple variables in QGIS?

texts:
I have a database of almost 10k crimes, which are divided in to 10 categories and plotted on a map

I also have a database of 12k addresses, which are also plotted on the map.

I aim to perform a multiple regression to determine whether there are correlations between types of crime and house price.

I have created a buffer of 250 meters around each of the addresses, and would like to list a count of each of the crimes which are located within this buffered zone for each of the addresses, however, I can only get QGIS to give me a count of all types of crime, rather than the specifics.

Does anyone know how to do a point to polygon analysis of this type? I am in the middle of breaking the crime databases to analyse each one separately, which is a bit difficult!


texts:
I am using Google Satellite imagery as a base map to capture the vector data in QGIS 3.2.2 for past couple of months. Yesterday i noticed that the google imagery has offset issue in some locations where i have already worked on. Since that location was already captured with reference to the google imagery but now it shows major offset(approx 70 meters) between the vector and google imagery. This location is belongs to California, USA. (-119.17907,34.23393).

I have attached the screen shot of the vector location which was captured with reference to the Google imagery(before imagery updates) and this location almost matching with bing imagery but it showing major offset while using recent Google imagery as a base map. I hope this explains well about this issue.

I also attached screen shot from the Google Earth Pro and its almost lined up with the vector data as similar to bing imagery.

Does google updated their imagery on these locations, if yes then how its possible they update imagery with major offset issue?
Does anyone can let me know what is the issue for this and how it will be resolved. Please help me on this.


texts: Delineate polygonal catchment areas from a raster using open-source GIS

texts:
I'm 90% sure this basic question has been asked before but looking through the other posts I simply do not see a useful answer.

I need to calculate polygonal(vector) catchment areas from a raster based on minimum driving distance/cost from a set of points. I've already got the raster, and for that matter, the nodes(intersections with cost attribute) used to interpolate it. What I need now are polygons for each point describing the area nearer it than any other point. This is basically like a voronoi diagram but based on network distance rather than euclidean space.

How exactly can I do this using only open-source tools? I'm comfortable with the command line if I just know what I'm looking for. I feel like GDAL is going to be the answer, but if so, what exactly is the command I'm looking for?

Thanks!


texts: Convert coordinates from degree unit to meters unit

texts: I have a CSV file including some coordinates in degrees. I need to have these coordinates in the meter unit. I tried to add it as a shapefile in QGIS, then I saved the features and set a projected CRS. Afterward, I read this new shapefile into QGIS but still, the longitude and latitude values are in degrees. How can I convert the longitude and latitude values to meters?


texts: Nearest point to another point with condition in QGIS 3

texts:
I have two different statewide layers of points. One set of points are locations of ballot dropboxes and the other points are residences. Both layers have a column for city. I'm trying to calculate the nearest dropbox to each residence and end up with a table of the closest dropbox to the residence that retains the residence unique id and all the attributes of the dropbox layer (geo, address string, note string). The catch is they must be in the same city. For example, in northern Boston, MA, the nearest dropbox point by distance is in southern Cambridge, MA but Boston voters cannot drop off their ballot in Cambridge; it must be a Boston dropbox if the residence is in Boston.
I'm trying to understand QGIS better and it seems like this is maybe something I could iterate over although I'm not sure how I'd go about something like residence.city == dropbox.city for each city. Or is it more of a batched process or something that would require a model?


texts: I have a network link KML which I'm downloading via Python script and then converting to shapefile in ArcGIS Desktop.  The problem I'm having is the max field length of 254.  It's truncating the PopupInfo field in my KML.  I need to convert all the data (longer than 254 characters!) to a feature layer or shapefile, so I'm seeking alternatives.  I have not been able to find a good answer here.


texts: QGIS Atlas export produces map in smaller scale then in composer

texts:
I've loaded in QGIS 2.1, OSM map with OpenLayers plugin and I've created two rectangles representing coverage layer for Atlas. Then I've set Atlas in composer, first picture below. This preview should correspond to final page of the exported map. However, when I export Atlas in any format the scale of the picture is much smaller, see the second picture.

Did I miss anything or is this normal behaviour?


texts: Dealing with custom projections (from prj file) in PostGIS (also GeoServer and OpenLayers)

texts:
Basically I want to perform spatial queries like distance calculations in PostGIS using data from shapefiles with different projections, some of which have no SRID code that I can find and also which the python osr package can't seem to handle. What do I do? I have little understanding of projections and how to deal with them.

Full story:

I'm creating a system using OpenLayers, GeoServer, and PostGIS. My data is from shapefiles that I loaded into PostGIS using shp2psql. Most of the shapefiles didn't have a .prj file associated with them and seem to just be long/lat / EPSG:4326. However for some of the files the coordinates have numbers in the millions which are clearly not degrees (fortunately there were .prj files here). In PostGIS the data appears as a geometry column, not geography. Until now I've been mostly concerned with displaying the data in OpenLayers, so I took care of projection in GeoServer by specifying an EPSG code when I added a new layer. Usually this was 4326, once it was a standard code I found using http://www.prj2epsg.org/search, and once I created a custom projection by adding the WKT to the epsg.properties files in the data directory. Everything seems to display fine in OpenLayers then; it's hard to be sure since I don't know the data personally but it looks like it all lines up. I haven't had to write any code on the client side (OpenLayers) relating to projections, everything is in terms of longitude and latitude and very easy.

Now, suppose I specify in the client (e.g. by clicking on the OpenLayers map) the coordinates of a point using longitude and latitude, and a distance such as 100 km, and I want to display all the features from a certain layer that are within said distance from said point. So, in PostGIS I'd like to say something like:

SELECT gid, geom
FROM layer
WHERE ST_Distance(geom, ST_GeomFromText('POINT(24 -30)')) < 100000;

Right now this returns everything (for most layers) since the geom column is a geometry type instead of geography. As far as I can see, I need to turn everything into a geography type, transforming into the EPSG:4326 where necessary. Then I should just be able to declare everything in GeoServer as EPSG:4326 and no longer worry about anything there, and ST_Distance should return results in metres calculated on the spheroid.

But for some of my data I can't see how to do this. For example, here is the text from the .prj file for one of the layers:

PROJCS["SA Albers Equal Area",GEOGCS["GCS_Hartebeesthoek_1994",DATUM["D_Hartebeesthoek_1994",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Albers"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",24.0],PARAMETER["Standard_Parallel_1",-24.0],PARAMETER["Standard_Parallel_2",-33.0],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]

This is a projection for South Africa, and none of the results from http://www.prj2epsg.org/search look relevant. So I want to define it as a custom projection in PostGIS. As far as I understand I can't use WKT for this and need the parameters to be expressed in proj4. I tried running it through this python script seen around here:

#!/usr/bin/env python

import osr
import sys

def main(prj_file):
    prj_text = open(prj_file, 'r').read()
    srs = osr.SpatialReference()
    if srs.ImportFromWkt(prj_text):
        raise ValueError("Error importing PRJ information from: %s" % prj_file)
    print srs.ExportToProj4()
    #print srs.ExportToWkt()

if __name__=="__main__":
    main(sys.argv[1])

and I got the following:

ERROR 6: No translation for Albers to PROJ.4 format is known.

Now I don't know how to proceed. I'd also like confirmation that my reasoning about everything else I've mentioned is sensible.


texts: How do I manage to get my swedish letters

texts:
I am trying to do the map in ´ggplot2´ in ´R´. I use this function from rgdal package:

mapa <- readOGR(dsn=path.expand("C:/~/SWE_adm_shp"),layer="SWE_adm2")

than I do:

fortify(mapa)

ggplot(mapa, aes(x = long, y = lat, group = group)) +
  geom_polygon(fill = "blue", col = "white") +
  coord_map()

It works fine but when I look at the

head(mapa)

I see that the names of the counties dont have the Swedish font (NAME_2)

ID_0 ISO NAME_0 ID_1         NAME_1 ID_2      NAME_2   HASC_2 CCN_2 CCA_2   TYPE_2    ENGTYPE_2
0  222 SWE Sweden    1 Ã–stergÃ¶tland    1   Ã–deshÃ¶g SE.OG.OD     0  <NA> Kommuner Municipality
1  222 SWE Sweden    1 Ã–stergÃ¶tland    2 Ã…tvidaberg SE.OG.AT     0  <NA> Kommuner Municipality
2  222 SWE Sweden    1 Ã–stergÃ¶tland    3     Boxholm SE.OG.BO     0  <NA> Kommuner Municipality
3  222 SWE Sweden    1 Ã–stergÃ¶tland    4   FinspÃ¥ng SE.OG.FI     0  <NA> Kommuner Municipality
4  222 SWE Sweden    1 Ã–stergÃ¶tland    5       Kinda SE.OG.KI     0  <NA> Kommuner Municipality
5  222 SWE Sweden    1 Ã–stergÃ¶tland    6  LinkÃ¶ping SE.OG.LI     0  <NA> Kommuner Municipality
  NL_NAME_2 VARNAME_2
0      <NA>      <NA>
1      <NA>      <NA>
2      <NA>      <NA>
3      <NA>      <NA>
4      <NA>      <NA>

How to work over that?


texts: How to create independent/fixed inset maps with leaflet?

texts: I am trying to create a leaflet map of the USA focusing on the contiguous 48 states and then add two inset maps of Alaska and Hawaii on the side. The problem is that addMiniMap() option creates an interactive inset mini map that just shows the location of the zoom window of the main map. Rather than this, I want to have inset maps with fixed position and zoom on Alaska and Hawaii, while the main map remains interactive.


texts: Failed to open raster - ArcGIS 10 SP2, Python 2.6.5, Windows 7

texts:
I have an arcpy / python script that is working fine with ArcMap 10 no Service pack (Build 2414) on windows XP, but consistently failing on using ArcMap 10 SP2 (Build 3200). Getting the ESRI 99999 error. In a nutshell, this script is supposed to reclassify and convert four raster files to polygon (shp).

Here is the code. Any help or ideas would be appreciated.

import arcpy, os, string
arcpy.CheckOutExtension("Spatial")
from arcpy.sa import *

#set output location
outLoc = "O:/gis/!KEYHOLE/DAILY_SHP"

#get file date, set data directory
fileDate = arcpy.GetParameterAsText(0)
dataFolder = "O:/gis/PRED_SVC/MET_PROJ/FOR_OBS_LTG_DATA"

#add file date message
arcpy.AddMessage("...using file date " + fileDate + "...")

#set precip amount file based on date input
prefixList = ["erc", "thz", "prpdur_", "prpamt_"]
for prefix in prefixList:
metName = prefix + fileDate
metFile = os.path.join(dataFolder, metName)
print metName
print metFile

#output workspace, turn overwrite output on
arcpy.env.workspace = "c:/tmp"
arcpy.env.overwriteOutput = 1

#reclassify
if prefix == "erc":
    remap = RemapRange([[-99, -1.282, 1], [-1.282, -0.842, 2], [-0.842, -0.524, 3], [-0.524, -0.253, 4], \
             [-0.253, 0, 5], [0, 0.253, 6], [0.253, 0.524, 7], [0.524, 0.842, 8], [0.842, 1.282, 9], \
             [1.282, 1.645, 10], [1.645, 1.881, 11], [1.881, 99, 12]])
    arcpy.AddMessage("...reclassifying ERC...")

elif prefix == "thz":
    remap = RemapRange([[-99, -1.645, 1], [-1.645, -1.282, 2], [-1.282, -0.842, 3], [-0.842, -0.524, 4], \
             [-0.524, -0.253, 5], [-0.253, 0.253, 6], [0.253, 0.524, 7], [0.524, 0.842, 8], \
             [0.842, 1.282, 9], [1.282, 99, 10]])
    arcpy.AddMessage("...reclassifying 1000HR fuels...")

elif prefix == "prpdur_":
    remap = RemapRange([[0, 0.009, 1], [0.009, 1.9999, 2], [1.9999, 4.9999, 3], \
                        [4.9999, 9.9999, 4], [9.9999, 25, 5]])
    arcpy.AddMessage("...reclassifying precipitation duration...")

else:
    remap = RemapRange([[0, 0.0009, 1], [0.0009, 0.09999, 2], [0.09999, 0.249999, 3], \
                        [0.249999, 0.499999, 4], [0.499999, 0.999999, 5], \
                        [0.999999, 1.999999, 6], [1.999999, 1000000, 7]])
    arcpy.AddMessage("...reclassifying precipitation amount...")

rcls = Reclassify(metFile, "Value", remap)

#make raster layer, convert to poly
metLyr = arcpy.MakeRasterLayer_management(rcls, "metLyr")
arcpy.env.workspace = outLoc
arcpy.env.overwriteOutput = 1
outPolyName = string.replace(prefix, "_", "") + "_poly.shp"
outPolyFile = outLoc + os.sep + outPolyName
arcpy.AddMessage("...converting " + metFile + " to poly...")
print outLoc
print outPolyName
print outPolyFile
arcpy.RasterToPolygon_conversion(metLyr, outPolyFile)
del rcls

arcpy.CheckInExtension("Spatial")


texts: No visable Difference between two Lists for Visualisation via Palette creating an Error

texts:
I want to Visualize a Landcover Classification

I have different success with seemingly identical Lists.

List 1:
Pallette Parameters inside the Code  like:

var landcoverPalette1 = [
'E6004D',   //Artificial surfaces > Urban fabric > Continuous urban fabric
'FF0000'    //Artificial surfaces > Urban fabric > Discontinuous urban fabric
]

List2:
Palette Obtained via Code

var landcoverPalette = ee.List(lc2012.get('landcover_class_palette'))

List1 and List2 are identical withing the print() statement(see Code), but List2 generates an Error

2012: Layer error: Image.visualize: Expected a string or list of strings for field 'palette'.

Code without Error

Code with Error


texts:
I saved some features from a WFS (http://geodienste.bfn.de/wfs/schutzgebiete - the theme "Landschaftsschutzgebiete"). We did this in the past without any problems.

Now, despite the features having attributes, I can't select a feature in the saved/downloaded data, nor use the identify-tool on them. Furthermore, when zooming in, the features disappear - unevenly. Some are visible only up to 1:1000000, others up to 1:500000, etc.

Has anyone encountered this behaviour, and how did you solve it?


texts: I have a postgis layer (layer.sql), which I am trying to add to my postgres database (mydatabase). The postgis layer is stored on S3 (https://s3.us.../layer.sql). I can add the layer to my database when it is stored on my local machine using psql mydatabase < ~/layer.sql but I'm failing to add the layer to my database, if I'm using the url. I tried psql mydatabase < https://s3.us.../layer.sql, which returns -bash: https://s3.us.../layer.sql: No such file or directory. Is there a different way command I can use when trying to add the layer with a url?


texts: Get list of equidistant points on polygon border PostGIS

texts:
I would like to get points from polygon border using PostGIS.

I've tried using ST_DumpPoints() function, like this:

CREATE OR REPLACE FUNCTION get_points_from_polygon(polygon geometry)
RETURNS SETOF geometry AS
$$
DECLARE
    point geometry;
BEGIN
    FOR point IN
        SELECT DISTINCT
            points.geom
        FROM (
               SELECT (ST_DumpPoints(polygon)).*
             ) AS points LOOP
        RETURN NEXT point;
    END LOOP;
END;
$$
LANGUAGE plpgsql ;

But it gets sometimes huge list of points (~100-500) and that is not what I need.
I need to get N equidistant points on the border, like on the screenshot attached.
Then, when I have these points I will analyze viewshed polygons from point as observer

I've found, how to identify is the point on the polygon border
But have no idea , how to create first them and get as a list of points geometry.

What would be your suggestion to do this?


texts: Is it possible to use .elemgeo dot macro on the pyspatialite?

texts:
I need to perform an explode, like PostGIS ST_DUMP on the QGIS plugin using pyspatialite.

There is a dot macro .elemgeo but I don't know if it is possible to call it from the pyqgis plugin via pyspatialite.

The version I use doesn´t have ElementaryGeometries table.

QGIS 2.8.1 Wien
Sqlite 3.7.17
SpatiaLite 4.1.1


texts: How can I change the coordinate system of gpx-points?

texts:
I'm working with qgis, v. 2.0.1. I have some points from a gps device (Garmin) and they already are in the .gpx format. Now I have two problems: First, I can't open the gps-Points with the "gpx-tools" plugin by simply choosing "load gpx file". I always have to use the gpx-converter, otherwise the points don't appear in the map window. But actually it should work with "load gpx file", shouldn't it?
The second problem is that I want to change the coordinate reference system of the points from WGS84-EPSG:4326 to CH1903/LV03-EPSG:21781. But as soon as I do this, the points disappear...What can I do?
Thanks.


texts: How to get my store locations which are withing 30 min drive distance using google map

texts:
I have the list of stores in my database.
I want to show only those stores to a customer who can reach that store within 30 minutes (configurable) of drive.
using google api


texts: Page becomes slow after adding some points at an interval

texts:
I have an interval that sends request to the server every 2 seconds and updates the map data.

The server sends json and in Javascript i iterate over and set some styles and features. However i noticed that the page becomes very slow after a while. Probably something has to be reseted, each time i call the posgres function.

See it here

   function Posgres(data){
    function iconStyle(iconSrc) {
        return new ol.style.Style({
            image: new ol.style.Icon({
                src: iconSrc
            })
        })
    }
    data.DEVICES.forEach(feature => {
        let popupData = p(feature)
        let iconName = feature.ICON
        let iconFeature = new ol.Feature({
            geometry: new ol.geom.Point([Number(feature.longitude), Number(feature.latitude)]).transform('EPSG:4326', 'EPSG:3857'),
            name: `<b>Name:</b>` + feature.name + `<br>` + popupData + `<br><b>Received at:</b>` + feature.received_at
        })
        iconFeature.setStyle(iconStyle(iconName))
        vectorSource.addFeature(iconFeature)
    })
}
function sendRequest() {
    $.post("MapMarkers.php", null, null, "json")
        .done(data => {
            Posgres(data)
        }).fail(function() {
            alert("error loading data")
        })

       var element = document.getElementById('popup')

    var popup = new ol.Overlay({
        element: element,
        positioning: 'bottom-center',
        stopEvent: false,
        offset: [0, -15]
    })
    map.addOverlay(popup)
}

sendRequest()

setInterval(() => {
    sendRequest()
}, 2000)


texts: Accessing Home directory of ArcMap using ArcPy?

texts:
I want to access to a directory and set that as a arcpy.env.workspace in my Python Script Tool. However, this workspace depends on users so I cannot just use my own folders.

Is there a way to access Home directory (where .mxd file is saved) in arcpy?

I tried retrieving the directory from os.getcwd but that refers to "c:\windows\system32" where ArcMap is accessing.


texts:
I have a layer with coord. ref. : NA.

How can I assign coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 to coord. ref. : NA?

class       : RasterLayer
dimensions  : 127, 146, 18542  (nrow, ncol, ncell)
resolution  : 0.01, 0.01  (x, y)
extent      : -122.39, -120.93, 36.72, 37.99  (xmin, xmax, ymin, ymax)
coord. ref. : NA
data source : C:\2012Tif\2012ASC5min_ppt_spas1306_0001_20121222_0005_UTC.asc.tif
names       : X2012ASC5min_ppt_spas1306_0001_20121222_0005_UTC.asc
values      : 0, 0.0782172  (min, max)

class       : RasterLayer
dimensions  : 180, 240, 43200  (nrow, ncol, ncell)
resolution  : 0.01, 0.01  (x, y)
extent      : -123.305, -120.905, 36.705, 38.505  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +to wgs84=0,0,0
data source : C:\Users\nahmlee\Documents\1998Test.tif
names       : X1998Test
values      : 0, 0.07884484  (min, max)


texts: Google Earth Engine: calculate percentile rank in ee.ImageCollection

texts:
I am trying to calculate the Vegetation Productivity Index sensu copernicus in Google Earth Engine, which requires estimating the percentile rank. I haven't found any lead so far on how to do so, does anyone know how to do it or recommend a way to explore?
Given the example below, how would you calculate the percentile rank of the NDVI values in img compared to those in ref? What I am expecting is a value from 0-1, where the value at each pixel is the rank compared to all historical values for the same pixel.
ref = ee.ImageCollection('MODIS/006/MYD13Q1').filterDate('2006-01-01', '2006-12-31').select('NDVI')
img = ee.ImageCollection('MODIS/006/MYD13Q1').filterDate('2007-01-01', '2007-01-10').select('NDVI').first()


texts: Linking Google Maps API for business from OpenLayers?

texts:
I just came across a web-page about supplying private key from OpenLayers to Google Maps for business API. But checking on the Google Maps API page, it is asking for signatures which is hash of the query with a secret key for each request. I think it is for maps image requests also, which is contradicting with what the web-page says.

Can the Google Maps for business API be used with OpenLayers and, if so, then how?


texts: R raster package problem with extract() function

texts:
I'm trying to extract pixel value using the extract() function from raster package in R and get the following error, bellow in my code bellow a well, where: df$PixelNumber is the "cell number value". The error :

  Error in (function (classes, fdef, mtable)  :    unable to find an
  inherited method for function ‘extract’ for signature ‘"character",
  "numeric"’

the code I'm using:

vec=as.vector(as.numeric(as.matrix(df$PixelNumber)));
rast = file.choose();
VlaueSlope=extract(x= rast, vec)

What could I be doing wrong in the extract() code that cannot get the pixel value in the variable VlaueSlope.


texts: Can't Find the Split selected features with selected line from another layer tool

texts:
I'm trying to split a selected polygon features with a selected line from another feature layer.

I, however, can't seem to get the tool for this to show up in my advanced digitize toolbar.

In other words, I need someone to tell me how to get my advanced digitizing toolbar to look like this:

  My current toolbar only contains the first half of these tools.
I need to gain access to the
  "Split selected feature with select line from another layer" tool and have it  show up at the top of my screen.


texts:
I have a layer with polygons in Openlayers 4 that I want to show as polygons when zoomed in and cluster them when I zoom out to a given level.
Ol.source.Cluster has a geometryFunction attribute where I translate the polygons into points which is necessary for the Cluster object (this seems reasonable to me, the distance between polygons are not intuitive)
The problem comes when I apply a Polygon Style to the layer, since the Cluster function transforms the polygons into points it doesn't show anything, is there a way to get around this using only one layer?
var zoomedInStyle = new ol.style.Style({
    image: new ol.style.Circle({
        radius: 10,
        fill: new ol.style.Fill({
            color: 'orange'
        })
    }),
    geometry: function(feature) {
        // Here I would need the original polygon
        return feature.getGeometry().getInteriorPoint();
    }
});

var zoomedOutStyle = new ol.style.Style({
    stroke: new ol.style.Stroke({
        color: 'orange',
        width: 3
    }),
    fill: new ol.style.Fill({
        color: 'rgba(255, 165, 0, 0.1)'
    }),
    geometry: function(feature) {
        var type = feature.getGeometry().getType();
        if (type === 'Polygon')
            return feature.getGeometry().getInteriorPoint();
    }
});

function featureStyleFunction(feature, resolution) {
    if (resolution > 0.5)
        return zoomedOutStyle;
    else
        return zoomedInStyle;
};

var source = new ol.source.Vector({
    features: (new ol.format.GeoJSON()).readFeatures(geoJson)
});

function clusterFunction(feature) {
    return feature.getGeometry().getInteriorPoint();
};

var clusterSource = new ol.source.Cluster({
    source: source,
    geometryFunction: clusterFunction,
    distance: 20
});

var layer = new ol.layer.Vector({
    source: clusterSource,
    style: featureStyleFunction
});


texts: R : How to build heatmap with the leaflet package

texts:
I read a post about interactive maps with R using the leaflet package.

In this article, the author createD a heat map like this:

X=cbind(lng,lat)
kde2d <- bkde2D(X, bandwidth=c(bw.ucv(X[,1]),bw.ucv(X[,2])))

x=kde2d$x1
y=kde2d$x2
z=kde2d$fhat
CL=contourLines(x , y , z)

m = leaflet() %>% addTiles()
m %>% addPolygons(CL[[5]]$x,CL[[5]]$y,fillColor = "red", stroke = FALSE)

I am not familiar with the bkde2D function, so I'm wondering if this code could be generalized to any shapefiles?

What if each node has a specific weight, that we would like to represent on the heat map?

Are there other ways to create a heat map with leaflet map in R ?


texts: City Engine - 3D zoning envelope & street orientation

texts:
Is there a way to change the front & back angles that generate the 3D zoning envelope according to street orientation? Ideally there should be 6 cases which describe angular ranges and for each case a different angle is assigned to front,back and sides of the lot.

Would this require scripting? I am a complete novice in this and any help would be appreciated.


texts: Add Photo Image to ArcGIS Online, need to format the url

texts: I would like to delete the last 10 letters of a URL to use in ArcGIS online to show photos in the popup. Any suggestions? I have seen people put JavaScript in the image configuration settings to get the photos to show in  Show attached images inside popup with ArcGIS Online?. Maybe deleting part of a URL is possible too.


texts: Getting value from another feature depending on value of field using QGIS Field Calculator

texts:
I'm a total newbie at coding.
I'm working with water catchments in QGIS. Runoff of some catchments feeds other catchments. I'm trying to get the cumulative run-off.

Each feature in the table has the following fields:

ID - each polygon's number
Fed by - the numerical ID of a polygon that feeds it
C - the (calculated) runoff capacity of that field

I would like to create a field C-in that pulls the value of C from another feature based on the ID of that feature if it is written in Fed by field.

I imagine it is something along the lines of: When Fed by equals the value of an ID, get C of that ID, for example:

In particular, we have:

ID: 19, Fed by: 9, C: 0.5
ID: 20, Fed by: 19, C: 0.7

So C-in of ID=20 should be: 0.5 (forget for a moment the C-in of ID=19)

So what do I need to write in the field calculator to get this magic to happen?


texts:
I have a point file of trees in my study area, which includes the attributes of each tree.   Each tree also has a crown radius associated with it.    I need to summarize the attributes of these points into a 5-meter raster file.    The 5-meter raster file placement is defined by a 5-meter fishnet.

I would like a suggestion as to which method is the most statistically robust.

One method I thought might work is:

For a given 5 X 5 meter cell,  take the area of each crown that occupies a cell and multiply it by the attribute of interest.   Do this for each crown and sum these values and divide by the total area of crown occupying the given cell.

i.e.
Let the shape i be have the value of some attribute of the Tree … say Crown bulk density,   (CBDi )    Then the value of the grid cell X = sum(CBDi * Ai) / sum (Ai)   where Ai is the area of each crown within the cell.

I don't know how to perform this method in ArcGIS though.

Another simple method:

Do a point to raster and calculate the mean value and snap the raster to the defined raster.
Choose either bilinear or nearest neighbor for the resampling method in the environment settings.


texts: Insert sequential row and column values automatically into a fishnet dataset

texts:
I have a grid that has about 30 rows with approx 15 columns each.

I want to label them as A1,A2,B1,B2 etc (Alphabet row, integer is column). Is there any way to do this programatically?

There is no fixed number of columns in a row, so we can not say that after every x columns the next item is B etc.

If necessary I guess we can compare the y on the top right of each grid and when it changes Row A becomes B but I am hoping that there is an easier way of doing this in SQL (Mapinfo) or QGIS.

UNFORTUNATELY I don't have access to ArcGIS grid tool that does this.

I used "QGIS--Vector--Research--Vector Grid" to make the grid.

Another way of asking this is, can this QGIS tool be modified to include grid numbering like in the ArcGIS equivalent Grid Index Features (Cartography) and also ideally the surrounding map numbers.

I am also happy to use a mapbasic program to do this...


texts: 'pyrcc5' is not recognized as an internal or external command

texts:
I installed QGIS 3.0 on Windows 10. I get this error everytime i run the pyrcc4 or pyrcc5 command:

'pyrcc4' is not recognized as an internal or external command

Even if I find a place with a pyrcc4.exe file or a place with apyrcc5.bat file and add the path to the environment variables I get the same error.
In the previous installation of osgeo4w i had a pyqtdev folder or something similar appear I think it was on C/. Was there supposed to be such a folder, do I need to download pyqt separately?


texts:
I was given an unprojected coordinate in the form of lat, lon and I am using pyproj to convert to projected coordinate as such:

from pyproj import Proj

# New Jersey
lat = 40.704895
lon = -74.094994

wgs84 = Proj(proj='utm', zone=18, ellps='clrk66')
x, y = wgs84(sample_data['lon'], sample_data['lat'])

# x = 576453.375, y = 4506180.871
# 576453 meters east and 4508180 meters north of the zone's origin

I'd like to know the origin (0, 0) of the UTM zone 18N. I tried using inverse conversion with

wgs84(0, 0, inverse=True)

I got (-79.48869479772156, 0.0) which means it's 79 meters south and right at the Prime Meridian which seems wrong.


texts: QGIS SVG-Fill pdf-output converted to EMF via Inkscape black borders

texts:
I imported a "wave" SVG fill for lakes. Exporting a PDF is nice. But when I convert the PDF to EMF (for Word import) with Inkscape the lakes have a black background around:

How can I get rid of this?


texts: ArcGIS Annotation Font Symbol Collection

texts:
Is there a way in Arcpy, to loop through the symbol collection to get back some attributes of each item in the symbol collection list (symbol ID, symbol name, font size, etc.)?

I need to do this in arcpy, but what I'm talking about is the information that you see in ArcCatalog's Annotation tab of the Feature Class Properties dialog. There is a list box with the title "Symbol Collection". The symbol ID appears in parenthesis and right next to it is symbol name.


texts: Point data location to be set by Polygon ID

texts:
This might be stupid but I am trying to put points collected on the field to match the location from where they were originally collected. AS the X and Y coordinate of the point might not match the location of the polygon. Let me explain:-

For example If I have a point with the Same ID (imported from a CSV) I want it to show up on top of this polygon. This point will have some additional information relevant to this polygon and is updated every month. I know I can use the X & Y coordinates of the point but the issue is that the X & Y coordinates of point are not reliable as they change every time the point is collected.

I am using QGIS.


texts: Converting layer file ArcGIS desktop 10.2 to 10.0 fails

texts:
I am having problems with converting an ArcGIS 10.2 layer file to a 10.0 layer file.
I used ArcGIS online to import webmaptiles into my ArcMap 10.2.2. desktop version. In order to use them offline, I saved them as a layer file.

My colleague uses ArcMap 10.0, so I saved the layer file as a 10.0 version.
However, when adding this layer in ArcMap 10.0 he gets following error message:

Adding the selected layer file failed.
Possible reasons include that it may have been created using a newer version of ArcGIS and can't be opened in the version you are using.

What causes this conversion to fail?

I tried importing the webmaptiles using ArcGIS online on his computer, but this did not work either.


texts:
Is there any way in SQL to supress errors in PostGIS functions, for instance if some of your geometries are invalid but you want to report a null value for the output or something?

To give an example, I have a table with a structure:

id, geom1, geom2

And I would like to run a query like select id, st_split(geom1, geom2) as geom from mytable, but whenever there is a problem with the geometry, I get ERROR:  Splitter line has linear intersection with input and the process stops. What if I want to discard this geometry and carrying on, giving a null result for the geometry?

I presume you can do this in PL-pgSQL, and I will post a solution here if I figure that out, but is there no way of doing it in plain SQL?


texts: Loading Sentinel 2 images into ArcGIS Mosaic Dataset?

texts:
I'm attempting to add Sentinel 2 images to a new mosaic dataset in ArcGIS 10.4.

AddRastersToMosaicDataset reports success, but there Errors reported in the messages and the result is that the images are not added to the mosaic dataset.

I've ensured that the path length for the files are <256 (https://geonet.esri.com/thread/173068):

  Please look to shorten the path by truncating part of the directory
  name that contains .SAFE at the end. All the details in the directory
  name are part of the metadata, so such a long directoy name is not
  required.

It appears to successfully find, but is unable to open the raster and removes it.

Messages
Executing: AddRastersToMosaicDataset C:\ssdvol1\GISData\s2test\sentinel_2016.gdb\sentinel_2016_md Sentinel-2 M:\S2A_MD_TEST\S2A_OPER_PRD_MSIL1C_PDMC_20160206T011622_R098_V20160204T175251_20160204T175251.SAFE\S2A_OPER_MTD_SAFL1C_PDMC_20160206T011622_R098_V20160204T175251_20160204T175251.xml UPDATE_CELL_SIZES UPDATE_BOUNDARY NO_OVERVIEWS # 0 1500 # *S2A*.xml SUBFOLDERS ALLOW_DUPLICATES NO_PYRAMIDS NO_STATISTICS NO_THUMBNAILS # NO_FORCE_SPATIAL_REFERENCE NO_STATISTICS #
Start Time: Tue Mar 22 13:41:23 2016
2016-03-22T13:41:25.799: Loading raster datasets
2016-03-22T13:41:25.829: Completed crawling 1 data source items. Added 1 mosaic dataset items.
2016-03-22T13:41:25.830: Synchronizing crawled data source items
2016-03-22T13:41:25.844: Synchronizing items associated with raster type instance 'Sentinel-2' [ID: 10].
2016-03-22T13:41:25.862: Error: 80042019: Could not build mosaic dataset item. [ID: 3, URI: 'M:\S2A_MD_TEST\S2A_OPER_PRD_MSIL1C_PDMC_20160206T011622_R098_V20160204T175251_20160204T175251.SAFE\S2A_OPER_MTD_SAFL1C_PDMC_20160206T011622_R098_V20160204T175251_20160204T175251.xml|MS']
2016-03-22T13:41:25.862: Error: 8004203f: Could not open source as a raster dataset. [M:\S2A_MD_TEST\S2A_OPER_PRD_MSIL1C_PDMC_20160206T011622_R098_V20160204T175251_20160204T175251.SAFE\Granule\S2A_OPER_MSI_L1C_TL_MTI__20160204T212051_A003241_T13SDD_N02.01\IMG_DATA\S2A_OPER_MSI_L1C_TL_MTI__20160204T212051_A003241_T13SDD_B01.jp2]
2016-03-22T13:41:25.865: Error: 8004205e: 1 mosaic dataset items could not be built.
2016-03-22T13:41:25.865: Removing invalid mosaic dataset items
2016-03-22T13:41:25.876: Removed 1 mosaic dataset items and 0 unreferenced raster type instances.
2016-03-22T13:41:25.948: Removing broken mosaic dataset items.
2016-03-22T13:41:25.958: Removed 0 mosaic dataset items and 0 unreferenced raster type instances.
2016-03-22T13:41:26.032: Completed synchronization: 1 items selected, 0 items synchronized.
Succeeded at Tue Mar 22 13:41:26 2016 (Elapsed Time: 2.66 seconds)

Is there something I need to do to the sentinel directory or files in order for ArcGIS to add it to mosaic dataset?


texts: Does using RAM Disk improve ArcGIS Desktop performance appreciably?

texts:
I configured a RAM disk on a Virtual PC running ArcGIS Desktop 10 and set my  HKCU\Software\ESRI\Output\TempPath to a directory on the RAM disk.  It’s hard to perceive any sort improvement with that sort of configuration and I’m curious to know if anyone else has done this on a non-virtual computer and what the results were.

Anecdotal answers are fine.

Also, putting the cache path on a ram disk is probably not a good idea, but is anyone doing this anyway?  (HKCU\Software\ESRI\Desktop10.0\Settings\Cache Path)


texts: Generate an accurate DTM with photogrammetry?

texts:
I generated imagery from a low flying drone (40 M AGL) equipped with a mica sense RedEdge multi-spectral camera and have processed the imagery in PIX4D and now am attempting to do some analysis in QGIS.

I am working to create CHM for several research plots and am noticing that the minimum height in DTMs generated from imagery on the same plot differ by about up to 2 meters. This is a big difference and will drastically skew results. Some of the plant heights are below 2 meters.

Any idea on how to generate a highly accurate DTM that can be used in CHM analysis?

I do have imagery from winter months and have considered generating a DTM from the leafless season but wonder about differences in ground level due to frost heaving and changes due to erosion. Any one have experience with canopy height modelling?


texts:
I had already formulated an entry recently in which I asked for the error message "msCGILoadMap (): Web application error. CGI variable" map "is not set." asked.

Now I think I have found a more accurate diagnosis. For our project we need several mapserv.exe, because we present different maps to different users (Mapfiles).
Every time I come up with a new copy of the mapserv.exe this error message follows. Does anyone know this problem? Does something have to be considered when copying the file?


texts: I'm using the "esri.virtualearth.VETiledLayer" in my map application as a basemap successfull,y but am restricted to 18 levels of detail.  BING's resultion is higher than this, but the ESRI widget seems restricted to only show to a certain level.  Is there a way around this?  I already tried setting LODs on the map object, but the map goes blank past level 18.


texts: In PostGIS: How to delete completely overlapping polygons?

texts:
I would like to eliminate overlapping polygons from a table that is often updated with new polygons, but only if a certain attribute has a higher value.  If this higher-priority polygon only partly overlaps the lower-priority polygon, only delete the part of the lower-priority polygon that intersects.

As an analogy (not actually my case), let's say we have a table that will contain both U.S. county and state polygons.  If I have county (low-priority) polygons for, say, Virginia, and then acquire a state-level (higher-priority) polygon, I would add the state and delete the counties.

On the other hand, if I have a regional (low-priority) polygon, say New England, and add a New Hampshire (high) polygon, punch a hole in the New England polygon, insert New Hampshire, and leave the rest of the New England polygon intact.


texts: ArcGIS 10 - need a tool/script to batch save multiple MXD files to layer package

texts:
I have a directory and sub-directory with many map documents (MXD). I need to save these MXDs to map packages to send to a client. This can be done through File/Create Map Package.

Are there any tools or scripts that will take all of these MXDs and produce map packages in batch mode, preferably in the background?


texts: Importing arcgis.rand error from python shell

texts:
Not sure what's going on here with this module.  Maybe a pathing issue?

If I import the module from the python window in ArcMap it works.
If I import it in a script, and use that script in a scriping tool the tool works.
If I import it in a script and run the script, it fails.  Error message says no module named arcgis.
If I call the script tool from a Python add-in tool it fails and gives the same error message.

This worked on my computer at home but does not work on my computer at work.

Any ideas?

I'm using ArcGIS 10.1 and Python 2.7.2.


texts: How to get several property values at once so there is only one call to .`evaluate`

texts:
I have a large FeatureCollection (100,000+ polygons) and want to get several values from it at once to display in a ui.Label with only one call to .evaluate.  I have found the following way to do it that works but want to check if this is the BEST way:

var selectedstr = (ee.Feature(selectedstate.first()).select(['NAME', 'GEOID', 'ALAND', 'AWATER'], null, false)).toDictionary().values(['NAME', 'GEOID', 'ALAND', 'AWATER']);
var selectedarray;
selectedstr.evaluate(function(result) {
    if (result) {
      selectedarray = result.toString().split(",");
        mylabel.setValue(selectedarray[0] + ", " + selectedarray[1]  + ", " + selectedarray[2]  + ", " + selectedarray[3]);
        // do other stuff
    }
});

Here's a working example on a much smaller layer.


texts: Seeking resources for Python programmer to learn ArcGIS?

texts:
I know that there are several resources recommendation questions, however, they are all from GIS people trying to learn Python (for example, What are some resources for learning ArcPy?).

My situation is the opposite: I have worked extensively in Python, but know nothing about ArcGIS. What is the best way to learn?

To be more specific, I want to learn the following:

Create maps from shapefiles
Incorporate other social science-y data on map. For example, conflict site, social network, economic activity, etc.
Most importantly, statistical modeling with a spatial component (network, hotspot analysis, etc.)


texts: Processing step is ignored in QGIS3 standalone python script

texts:
In my python script based on  Qgis 3.0.3-Girona, the processing.run() is completely ignored. There are similar questions that are marked as duplicate - but those do not address or solve the issue that I list here. Here is my script.... *

from qgis.core import (
     QgsApplication,
     QgsProcessingFeedback,
     QgsVectorLayer
)
from qgis.core import QgsApplication
# See https://gis.stackexchange.com/a/155852/4972 for details about the prefix
QgsApplication.setPrefixPath('e:/Qgis3/apps/qgis', True)
qgs = QgsApplication([], False)
qgs.initQgis()
# Append the path where processing plugin can be found
sys.path.append('e:\\QGIS3\\apps\\qgis\\qtplugins')
sys.path.append('e:\\QGIS3\\apps/Qt5\\plugins')
sys.path.append('e:\\QGIS3\\apps\\qgis\\python\\plugins')
sys.path.append("e:\\QGIS3\\apps\\Python36\\Lib\\site-packages")
from processing.core.Processing import Processing
Processing.initialize()
print("Finished importing")
params = {
    'INPUT_RASTER':'E:\\arp\\LEAP\\s823_actfilt.tif',
    'RASTER_BAND':1,
    'INPUT_VECTOR':'E:\\arp\\LEAP\\test.shp',
    'COLUMN_PREFIX':'actfilt_',
    'STATS':[0,1,2],
}
feedback = QgsProcessingFeedback()
print("Starting zonalstats")
processing.run("qgis:zonalstatistics", params,feedback)
print("Finished zonalstats")

*
The output prints "Starting zonalstats" ; "Finished importing"  but never prints "Finished zonalstats". I tried several combinations including capturing the output to res variable etc. (as explained in other similar posts) but none of them work - the processing.run() is simply ignored. Not sure if it is a bug or something I'm doing wrong.


texts: plpgsql combine ST_Overlaps and ST_Intersects

texts:
I have a question regarding creating a trigger function in postgres with two conditions:

When creating a linestring in QGIS there must be either overlap nor intersect.  So I wrote two functions, each working with either ST_Overlaps or ST_Intersects.

CREATE OR REPLACE FUNCTION check_linesintersecting() RETURNS TRIGGER AS
$BODY$
DECLARE
BEGIN
IF TG_OP = 'INSERT'
THEN
IF
(SELECT COUNT(*)
FROM
(SELECT gid
FROM line_layer AS t
WHERE st_intersects(NEW.geom, t.geom)) AS foo)
THEN
RAISE EXCEPTION 'Speicherung abgebrochen: Ueberschneidende Linien!';
RETURN NULL;
ELSE
RETURN NEW;
END IF;
ELSIF TG_OP = 'UPDATE'
THEN
IF
(SELECT COUNT(*)
FROM
(SELECT gid
FROM line_layer AS t
WHERE st_intersects(NEW.geom, t.geom)
AND (t.gid <> OLD.gid)) AS foo) > 0
THEN
RAISE EXCEPTION 'Speicherung abgebrochen: Ueberschneidende Linien!';
RETURN NULL;
ELSE RETURN NEW;
END IF;
END IF;

END;
$BODY$
LANGUAGE plpgsql;

All I had to to in order to make it detect overlapping was to change the ST_Intersects into ST_Overlaps.  But how can I combine ST_Overlaps and ST_Intersects?

I've tried using the same loop right after each other but only the the first function then is working. Any clues?


texts: Google Earth Engine ee.data.authenticateViaPrivateKey in Angular app not working

texts:
I have an Angular app, developed in Angular 8 using Google Earth Engine for NDVI mapping.

Below is the basic code tried initially:

ee.data.authenticateViaPrivateKey(PRIVATE_KEY_JSON_FILE_PATH);

Since, it showed an error so tried to call this file with microservice API call

this.googleapi.authGoogleEarthEngineService().subscribe(
    (data) => {
        ee.data.authenticateViaPrivateKey(data);
        // HERE data HAS RECEIVED THE JSON FILE CORRECT. CHECKED WITH CONSOLE.LOG FUNCTION
        return data;
    },
    (err: HttpErrorResponse) => {
        // ERROR HANDLING
    }
);

I tried microservice approach because a connection of mine suggested to use server-side call to this privatekey JSON file.

Also, tried this code to see if something else isn't the issue:

        ee.data.authenticateViaPrivateKey(
            PRIVATE_KEY, () => {
                console.log('Authentication succeeded!');
                // Attempt to initialize Earth Engine.
                ee.initialize(
                    null, null,
                    // On a successful initialize
                    () => {
                        console.log('Successfully initialized the EE client library.');
                        app.listen(PORT);
                        console.log(`Listening on port ${PORT}`);
                        res.send("Authenticated and initialized!\n");
                    },
                    // On a failure to initialize
                    (err) => {
                        console.log(err);
                        console.log(`Initialization failed.`);
                    }
                );
            },
            // When the authentication failed.
            (err) => {
                console.log(err);
                console.log('Authentication failed');
            }
        );

In either case, the error is same i.e.

  ERROR Error: Use of private key authentication in the browser is insecure. Consider using OAuth, instead.

With ee.data.authenticateViaOauth(this.CLIENT_ID), it logs permission denied error when any email ID other than email ID associated with Google Earth billing is used. So, that approach also doesn't work out for me. However, I've checked Google cloud console and billing. Everything seems good.

Any guidance or resource in the right direction?

Edit:

After considering Daniel's answer, if I try with OAuth then it authenticates only with the email ID associated in the billing. Is this a normal behavior or am I missing something in the process?


texts:
I need a special projection for the purpose of effective presentation on a poster. I'd like to present the study area (Czech republic) in a context of whole Europe in such a way, that the Czech republic takes the most area (like 70%) in the center, and the rest of Europe is around it, but ridiculously small. I.e. like massive distortion in favor of the Czech republic, like a fish-eye lens, when you put them close to some object and they will massively magnify it compared to the surroundings.

Is there any projection I can use for this purpose in Quantum GIS? EPSG or the formula in standard (which can be used in QGIS) is welcome.


texts: Tracking in browser while the phone is locked with ArcGIS (ESRI)

texts:
The Question

How can I track a user's movements via ESRI's JavaScript API in a web browser when their phone is shut/locked?

This would ideally work for all web browsers (i.e. Android and Apple versions of Chrome, Firefox etc..). However, any platform specific suggestions would also be very welcome.

What I've tried

I've thought about playing a video on loop (hidden) on the maps page so the phone doesn't automatically lock, keeping the browser open. However, this would potentially be a huge waste of battery for users and false clicks may get recorded while in the user's pocket so I'd like to avoid this solution if possible.

Some more context:

I'm developing a charity webapp for a community group trying to track their movements when out on patrols which can last between 30 mins and a couple of hours.

The tracking should work similarly to strava and should record their movements at regular intervals (I've decided on every 10 secs but this is somewhat arbitrary).

I've already written the majority of the app including the backend/database code and deployment scripts so I really wouldn't want to start again creating a stand-alone Android and iOS app.