texts: How to get value of uvIndex for a particular latitude and longitude from a grib2 file?

texts:
I'm trying to find out the uvIndex of a particular set of lat,lng. I'm getting the data file from this link. I'm using pygrib to read the file.

import pygrib
grbs = pygrib.open('uv.t12z.grbf01.grib2')
grb = grbs.select(name='UV index')[0]
print grb.data(23.5,55.5)

This prints 3 numpy arrays. But I'm not able to figure out how to extract data from it.


texts: Will Topographic Position Index (Tpi) not work in scattered raster layers?

texts:
I have a non-continuous elevation raster layer and want to use tpi in the
spatialEco package in R.

The elevation raster is in EPSG:5677 and looks like this:

I always get this (to me) rather cryptic error:

library(spatialEco)
tpi(elevation.raster,scale=150)
Error in .local(x, ...) : w must have uneven sides

I wonder if this is due to the raster layer being uncontinuous?


texts: Connecting Geoprocessing Service to ArcSDE?

texts:
Recently we updated from ArcMap 10.1 to 10.2.2. ArcGIS Server is still on 10.1. We have numerous GP and Mapservices that were published on ArcMap 10.1 and are still running without problems.
I also published a Map service from ArcMap 10.2, with data in an Enterprise GDB that is registered, without any problems.
Now I have a python script that uses data from the same DB and works perfectly in ArcMap 10.2. I can publish it but when I run it, it returns false on the first

arcpy.Exists(sdeconnection\\FeatureClass).

According to this link I should however be able to publish form 10.2.2 to 10.1 without too many problems.

Any ideas on what could be wrong? I guess it should be the data connection. But my Map service uses the same SDE, so my data should still be accessible.

In the original file, my workspace was defined as follows:

arcpy.env.workspace = "\\\\v057GIS\\Data\\Connection to v057gis.sde"

After publishing, this line changed suddenly to:

arcpy.env.workspace = arcpy.env.packageWorkspace

I can open the "C:\arcgisserver\directories\arcgissystem\arcgisinput\GPservice.GPServer\extracted\v101" file in python IDE, change the line and run it correctly. I'm however not able to save it.
I read Modifying script published as Geoprocessing Service in ArcGIS for Server? about how to do it, but am a bit afraid of doing it since we only have production server.

Is there a workaround for that? Maybe directly in packageWorkspace?


texts: Problem using qgis when I upgraded to osx 10.8.1

texts:
I have been running qgis in Mountain Lion for a few weeks now, but when I upgraded to OS X 10.8.1, qgis would no longer work.  I get an error saying:

QGIS cannot be opened because of a problem.
Check with the developer to makes sure QGIS works with this version of OS X. You may need to reinstall the application.  Be sure to install any available updates for the application and OS X.

I have tried re-installing qgis (I have even tried installing the nightly build) and I have also updated all of the packages that qgis requires.

Has anyone had a similar problem or more importantly, does anyone have a solution?

Thank you!


texts:
I have inherited data that is in need of the lines to be in the same place. Unfortunately this is not how the data was made. I have been moving these by hand.
Is there a tool that will arbitrarily pick the center of the options and make one cohesive line here?
I have tried the snapping tool but it is not fixing lines like these.


texts: How to set QGIS Python path so that I can use system packages on Windows?

texts:
Several Python packages are best installed with a Windows installer. But those installers install to the system-wide Python directory. How do I set the Python path in QGIS?

I tried setting PYTHONPATH as a Windows environment variable. It can import netCDF4 (which is installed globally), but the import doesn't succeed because it cannot load a DLL. So perhaps also some DLL load path needs to be adjusted?


texts:
I need the Evapotranspiration data for my study (The study area is located in Australia).
Somebody told me maybe I can use MOD 16 MODIS actual ET estimates. I could not find this data on earth explorer. Instead  I have found the following link http://www.ntsg.umt.edu/
which provide MOD 16 MODIS data. There are several folders in data folder which contain *.hdf files. Does anyone knows which folder I have to choose and how I can open the *.hdf files?(I tried to open one of *.hdf files in Arcmap. 3 bands are blank and one band is completely grey) .


texts: Processing algorithm doesen't work in standalone script - PyQGIS 3.4

texts:
I am trying to execute the processing algorithm extractbylocation. Anyway, firstly I run that algorithm in QGIS3 and then just copy the code from the processing history and changes paths to variables. I have some features with invalid geometries (in QGIS in the settings it is set for processing to ignore that features),so that is maybe my problem. In order to avoid that I deleted all the features with invalid geometries and the problem is still there.
So my code looks like this:

input_path = 'C:/Users/file.shp'
output_path = 'C:/Users/file_out.shp
crs = QgsCoordinateReferenceSystem("epsg:4326")
input = QgsVectorLayer(input_path, "links", "ogr")
input.setCrs(crs)

#deleting features with invalid geometries
for f in input.getFeatures():
    geom = f.geometry()
    if geom.isGeosValid():
        print('valid')
    else:
        input.deleteFeature(f.id())
        print('not valid')

sys.path.append('C:/Program Files/QGIS 3.4/apps/qgis/python/plugins')
import processing
from processing.core.Processing import Processing
Processing.initialize()
QgsApplication.processingRegistry().addProvider(QgsNativeAlgorithms())

processing.run("native:extractbylocation",
               {'INPUT': input, 'PREDICATE': [1],
                'INTERSECT': 'C:/Users/points.shp', #I also tried to put here the QgsVectorLayer object instead of path
                'OUTPUT': outputh_path})

I got the error next to the output line in the processing algorithm.

 File "path_to_my_file_py", line 49,    'OUTPUT': outputh_path})
 File "path_to_my_file_py", line 84, in run
return Processing.runAlgorithm(algOrName, parameters, onFinish, feedback, context)
 File "path_to_processing_file_py", line 194, in runAlgorithm
raise QgsProcessingException(msg)
_core.QgsProcessingException: There were errors executing the algorithm.

And btw this code works at python console in qgis.


texts: Inserting custom latitude-longitude data in OpenLayers Heatmap

texts:
I need to be able to insert custom latitude and longitude data into my heatmap, and I thought I found a way by following this previous answer.

However, the point is always centered at the same location, regardless if I change the point's lat-lon. Below is a snippet of my code:

// Map
var map = new ol.Map({
  layers: [
    new ol.layer.Tile({
      source: new ol.source.OSM()
    }),
  ],
  target: 'map',
  controls: ol.control.defaults({
    attributionOptions: /** @type {olx.control.AttributionOptions} */ ({
      collapsible: false
    })
  }),
  view: new ol.View({
    center: [0, 0],
    zoom: 2
  }),
  interactions: ol.interaction.defaults().extend([
    new ol.interaction.DragRotateAndZoom(),
    new ol.interaction.PinchRotate(),
  ])
});

// created for owl range of data
var data = new ol.source.Vector();
var coord = [2.1734, 41.3851];  // Barcelona, Spain
var lonLat = new ol.geom.Point(coord);
var pointFeature = new ol.Feature({
   geometry: lonLat,
   weight: 20 // e.g. temperature
});
data.addFeature(pointFeature);

// create the layer
heatMapLayer = new ol.layer.Heatmap({
   source: data,
   radius: 50
});

// add to the map
map.addLayer(heatMapLayer);

I cannot seem to figure this out. I was able to create a proper heatmap following the example provided by OSM here, but that would require me to load a KML file. However, I need to be able to use entries in my database to populate this heatmap instead. Below is a picture of what the heatmap looks like:

Any ideas? Does it have something to do with needing to transform the map view so the coordinates are read correctly? I've been wracking my brain over this for days.


texts: Gathering data for large number of census block groups?

texts:
I would like to collect some census data for every block group in New England.  I'm primarily interested in median income and median age, and ideally I'd like a large .csv file with one row per census block group and one column per variable (e.g. id variables, median_income, etc.).

American FactFinder seems almost useless for this kind of thing.

Does the census have a "bulk data export" option?


texts: How can I use WFS from Geoserver with OpenLayers 3?

texts:
I have tried with different setups and formats, but no success.

With firebug, I can see that the json data is loaded, but nothing shows up on the map.

Have anyone experience with geoserver and wfs in openlayers 3?

var layer_ojd_vector = new ol.layer.Vector({
 source: new ol.source.Vector({
  parser: new ol.parser.GeoJSON(),
  url: http://SERVER:8080/geoserver/NS/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=NS:LAYER&maxFeatures=50&outputFormat=json
 }),
 style: new ol.style.Style({rules: [
  new ol.style.Rule({
   filter: 'geometryType("point")',
   symbolizers: [
     new ol.style.Shape({
       fillColor: '#013',
       size: 40
     })
   ]
 })
]})
});


texts:
I'm trying to upload a .pdf file in the content of my Portal for ArcGIS with a python script. I'm trying with the "additem" REST operation. My code is as follows:

import urllib
import urllib2
import json
import argparse
import os
import json

sourcePortal='https://xxx/arcgis'

def generateToken(username, password, portalUrl):
    '''Retrieves a token to be used with API requests.'''
    parameters = urllib.urlencode({'username' : username,
                                   'password' : password,
                                   'client' : 'referer',
                                   'referer': portalUrl,
                                   'expiration': 60,
                                   'f' : 'json'})
    response = urllib2.urlopen(portalUrl + '/sharing/rest/generateToken?',
                              parameters).read()
    try:
        jsonResponse = json.loads(response)
        if 'token' in jsonResponse:
            print("Token generated OK")
            return jsonResponse['token']
        elif 'error' in jsonResponse:
            print jsonResponse['error']['message']
            for detail in jsonResponse['error']['details']:
                print detail
    except ValueError, e:
        print 'An unspecified error occurred.'
        print e

def uploadFile(token):

    '''Creates a new item in a user's content.'''
    parameters = urllib.urlencode({'title' : 'myreport',
                                   'overwrite': 'true',
                                   'type' : 'PDF',
                                   'token' : token,
                                   'async' : 'true',
                                   'multipart' : 'true',
                                   'filename' : 'newfile',
                                   'f' : 'json'
                                   })
    response = urllib2.urlopen('https://xxx/arcgis/sharing/content/users/admin/additem',
        parameters).read()
    return response

sourceToken = generateToken(username='xxx', password='xxx',
                            portalUrl=sourcePortal)

reply=uploadFile(sourceToken)
itemID = json.loads(reply)["id"]

pdfFile = open('C:\\xxx\\xxx.pdf')

parameters = urllib.urlencode({
    'file' : pdfFile,
    'partNum': '1',
    'token' : sourceToken,
    'f' : 'json'})
print(urllib2.urlopen('https://xxx/arcgis/sharing/content/users/admin/items/'+ itemID + "/addpart",parameters).read())

parameters = urllib.urlencode({
    'token' : sourceToken,
    'f' : 'json'})
print(urllib2.urlopen('https://xxx/arcgis/sharing/content/users/admin/items/'+ itemID + "/commit",parameters).read())

I have a problem with the "addpart". Here is the message :

{"code":400,"message":"Unable to add part.","details":["Invalid or missing 'file' parameter."]}

Does anyone have an idea to figure this out?


texts: Downloading or caching web based map service for faster access?

texts: The Western Australia Department of Mine and Petroleum allows access to one of their online systems called Tengraph Web with which I can view tenements, major infrastructure, and topography for the state of Western Australia. There is an option to overlay "Landgate Aerial Imagery" which is basically a satellite view of the same area however none of the areas I scroll over are cached thus every time I move my window, the aerial view has to load it again. With poor Australian internet, I was wondering if anyone had ideas as to either download or cache the aerial view data so I can navigate more easily without having to wait a few minutes for every kilometre of land to load.


texts: Legend load event in ArcGIS JavaScript API

texts: I am working on arcgis map JavaScript API and i want to show map layer legend with close icon but i am unable to find legend load event to show the same. Legend along with close icon is shown in below image.


texts: Is there a way to disable code autocomplete in the Python window in ArcGIS 10.0

texts:
Autocomplete slows down my workflow in the Python window, especially if I am remoted into a machine.

Is there a way to disable this?

Alternatively, a way to disable this and then invoke as needed (similar to Intellisense options in Visual Studio)?


texts: Creating Complex Polygon from Point Layer using only Boundary Points in ArcGIS Desktop

texts:
I need to convert a point layer to a polygon, using the boundary points from a complex grid to define the polygon's edges.

I need to incorporate this into a ModelBuilder framework in ArcGIS Desktop 10.3. Iterating the process will be required (if possible) because of lots of incoming data.

The point layer is gridded over a river segment, and I need to determine the rivers boundary  points, and connect them to create a polygon layer of the river segment.

Convex hull doesn't seem to work with how the rivers meander, I need a clean tight boundary, not a containment like convex hull does. I do have layers for just the boundary points, but I don't know how to connect them to get to a polygon.


texts: Spatial Query: The query from "LayerName" using "FID" in field not possible

texts:
In QGIS Wien, I am performing a spatial query with two layers (points and buffer areas) and receive the following error when creating a new layer with intersecting features: The query from "LayerName" using "FID" in field not possible.

Spatial query settings:
- where the feature: intersects
- and use the result to: create new selection

In the same file, I have performed the same spatial query with other layers (points and buffer areas) without running into this problem.

How should I go about fixing the error?


texts: Different raster shapes even though I used the same clipping polygon (Python)

texts:
I have a bathymetry raster image that I would like to use for masking regions where elevation > -1. I have started by clipping the bathymetry image and another target image using the same rectangular polygon. However, when checking the shape of clipped images, results seem to be slightly different:
img_mask = rasterio.open('\Bathymetry\spainMask\mask_1_ESP.tif')
img_mask.shape

output: (972, 1159)

img = rasterio.open("\Rw490_ESP.tif")
img.shape

output: (972, 1160)

This slight difference in shape makes it impossible to apply the bathymetry mask subsequently. Is there a way to avoid this problem or to detect the pixels that do not exist in both images and remove them?
Links to clipped images:
https://drive.google.com/file/d/1C9SLT3mCDuLZC8ce508w70EXXguJCysu/view?usp=sharing
https://drive.google.com/file/d/1L0wlljxBXRsnjDf5rFClqEI4hfsMqOeL/view?usp=sharing
The code used for clipping the images:
ESP = gpd.read_file("\Vector\spainVector.shp")
geometry_ESP = ESP["geometry"]

clipImg, clipImg_Transform = mask(img, geometry_ESP, crop=True, invert=False)
clipImg_Meta = img.meta.copy()
clipImg_Meta.update({"driver": "GTiff",
                    "height": clipImg.shape[1],
                    "width": clipImg.shape[2],
                    "transform": clipImg_Transform})
outputDirectory = '\bathymetry_ESP.tif'

with rasterio.open(outputDirectory, 'w', driver='Gtiff',width=clipImg.shape[2],height=clipImg.shape[1], count=1,
                    crs=img.crs,transform=clipImg_Transform, dtype='float32') as dst:
    dst.write(clipImg)


texts: AssertionError 25 on getExtent, Openlayers?

texts:
When doing this.map.getView().fit(this.layer.getSource().getExtent()) with my vector layer I get an assertion error as following:

  AssertionError {code: 25, name: "AssertionError", message: "Assertion
  failed. See https://openlayers.org/en/v5.3.3/doc/errors/#25 for
  details."} code: 25 message: "Assertion failed. See
  https://openlayers.org/en/v5.3.3/doc/errors/#25 for details." name:
  "AssertionError"

This is the value received from the server and applied to the vector source:

 {
   "type":"FeatureCollection",
   "totalFeatures":1,
   "features":[
      {
         "type":"Feature",
         "id":"CATAST_Pol_Municipio.201",
         "geometry":{
            "type":"MultiPolygon",
            "coordinates":[
                //Long text
            ]
         },
         "geometry_name":"the_geom",
         "properties":{
            "FEATURE":200020,
            "CMUNICIPIO":201,
            "MUNICIPIO":"Pamplona / Iruña",
            "MUNINOAC":"Pamplona / Iruña",
            "GEOM_AREA":2.509439497E7,
            "GEOM_PERI":36044.89,
            "BEGINLIFE":"04/10/2019"
         }
      }
   ],
   "crs":{
      "type":"name",
      "properties":{
         "name":"urn:ogc:def:crs:EPSG::25830"
      }
   }
}

If I do this.layer.getSource().getExtent() I get:

[Infinity, Infinity, -Infinity, -Infinity]

  How can I solve this error or get the correct extent with the received
  values from the server?


texts: QGIS - does a tool exist to calculate the (PostGIS) polygons containing a given point?

texts:
This is the converse of the analysis provided by the fTools plugin "Points in Polygon":

Given a single point and a list of vector layer polygons (PostGIS in my specific case), is it possible to list only those layers containing the point? I have a number of 'zones' mapped as PostGIS vector layers and need a simple tool to determine which zones any given point is in.

Question #172237 poses the general question for a Java or Python solution and an answer to question #27005223 over on regular Stack Overflow provides the SQL to do it for PostGIS data.

I'd guess that this is a very commonly required spatial analysis, so I thought I'd check it didn't already exist as a QGIS plugin before going down the SQL/coding route myself and potentially reinventing the wheel.


texts: Find Minimum Bounding Box for a set of coordinates

texts:
I am working with GPS tracks (lat/lon values) and my goal is to compute the minimum bounding box for a small geospatial library I am working on. The datasets are not huge, maybe a couple of tens of thousands of points. I could not find an algorithm that deals with tracks/datasets crossing from -180 to 180 lat or lon. Does anybody have pointers to algorithms?

To sum it up very simply, what I am looking for is:

Input: (lat/long)* [possibly tens of thousands of GPS points, potentially crossing 180/-180 boundary)]
Output: (lower-left, upper-right) of the minimal box containing all points

The background for me is that I need the bounding boxes to initialize map viewers, in particular when there is more than one track displayed on a map.


texts: I have 2 datasets of raster footprints (Luc_1,  Luc_2, Luc_3……. L_250) and  (Dir_1, Dir_2, Dir_3 …. Dir_250). I need to combine Luc_1 with Dir_1, Luc_2 with Dir_2, Luc_3 with Dir_3, and so on. I have created a Model (using the Iterate Raster function) combining Luc_1 with all Dir rasters, Luc_2 with all Dir rasters, Luc_3 with all Dir rasters. It's a bit more time-consuming and in the end I have to delete the 249 rasters not needed. Is there any way I can specifically pair only the rasters I am after?


texts: VB Script or Python for Field Calculator when using greater/less than?

texts:
What I have:
I have 4 fields for my address ranges.  L_ADD_FROM, R_ADD_FROM, L_ADD_TO & R_ADD_TO (all in DOUBLE field type).

What I need:
The lesser value from L_ADD_FROM and R_ADD_FROM to populate in a field called FROMADDRESS.

I've searched extensively, trying Python and VB Script without much luck.  What's the best method for this?


texts: Layout misshapen when export to PDF in Arc GIS 10

texts:
I want to include an Excel table in my ArcMap Layout which is working fine as long as I won't print to PDF. If I export the map to PDF the table looks like this:

There are a lot of weird characters and I have no clue why it looks so different to the original table which is still correct in ArcMap. Does anybody knows a solution for this problem?


texts: How to close open files after processing.runlag() command (Qgis: Error No 24)

texts:
I'm creating a pyqgis script. In this script I use the processing.runalg command several times. The problem is that these processes are within a for loop. After several iterations, the error number 24 "too many open file" appears. How can I close the files generated by the processes? Here is an example of the code:

  for i in range (1, 13000):

path1 = '/path/to/file1_' + str(i) + '.shp'
path2 = '/path/to/file2_' + str(i) + '.shp'
path3 = '/path/to/file3_' + str(i) + '.asc'
path4 = '/path/to/file4_' + str(i) + '.asc'
path5 = '/path/to/file5_' + str(i) + '.shp'
path6 = '/path/to/file6_' + str(i) + '.shp'
path7 = '/path/to/file7_' + str(i) + '.asc'
path8 = '/path/to/file8_' + str(i) + '.asc'

processing.runalg('qgis:reprojectlayer', path1,'EPSG:3004',path2)
processing.runalg('qgis:reprojectlayer', path5,'EPSG:3004',path6)

processing.runalg('saga:cubicsplineapproximation', path2,'field',3.0,20.0,5.0,140.0,'2596187,2675187,4321914,4403914',1000.0,0,path3)
processing.runalg('saga:cubicsplineapproximation', path6,'field',3.0,20.0,5.0,140.0,'2596187,2675187,4321914,4403914',1000.0,0,path7)

processing.runalg('gdalogr:cliprasterbyextent', path3,'-9999','2611187.7500000000000000,2660187.7500000000000000,4336914.1990000000223517,4388914.1990000000223517',5,4,75.0,6.0,1.0,False,0,False,None,path4)
processing.runalg('gdalogr:cliprasterbyextent', path7,'-9999','2611187.7500000000000000,2660187.7500000000000000,4336914.1990000000223517,4388914.1990000000223517',5,4,75.0,6.0,1.0,False,0,False,None,path8)
...     ...


texts: How to update the OSM data in a PGRouting Database built with osm2pgrouting?

texts:
I'm using a PGRouting database originally built with osm2pgrouting from a *.osm file downloaded from Mapzen.

If I want to update the OSM data in that database, ie. re-build my PGRouting Database, what's the best workflow to accomplish this?

Also: I'm not married to the manual download of the *.osm file, so a solution involving osm2pgsql would be fine with me too.


texts: Using QGIS (or PostGIS), determine which network featureIds are shortest path between two points

texts:
I have a network layer, and a layer of target start/end points (both as coordinate sets and as start feat_id and end feat-id).
like:
featid, geom
and
start_coord_x, start_coord_y, start_featid, end_coord_x, end_coord_y, end_featid
(although really "start" and "end" don't matter for this case, and I could use either the input coord sets or the respective featids)
I'm trying to determine the shortest path featid's that connect the two points. I have found many QGIS tools that will give me the shortest path as associated with cost, but none I have tried report back the exact featids used to create that shortest path.
Has anyone found a tool that will do this?


texts: Changing vertex order for polygons in QGIS?

texts:
I have this intersection of a polygon, see image below

I would like to redraw the lines so they don't cross each other.

So, in the end, I would like to have a polygon without self-intersecting parts. How can I accomplish this?


texts: How to create bilingual composite address locator in ArcGIS 10?

texts:
I'm fairly new to geocoding in ArcGIS and stumbled upon first problem.

I'm working with GeoPost data from Swiss Post. As long as I'm using their default locator everything works flawlessly.

I also received additional table with alternative street names (French <-> German) that I can link to main data by unique street ID. They exist for small subset of main dataset. However, there is no way to add alternative street names using the Swiss Post locator.

I browsed through ESRI's docs but it didn't help much in this case.

I've tried building two separate locators (one with German names, one with French) and wanted to combine them into composite locator. Unfortunately I get general error message without any further pointers :/

Is there a way I could use this additional data to build composite address locator that will use both language versions of street names? Is there any other potential solution?


texts:
I would like to know some examples of GIS projects that have gone wrong.

i.e. Bad decisions would/were made, based upon the GIS output from that project.

I am not looking at naming and shaming anyone, so if it is a company you have worked for, please keep any names out of it, and disguise the story where necessary.  I am sure there must be quite a few major public examples?

e.g. A disaster mgmt project, gone wrong due to bad or out-of-date data.

I admit, it is a fairly subjective question with no right answer.
Therefore, I will let the votes do the talking, and if someone could convert to c-wiki if need be.

Personally, I will get great benefit from this, and I will be passing the examples on when teaching GIS.


texts: Exporting 2D vector graphic from ArcScene to Illustrator

texts:
How can I export 2d vector graphic from ArcScene so I can edit it with Adobe Illustrator?
I've tried changing the options when I export 2D scene as AI file or PDF file, but they always remain as a unique image that does not allow me to edit the lines and polygons.


texts: The relation between the Feature's Geometry and its externalGraphic

texts:
I wonder what is the relation between the Geometry.Point one sets for a Feature.Vector and the image provided as externalGraphic for the style, in OpenLayers.

Is the Geometry.Point of a feature for the image:

The centre point of the image,
or the left top most corner of the image.

In case it is 2 then please advice how would I set the values of offsets if, for example, the image is 32x32 (WxH). Is it simply 16 (W/2) for the value of graphicXOffset, and accordingly for the graphicYOffset?

I need to know this information to set the image centre to be the point of a feature. Sadly the documentation doesn't say a word about it, or I am unable to find it.


texts:
I just noticed some strange behavior in PostGIS, with respect to WKT parsing.

A GeometryCollection containing 1 or more empty geometries

postgis_test=# SELECT ST_GeomFromEWKT('GEOMETRYCOLLECTION(POINT EMPTY, LINESTRING EMPTY)');
 010700000000000000

is the same an empty GeometryCollection

postgis_test=# SELECT ST_GeomFromEWKT('GEOMETRYCOLLECTION EMPTY');
 010700000000000000

This doesn't seem right. But here's the truly weird part. If you throw a non-empty geometry into the collection, the empty geometries are suddenly represented in the collection:

postgis_test=# SELECT ST_GeomFromEWKT('GEOMETRYCOLLECTION(POINT EMPTY, LINESTRING EMPTY, POINT EMPTY)');
 010700000000000000
postgis_test=# SELECT ST_GeomFromEWKT('GEOMETRYCOLLECTION(POINT EMPTY, LINESTRING EMPTY, POINT (1 2))');
 0107000000030000000104000000000000000102000000000000000101000000000000000000F03F0000000000000040

It seems strange that empty geometries are only recognized if the collection contains at least one non-empty.

Is this is a bug? Is this behavior documented somewhere?


texts: Interpolation with z-values taken from a related table

texts:
This is a follow-up question to organizing attribute table: multiple sets of variables per point.

I am working with ArcGIS 10.1, advanced license. My aim was to build a points layer, each point being a weather station in the study area. For each weather station, some data was available, consisting in the winds speed, average direction, air temperature, and the like.

I have created two related tables: one with the actual station points, the other with the data (i.e., variables) organized per month.

Now, I wish to perform some sort of interpolation on, say, the wind speed using, e.g., IDW. I am required to feed into the arcgis' function the feature dataset and Z value field. The problem is that the latter is stored in the related table, which is different from that of the feature class (i.e., station points).

How am I supposed to go on this?


texts: Calculating km² represented by same pixel value in a 1km² resolution raster?

texts:
From a raster resolution of 1km² we want to calculate the area that represents each pixel of the same color. The calculation must be done in QGIS.

Update : what parameters should I insert here:


texts: FWTools ogr2ogr command does not seem to recognise dim parameter

texts:
I am trying to convert a 3D shapefile to 2D via FWTools shell.

This is the command I am trying

ogr2ogr REGC2013_GV_Clipped_lonlat_2d.shp REGC2013_GV_Clipped_lonlat.shp -dim 2

Which returns the Usage - there is no mention of the dim parameter in the usage. It does not output a new shapefile. I have managed to convert a shapefile from map projection to lonlat coordinates using this

ogr2ogr REGC2013_GV_Clipped_lonlat.shp -t_srs EPSG:4326 REGC2013_GV_Clipped.shp

but if I add the -dim 2 it does not work and I get the usage returned.

I am using Windows 7 and this is the version info

C:Users\xxxxx\Desktop\Plotting>gdalinfo --version
GDAL 1.7.0b2, FWTools 2.4.7, released 2010/01/19


texts: Automatically calculate average slope of a street segment

texts: I'd like to be able to automatically calculate the average slope of a street segment in ArcMap (start elev - end elev / length). My map looks something like this, with street centerlines in black and 2-foot contour lines in green. The blue polygons represent street drainage areas I have drawn, which I'd eventually like to assign the slopes of the adjoining centerlines to as an attribute, but can't even figure out how to automatically find the centerline slope in the first place. Any thoughts or past experience?


texts: ArcMap "datum transformation is expected, but the parameter is empty" projection error

texts:
I have a feature class in GCS North American 1983 that I need to convert to US National Atlas Equal Area. When I try to use the Project feature in ArcToolbox, after selecting US_National_Atlas_Equal_Area as the output coordinate system, I get this error:

  Warning 001839: A datum transformation is expected, but the parameter
  is empty.


texts: Using Leaflet: How do I get it to display one popup based on a condition?

texts:
Using Leaflet: How do I get it to display one popup based on a condition.  See below.

    basemap_0.addTo(map);
    var layerOrder=new Array();
    function pop_thunder(feature, layer) {

    if (feature.properties.CODE = 0) {

    var popupContent = feature.properties.html1_exp;
    }
    else {
    var popupContent = feature.properties.html_exp;
    }
        layer.bindPopup(popupContent);
    }


texts:
I set up a script for supervised following the model from the EE tutorial -- Only thing I changed was to substitute Landsat 5 with Sentinel 2. The classifier does not give me any result -- Any idea why or how to fix it?

var point = ee.Geometry.Point([18.31, 20.85]);
var s2 = s2toa.filterBounds(point) //filter all scenes using point geometry from above (i.e. limit to Friuli)
  .sort('CLOUDY_PIXEL_PERCENTAGE', true).limit(300);
var s2Time1 = s2.select('B[2-8]');
var input = ee.Image(s2Time1.first()).clip(studysite);
print(input)
Map.centerObject(studysite, 10);
Map.addLayer(input, {
  min: 750,
  max: 2100,
  bands: ['B8', 'B3', 'B2']
}, 'input image');
var training = basalt_train.merge(Sand_train);
var training = input.sampleRegions({
  collection: training,
  properties: ['class'],
  scale: 30
});
var classifier = ee.Classifier.smileRandomForest(100).train(training, 'nd'); //set the number of trees in random forest to 100
var classified = input.classify(classifier);
print(classified)
Map.addLayer(classified, {
  min: 0,
  max: 1
}, 'classified image');

https://code.earthengine.google.com/751ac9cfdb5c33ad5790bf69bae87003


texts: Complexity restrictions for Google Maps kml

texts:
What is meant by "Maximum number of total document-wide features = 1000"?
http://code.google.com/apis/kml/documentation/mapsSupport.html

What is a feature in this regard?


texts: How to set -srcwin with python gdalbinding and driver.createcopy

texts:
I am using gdal_translate currently to cut a 60000x60000 geotiff into smaller 1000x1000 geotiff files. This works but it takes about 24 hours(and I have 44 of those big images to compute).

Im not sure if it would be faster but my theory is that using the gdal binding would be faster if it just loads the image into memory 1 time, then creates the new geotiff from that.(instead of running the gdal_translate command about 3600 times in series).

So the real question is, can I do that with the driver.createcopy command? Or do I have to use the driver.create method ? And would this be faster then actually starting gdal_translate 3600 times?

Currently I'm using gdal_translate like this from a batch file:

gdal_translate -of GTiff -srcwin 0 0 1000 1000 -co TFW=YES -co COMPRESS=JPEG -co PHOTOMETRIC=YCBCR -co JPEG_QUALITY=90 sourcefile.tif outfile.tif

I have seen that you can set the creation options like TFW=YES and PHOTOMETRIC=YCBCR with driver.createcopy but havent found any info on selecting a subwindow.


texts: Using grid in globe map indicator of ArcGIS Pro?

texts:
I'm trying to reproduce a map indicator such as:

But in ArcGIS Pro, grids are only available to 2D maps. So far I've got this:

Is there a way to add a grid in 3D maps without adding an extra layer (such a shapefile with CRS's grid net) or a way to visualize 2D maps as a globe?


texts: PG-Strom and spatial queries in PostgreSQL

texts:
I would like to start using the GPU on my server in my PostgreSQL queries. I cannot find much information about the use of PG-Strom for spatial queries. The most I have found is that it may be coming in PostgreSQL 9.6. We are beyond that now. So? What's happening with this? Where are the resources?

http://on-demand.gputechconf.com/gtc/2015/presentation/S5276-Kohei-KaiGai.pdf

https://github.com/heterodb/pg-strom

I am working with PostgreSQL 12.1, PostGIS 3 with a VOLTA Tesla V100 16GB


texts: Automate Projection of multiple shp files using ArcPy

texts:
I have been trying to create a Python script that will convert all shape files within a folder from GCS WGS 84 to Web Mercator.  Below is the code I have so far but what I'm struggling with  is how to iterate thru each of the files while creating a new name for the output (the new projected file in Web Mercator).  Any ideas on how I could change this code to make this work? I know I don't have the output parameter right but I think everything else is right.

import arcpy
import glob
inFCFolder = "C:\Data"
def projectToWebM():
    for f in glob.glob(r""+ inFCFolder + "\*.shp"):
        try:
            arcpy.Project_management(f, f + r"\Project\*.shp","PROJCS['WGS_1984_Web_Mercator_Auxiliary_Sphere',GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Mercator_Auxiliary_Sphere'],PARAMETER['False_Easting',0.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',0.0],PARAMETER['Standard_Parallel_1',0.0],PARAMETER['Auxiliary_Sphere_Type',0.0],UNIT['Meter',1.0]]","#","GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]")
        except:
            print 'Error: Unable to project ' + f

if __name__== '__main__':
    projectToWebM()


texts: Geocode from mile marker/milepost to latitude/longitude?

texts:
I have a dataset that describes a few million incidents along Texas roadways. The vast majority of incidents only list a highway name and a milepost. For example, milepost 120 on I-20 or milepost 263 on US 60.

How would I convert these milepost/highways to latitude and longitude?

This is for a low-budget academic project, so trying to stick with free tools as much as possible.

My Google searching hasn't turned up much fruit, although some indications are that TxDOT may have data I could request through the Texas Public Information Act, but not sure yet.


texts: Scripting color palette to align with corresponding classes from supervised classification in Google Earth Engine

texts:
When I display my classification map, the colors don't line up with the variable that I want them to. For example, my forest category is represented as red (instead of my goal color-green), or my water category shows up as gray (instead of 'blue').

How do I order my color palette such that it corresponds with the classes appropriately?

This is my script so far in Google Earth Engine:

//Create a classification map. Classify the input imagery
var classified = composite.select(bands).classify(classifier);

//Define the palette for the land use classification

var palette = [
  '0000FF', // water (0)  -- blue
  'FFFF00', // cultivated (1)  -- yellow
  'DCDCDC', //  Settlement (2) -- gray
  'F0E68C', // Grasslands (3)--khaki
  '228B22', // forest  (4)-- forest green
  '000000', //clouds (5)-- black
  '808000', //intact veg (6)-- olive
  'FFA500', //impacted_veg (7)-- orange
  'DC143C', //mines (8)-- crimson
];

// Display the classification result and the input image.
Map.centerObject(K2C_aoi);
Map.addLayer(classified, {min: 0, max: 2, palette: palette}, 'Land Use Classification');


texts: ArcGIS ModelBuilder output data from collect values as table or csv

texts:
I have 90 point shapefiles. I would like to export the largest value of one of the fields in the attribute table. I have sorted the field so that the largest number is at the top. I have also added in another field, which is the name of the shapefile so I can export both fields and still know the name of the file.
I am now trying to use 'Get Field Value', 'Collect Values' and an array of tools to get the data into a table or csv or similar which I can work with. see image: I don't have to do both at the same time, I just want to be sure which number comes from which file.
I am not great with Python so if this is possible within ModelBuilder can you help?


texts:
I'm running QGIS 2.12 and have the same issue on 2.8, but it is now worse in this latest version. The issue is that not underlying raster layers (satellite photos loaded in as .tif) are being displayed. However moving the map, or changing the zoom level, can usually fix these particular layers but then others are not displayed.(see screenshots) This seems to me like a limit to the number of raster layers that can be shown at any time. However, this issue never occurred with older versions of QGIS.

Any suggestions how to fix this?

EDIT: Screenshots here: http://imgur.com/a/nzIT7


texts: How is the pground metric (to measure percent ground hits) implemented into a function for grid_metrics in lidR?

texts: I'm having difficulty with the method to obtain the percentage of points classified as ground using grid_metrics(). I would like to implement pground (one metric) like stdmetrics (a set of metrics), but can't see how. Does it have to be set up in a custom function instead? (reference I'm using: https://www.rdocumentation.org/packages/lidR/versions/3.1.0/topics/stdmetrics). thanks.


texts:
I've used QGIS for many years and made great flow maps.

But within QGIS 3 I can't find any flow making plugin like the ones in QGIS 2.X versions.

Do you know how to do?


texts: Error in merging two shapefiles of different geometry types

texts:
I am trying to merge two shapefiles...

Here are the links of two shapefiles

Australian State & Territory

State and Territory ASGC Ed 2011 Digital Boundaries in ESRI Shapefile Format (http://www.abs.gov.au/ausstats/subscriber.nsf/log?openagent&1259030001_ste11aaust_midmif.zip&1259.0.30.001&Data%20Cubes&6E45E3029A27FFEFCA2578CC0012083E&0&July%202011&14.07.2011&Latest)

Australian Suburbs

State Suburbs ASGS Non ABS Structures Ed 2011 Digital Boundaries in ESRI Shapefile Format
(http://www.abs.gov.au/AUSSTATS/subscriber.nsf/log?openagent&1270055003_ssc_2011_aust_shape.zip&1270.0.55.003&Data%20Cubes&D68DFFC14D31F4E1CA2578D40013268D&0&July%202011&22.07.2011&Previous)

When I try to merge with QGIS extension MMQGIS.

But it is giving me following error. I try to Google it but can't find something that a newbee can understand.

Error!! Merged layers must all be same type of geometry (polygon != polygon 2.5D)

Can someone explain what exactly this error is and how to fix the problem and merge two files?


texts: Determining the area in Sq Miles of 1 polygon that are contained completely within another polygon

texts:
I am looking to determine how much of the area of 1 polygon (Lets say a body of water for instance) falls completely with in another polygon (Call it a state or county for instance) i have a layer that contains all of the area covered by water and another layer that has all of the county boundaries.  I have selected the specific county that I am interested (creating a new "table" called selection). I cannot get the syntax correct on my SQL query in MapInfo.

I believe I am looking to count the area of the object in the first layer that falls within the selection object but the syntax evades me.

Picture:


texts: Exporting XSD from GML using QGIS or ArcGIS for Desktop?

texts:
I want to export a XSD from a GML.
I already have the GML and I want to export a XSD that fits the GML

For an ETL tool we've programmed, people have to deliver a GML, XML and XSD. Those people already have a GML and XML, the ETL tool also requires an XSD. So I'm looking for an easy way to export an XSD that fits the GML

Is it possible to make an export in QGIS or ArcMap?


texts: How to make raster DEM from SHP with FME?

texts: I have a Point Featureclass in an ESRI FGDB with 22 Million points. The points have a distance of 10m from each other.I would like to create a raster DEM using SAFE FME 2014. Looking at the transformers, I can´t see a logical way to do this. Is there any standard methods to do this, or is FME really not suitable for this? I do have Spatial Analyst at my disposal but I wanted to try FME because it is usally 100 times quicker than ArcGIS.


texts: Autofilling a form field using a spatial query

texts:
I am trying to autofill values of fields in a form when adding features to a database. When the point is added there are several fields that relate to polygons that the point sit within, such as within a county shape or a postcode district.

I would like to be able to set these fields automatically with the name of the intersecting polygon by running a spatial filter with the layers that contain the polygons and the point that is being added.

I have tried to do this by using the Relative Relation option in the fields widget when adding the point, but this still requires a click on the canvas in the correct polygon to fill the box. There are 3 fields that would need to be updated this way and if several points are being added it is a lot of extra mouse clicks that it would be preferable to avoid. This method could be error prone if the user clicks the wrong polygon. Also, if I looked at the info for an existing point, these fields would be empty and would need to be set again.

I thought that I may be able to set it using Value Relation, but I cannot figure out how to make this work. I assume that I would need to the geometry of the shape that I would be testing against the point, which I would not know beforehand.

I tried to script something to add as a Python Init function (after reading this article) , but I could not get the code to work properly in the Python Console. I would get the value for every point within the layer and I cannot see how to identify the point that I am attempting to add.

Is there a better way of doing this? Or is there some code that I could use that gets the coordinates of the point that is being added to do the intersection test?


texts: Installing GeoNetwork on Debian with Tomcat?

texts:
How can I install GeoNetwork on a Debian server that already has tomcat7.

Could you give me help using the guide, please?

I followed the steps but I get the following error:

  Error                Raised exception while starting appl handler.
  Skipped. Exception        java.lang.IllegalArgumentException: Repository
  uses database URL of
  'jdbc:h2:~/.geonetwork/geonetwork;LOCK_TIMEOUT=20000;DB_CLOSE_ON_EXIT=FALSE;MVCC=TRUE'
  which does not match current database URL
  'jdbc:h2:~/.geonetwork/geonetwork;LOCK_TIMEOUT=20000;MVCC=TRUE'.
  Modify the svn property geonet:dburl on the root of the subversion
  repository at
  /var/lib/tomcat7/webapps/geonetwork/WEB-INF/data/data/metadata_subversion
  or specify a different subversion repository

  Message       Repository uses database URL of
  'jdbc:h2:~/.geonetwork/geonetwork;LOCK_TIMEOUT=20000;DB_CLOSE_ON_EXIT=FALSE;MVCC=TRUE'
  which does not match current database URL
  'jdbc:h2:~/.geonetwork/geonetwork;LOCK_TIMEOUT=20000;MVCC=TRUE'.
  Modify the svn property geonet:dburl on the root of the subversion
  repository at
  /var/lib/tomcat7/webapps/geonetwork/WEB-INF/data/data/metadata_subversion
  or specify a different subversion repository

  Handler       org.fao.geonet.Geonetwork

  Stack     java.lang.IllegalArgumentException: Repository uses database
  URL of
  'jdbc:h2:~/.geonetwork/geonetwork;LOCK_TIMEOUT=20000;DB_CLOSE_ON_EXIT=FALSE;MVCC=TRUE'
  which does not match current database URL
  'jdbc:h2:~/.geonetwork/geonetwork;LOCK_TIMEOUT=20000;MVCC=TRUE'.
  Modify the svn property geonet:dburl on the root of the subversion
  repository at
  /var/lib/tomcat7/webapps/geonetwork/WEB-INF/data/data/metadata_subversion
  or specify a different subversion repository at
  org.fao.geonet.kernel.SvnManager.(SvnManager.java:179) at
  org.fao.geonet.Geonetwork.start(Geonetwork.java:345) at
  jeeves.server.JeevesEngine.initAppHandler(JeevesEngine.java:600) at
  jeeves.server.JeevesEngine.init(JeevesEngine.java:206) at
  jeeves.server.sources.http.JeevesServlet.init(JeevesServlet.java:91)
  at javax.servlet.GenericServlet.init(GenericServlet.java:160) at
  org.apache.catalina.core.StandardWrapper.initServlet(StandardWrapper.java:1266)
  at
  org.apache.catalina.core.StandardWrapper.loadServlet(StandardWrapper.java:1185)
  at
  org.apache.catalina.core.StandardWrapper.load(StandardWrapper.java:1080)
  at
  org.apache.catalina.core.StandardContext.loadOnStartup(StandardContext.java:5026)
  at
  org.apache.catalina.core.StandardContext.startInternal(StandardContext.java:5313)
  at
  org.apache.catalina.util.LifecycleBase.start(LifecycleBase.java:150)
  at
  org.apache.catalina.core.ContainerBase.addChildInternal(ContainerBase.java:899)
  at
  org.apache.catalina.core.ContainerBase.addChild(ContainerBase.java:875)
  at
  org.apache.catalina.core.StandardHost.addChild(StandardHost.java:618)
  at
  org.apache.catalina.startup.HostConfig.deployWAR(HostConfig.java:963)
  at
  org.apache.catalina.startup.HostConfig$DeployWar.run(HostConfig.java:1600)
  at
  java.util.concurrent.Executors$RunnableAdapter.call(Executors.java:471)
  at java.util.concurrent.FutureTask$Sync.innerRun(FutureTask.java:334)
  at java.util.concurrent.FutureTask.run(FutureTask.java:166) at
  java.util.concurrent.ThreadPoolExecutor.runWorker(ThreadPoolExecutor.java:1146)
  at
  java.util.concurrent.ThreadPoolExecutor$Worker.run(ThreadPoolExecutor.java:615)
  at java.lang.Thread.run(Thread.java:701)


texts:
I know that this is not an uncommon question, but I've struggled to find an answer that works for my application.

I'm trying to find the number of nearby permits for any given property. In order to do this, I've used a SpatiaLite index to create a buffer around each point and check the other points for intersection with that buffer.

The problem is that for some reason--regardless of the CRS that the parent layer of points is set to--the buffer will not project in meters or feet (depending on the projection). As you can see, I've set the CRS when creating the layer (and checked that it was correct after layer creation).

filepath = "filepath/permits.csv"
uri = filepath+"?delimiter=%s&xField=%s&yField=%s&crs=%s" % (",", "Latitude", "Longitude", "EPSG:2264")
vlayer = QgsVectorLayer(uri, "Permits", "delimitedtext")

But the buffer still seems to be using decimal degrees for its projection. Any tips on how to get the buffer layer to project in feet or meters?

import sys

from qgis.core import *
from qgis.gui import *
from qgis.utils import iface

from PyQt4.QtGui import QApplication
from PyQt4.QtCore import QVariant

import itertools

# Initialize Application
app = QApplication([] , False)
QgsApplication.setPrefixPath("/Applications/QGIS.app/Contents/MacOS", True)
QgsApplication.initQgis()

# Create vector point layer from CSV
filepath = "filepath/permits.csv"
uri = filepath+"?delimiter=%s&xField=%s&yField=%s&crs=%s" % (",", "Latitude", "Longitude", "EPSG:3857")
vlayer = QgsVectorLayer(uri, "Permits", "delimitedtext")
print vlayer.crs().authid()

#Use loaded vlayer as active vlayer to set the CRS
QgsMapLayerRegistry.instance().addMapLayer(vlayer)

# ID of each feature
_NAME_FIELD = 'PermitNum'

# Names of the new field to be added
_NEW_NEIGHBORS_FIELD = 'WithinMile'

# Create 2 new fields in the vlayer that will hold the # of neighbors
vlayer.startEditing()
vlayer.dataProvider().addAttributes([QgsField(_NEW_NEIGHBORS_FIELD, QVariant.String)])
vlayer.updateFields()

# Create a dictionary of all features
feature_dict = {f.id(): f for f in vlayer.getFeatures()}
print vlayer.crs().authid()

# Build a spatial index
index = QgsSpatialIndex()
for f in feature_dict.values():
    index.insertFeature(f)

# Loop through all features and find features that touch each feature
for number, f in enumerate(feature_dict.values()):
    # Check to see if geometry is NoneType
    if f.geometry() is not None:
        print 'Working on %s' %  f[_NAME_FIELD]
        # Create buffer of given distance
        geom = f.geometry().buffer(5280,-1)
        intersecting_ids = index.intersects(geom.boundingBox())
        neighbors = len(intersecting_ids)
        print neighbors

vlayer.commitChanges()

print 'Processing complete.'


texts: Error using UpdateCursor to copy a GlobalID field

texts:
I am trying to use an Update Cursor to copy a GlobalID field into a string field.  However, I keep getting an error saying "Objects in this class cannot be updated outside an edit session [Medians]".

Aren't I initiating as edit session by using arcpy.da?

The error occurs on the fourth line.

workspace = "x/y/z"
edit=arcpy.da.Editor(workspace)
cursor = arcpy.da.UpdateCursor(medians, ("GlobalID", "GlobalID2",))
for row in cursor:
    globalID = row[0]
    globalID = globalID[1:-1]
    row[1]=globalID
    cursor.updateRow(row)


texts: QGIS 3.2 - How do I only display only one of my overlapping rasters at a time?

texts: I have a series of overlapping rasters and I only want one to display at a time so that the borders of neighbouring rasters are not obscuring anything on the centre raster. Is there a way to tell QGIS to only display the centre-most raster? Preferably this would only occur at a certain scale so one can still see the full extents of the raster when zoomed out.


texts:
QGIS 1.8.0 Mac OS 10.6.8: I want to label each polygon in a layer in a similar way to this image (done in WP Presentations) and am totally confused by the manual. Obviously I have to find a  list of polygons in the layer, then populate a field with the labels. Help getting started will be greatly appreciated.

http://johnsankey.ca/qgis1.jpg
http://johnsankey.ca/qgis0.png


texts: Has anyone had success in mapping WaterML data in cartodb? (WaterML). I can't get it to work and not sure what I am doing wrong. I also used this tool to generate a REST URL but still couldn't get it to work.


texts:
I have a code that wrote with C#.But I need to convert it to VBA and then use from it into GIS.
how can I convert it be VBA?


texts: retrieve rectangle from a latitude and longitude coordinates in FME

texts:
I wanted to retrieve a rectangle with FME from two points (point 1: latitude and longitude and point 2: latitude and longitude) that are stored in the postgres database in text fields:
text field 1 = latitude point 1
text field 2 = longitude point 1
text field 3 = latitude point 2
text field 4 = longitude point 2

  what transformer used on FME ???


texts: Averaging accumulated cost surfaces/shortest paths [r] [gdistance]

texts:
I'm trying to find a way to be able to average two accumulated cost surfaces so I can find the optimal route from both A to B and B to A. Although the difference is cost surface is due to the travel direction, the travel route is the same for the origin and goal point, therefore I need to average these.

As can be seen, the accumulated cost surfaces are different based on whether the travel is from A to B or B to A. And so, the optimal route calculated differs.

Is there a way to average the two accumulated cost surfaces to create a new cost surface that is based on both travels?

reproducible code:

library(gdistance)
r <- raster(system.file("external/maungawhau.grd", package="gdistance"))

altDiff <- function(x){x[2] - x[1]}
hd <- transition(r, altDiff, 8, symm=FALSE)

slope <- geoCorrection(hd)

adj <- adjacent(r, cells=1:ncell(r), pairs=TRUE, directions=8)
speed <- slope
speed[adj] <- 6 * exp(-3.5 * abs(slope[adj] + 0.05))

Conductance <- geoCorrection(speed)

A <- c(2667670, 6479000)
B <- c(2667800, 6479400)
AtoB <- shortestPath(Conductance, A, B, output="SpatialLines")
BtoA <- shortestPath(Conductance, B, A, output="SpatialLines")

AtoBacc <- accCost(Conductance, A)
BtoAacc <- accCost(Conductance, B)

par(mfrow=c(1,3))
plot(AtoBacc, main="A to B")
plot(BtoAacc, main="B to A")

plot(r, xlab="x coordinate (m)", ylab="y coordinate (m)", legend.lab="Altitude (masl)", main="Both")
lines(AtoB, col="red")
lines(BtoA, col="blue")
text(A[1] - 10, A[2] - 10, "A")
text(B[1] + 10, B[2] + 10, "B")

I understand that the accumulated cost surface is a raster layer, and so I could just do a raster calculation. However, I would need to then convert it back into a transition layer, which would change the values.


texts: Cause of ERROR 000732: Input Layer: Dataset plan does not exist or is not supported?

texts:
I have several mxds (located in C:\Project ) that have 250 layers each. All the layers as the same name- "plan" and they all in group layers like that (added picture for example):

I try to change all symbology layers of "plan" to the same symbology of layer "landuse":

All "plan" layers located one folder that divided to a lot of sub folders and sub-sub folders and their names (in the folders) are "kavim_kchulim".
i try this code:

import arcpy,os,sys,string
import arcpy.mapping
from arcpy import env

env.workspace = r"C:\Project"
symbologyLayer = r"C:\Project\layers\landuse.lyr"
for mxdname in arcpy.ListFiles("*.mxd"):
    print mxdname
    mxd = arcpy.mapping.MapDocument(r"C:\Project\\" + mxdname)
    dfList = arcpy.mapping.ListDataFrames(mxd, "*")
    for df in dfList:
        print df.name
        for lyr in arcpy.mapping.ListLayers(mxd, "", df):
            if lyr.name == "plan":
                arcpy.ApplySymbologyFromLayer_management("plan", symbologyLayer)
                print 'changed'
    mxd.save()
del mxd

but get an error:

>>>
Project.mxd
Layers

Traceback (most recent call last):
  File "C:\Users\yaron.KAYAMOT\Desktop\python.py", line 16, in <module>
    arcpy.ApplySymbologyFromLayer_management("plan", symbologyLayer)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.3\ArcPy\arcpy\management.py", line 6470, in ApplySymbologyFromLayer
raise e
ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000732: Input Layer: Dataset plan does not exist or is not supported
Failed to execute (ApplySymbologyFromLayer).


texts:
I am trying to retrieve data from 324 eopatches I sampled beforehand. Patches are basically numpy arrays as explained in this documentation EOPatch.
This is an example of the content of a patch after I resampled it :

EOPatch(
  data: {
    FEATURES: FeatureIO(/data/FEATURES.npy)
    FEATURES_SAMPLED: FeatureIO(/data/FEATURES_SAMPLED.npy)
  }
  mask: {}
  scalar: {}
  label: {}
  vector: {}
  data_timeless: {}
  mask_timeless: {
    LULC: FeatureIO(/mask_timeless/LULC.npy)
    LULC_ERODED: FeatureIO(/mask_timeless/LULC_ERODED.npy)
    LULC_ERODED_SAMPLED: FeatureIO(/mask_timeless/LULC_ERODED_SAMPLED.npy)
  }
  scalar_timeless: {}
  label_timeless: {}
  vector_timeless: {}
  meta_info: {}
  bbox: BBox(((500077.09501641133, 5095402.981379905), (501746.98615037295, 5097085.523204274)), crs=CRS('32633'))
  timestamp: [datetime.datetime(2017, 1, 1, 0, 0), ..., datetime.datetime(2017, 12, 19, 0, 0)], length=23
)

To retrieve certain data from the patch as my training and testing data. I first loaded the 324 patches like this:

# load sampled eopatches
eopatches = []
path_out_sampled = './eopatches_sampled_Slovenia'

for idx in range(len(patchIDs)):
    eopatches.append(EOPatch.load(f'{path_out_sampled}/eopatch_{idx}', lazy_loading=True))

eopatches = np.array(eopatches)

And then I used :

labels_train = np.array([eopatch.mask_timeless['LULC_ERODED_SAMPLED'] for eopatch in eopatches[patchIDs]])
features_train = np.array([eopatch.data['FEATURES_SAMPLED'] for eopatch in eopatches[patchIDs]])

However, it was taking forever executing and it was taking up so much RAM.

patchIDs is a list of numbers ranging from 0 to 323.

I tried to use this instead:

def data_retrieval(eopatch):
  feature_data=[]
  label_data=[]
  for eopatch in eopatches:
    f=eopatch.data['FEATURES_SAMPLED']
    feature_data.append(f)
    l=eopatch.mask_timeless['LULC_ERODED_SAMPLED']
    label_data.append(l)
  feature_data=np.array(feature_data)
  label_data=np.array(label_data)
  return feature_data, label_data

I just ended up with the same result.

I have also reinitialized the environment multiple times but nothing changed and I am using GPU.
Using h5py didn't help as well. I was able to load only 20 patches out of the 324.


texts:
In PostgreSQL, I want to use ST_Within from PostGIS to attribute a polygon reference to points.

Upon polygon modifications, I catch OLD.geom to erase the points feature in a before trigger. (Then, I catch NEW.geom to rename the new points after modifications.)

BEFORE UPDATE :

IF TG_OP = 'UPDATE'
THEN

    UPDATE pointstable
    SET polyid = NULL
    WHERE pointstable.polyid = OLD.polyid
      AND ST_Within(pointstable.geom, OLD.geom);

AFTER UPDATE :

IF TG_OP =  = 'UPDATE'
    THEN
    UPDATE pointstable
    SET polyid = NEW.polyid
    WHERE ST_Within(pointstable.geom, NEW.geom);

When it runs, it seems to catch the centroid of the modified points and sometime yield the wrong results. Some points polyid gets erased and not reassigned.

My main question therefore is : How can I iterate over the values in NEW or OLD in plpgsql. Must it be declared then defined (I am more or less new on the matter).

My second question is : Is there anyone familiar with the behavior of NEW and OLD, the SQL language and ST_Within that can explain to me what is happening exactly from the resulting query. As said, sometime it works, others half the points are renamed and the rest is not.

Edit :

I am repeating almost the same modification over and over. It works 90% of the time but once in a while, the points don't get renamed, any one has an idea?

I have a weird case where half the points gets renamed but the rest does not.
It seems to make no sense since the trigger is on polygons modifications. Therefore all points should be renamed or nothing. Because of that, I suspect it would have to do with ST_Within.


texts:
NZ Farmer gave a solution how to use a control function for a customised attribute form from within your own plugin. It seems that the control function accompanying a customised UI Form needs to be a py file somewhere located (addressed with layer.setEditFormInit), but I would like it to be part of my running plugin.

NZ Farmers'solution:
In the init method of your plugin:
Firstly add the plugin directory to the system path list (or a subdirectory of it) Secondly, initialise a static variable with the Plugin's instance (self)

However, I cannot figure out NZ Farmer's solution, i.e. how to put that solution in Python code and where the code should be placed. Or perhaps there is a different solution. Any help?


texts: In a software application displaying latitude & longitude, should N/S/E/W be localised?

texts:
In English software, it's conventional to render WGS 84 co-ordinates in a manner similar to the following:

  38.8951° N, 77.0364° W

or, encoding the directionality into the numeric values:

  38.8951°, -77.0364°

If I'm using the former, should I localise "N", "S", "E" & "W"? Or does WGS 84 (or any other widely adopted standard) allow for the universal use of these English abbreviations? What is the convention here?

Wikipedia's localised articles do alter these suffixes (e.g. "90°G a 90°D" on the Welsh-language "Geographic coordinate system" article), but I can't tell whether that's an isolated example or a manifestation of common global practice.

The same question of course applies for other representations such as 38° 53' 42.36" N.


texts: QGIS how to troubleshoot few SVG symbols that are displayed as a question mark?

texts: I created SVG symbols for points and pathed them correctly. When I export the project to test locally some points show up as a question mark symbol. I made sure the paths were correctly assigned and still have this issue. I selected the question marks and found that its a random few that have this issue. Has anyone experienced this before?


texts: Boundary of output raster for kernel density tool not congruent w shapefile boundary

texts:
I am using the Kernel Density tool in ArcGIS Desktop 10.4.1 to visualize population density around two metro area transit lines. According to ESRI: "The Kernel Density tool uses an algorithm to calculate density of point features around each neighborhood. The algorithm determines the default search radius (bandwidth), which allows for better weighting of highly dense points and smoother outputs. The following ArcGIS Help document, How Kernel Density works, further describes the tool." - source

The metro areas are Portland, OR and Charlotte, NC.

The tool input is a census tract shapefile joined with a population data table.

The resulting raster has the same boundary as a tract shapefile for one metro but is just a rectangle for other.

Any ideas what could cause this?

Kernel Density Tool results:

Portland, OR metro census tracts

Charlotte, NC metro census tracts


texts: Using vector affine transformation in QGIS?

texts:
I'm trying to use the QGIS plugin vector affine transformation to translate a shapefile with wrong coordinates using the right coordinates of another shapefile.
It looks like this

I read other answers to this problem but still my transformation does not work.
I take the coordinates of 3 points of the original shapefile:

396379.157,4992456.237
396434.636,4992428.498
396494.740,4992398.908

I take the coordinates of the same 3 points of the target shapefile:

396373.609,4992445.141
396431.169,4992416.939
396492.312,4992386.888

Then from online solver I get parameters to include into the transformation matrix:

U = -0.276
V = -2.628
W = 13625340.309
X = 19.103
Y = -197.695
Z = 380203686.883

or using another solver:

q = -0.279
r = -2.634
s = 13657434.35
t = 48.227
u = 158.434
v = 177246147.891

then I complete the QGIS transformation matrix

 -0.279 -2.634 13657434.35
48.227 158.434 177246147.891

But what happens is not a translation and it looks like this

What am I doing wrong?

Is it because my transformation does not include rotation?

Ideas about how can I solve this problem?


texts: QGIS 2.2.0.1 Valmiera, Windows 7 64 bit - copying out of GDALInfo-Window in QGIS not possible

texts:
When trying to copy information given out by GDALInfo with different kinds of raster data (DEM, normal georeferenced images etc.) through the interface of QGIS 2.2  by right-clicking and choosing copy (I want it in the clipboard) a Python error message keeps appearing.
The steps to get this error are as follows: I go over Raster - Miscellaneous - Info in the menu bar. Then I right click the information that's shown in the resulting window and choose "copy all" or "copy".
It's basically the equivalent of the QGIS GUI to entering gdalinfo and the location of the raster in the OSGEO4W shell.
Only solution is to either enter all data manually in the program i want it in, or to use the OSGEO4W shell and copy the info about the raster from there (since it's not possible to do that out of properties of raster data in QGIS). This is the error message that shows up every time I try, was like this in QGIS 2.0, and it hasn't changed in 2.2:

Fehler bei der Ausführung von Python-Code:

Traceback (most recent call last):
File "C:/PROGRA~2/QGISVA~1/apps/qgis/./python/plugins\GdalTools\tools\doInfo.py",
line 73, in doCopyLine
output.append( r.text() + "\n" )
AttributeError: 'str' object has no attribute 'append'

Python-Version:
2.7.4 (default, Apr  6 2013, 19:54:46) [MSC v.1500 32 bit (Intel)]

QGIS-Version:
2.2.0-Valmiera Valmiera, c3a2817

Python-Pfad: ['C:/PROGRA~2/QGISVA~1/apps/qgis/./python/plugins\\processing',
'C:/PROGRA~2/QGISVA~1/apps/qgis/./python', u'C:/Users/Lars/.qgis2/python', u'C:/Users
/Lars/.qgis2/python/plugins', 'C:/PROGRA~2/QGISVA~1/apps/qgis/./python/plugins',
'C:\\PROGRA~2\\QGISVA~1\\bin\\python27.zip', 'C:\\PROGRA~2\\QGISVA~1\\apps\\Python27
\\DLLs', 'C:\\PROGRA~2\\QGISVA~1\\apps\\Python27\\lib', 'C:\\PROGRA~2\\QGISVA~1\\apps
\\Python27\\lib\\plat-win', 'C:\\PROGRA~2\\QGISVA~1\\apps\\Python27\\lib\\lib-tk',
'C:\\PROGRA~2\\QGISVA~1\\bin', 'C:\\PROGRA~2\\QGISVA~1\\apps\\Python27', 'C:\\PROGRA~2
\\QGISVA~1\\apps\\Python27\\lib\\site-packages', 'C:\\PROGRA~2\\QGISVA~1\\apps
\\Python27\\lib\\site-packages\\PIL', 'C:\\PROGRA~2\\QGISVA~1\\apps\\Python27
\\lib\\site-packages\\win32', 'C:\\PROGRA~2\\QGISVA~1\\apps\\Python27\\lib\\site-
packages\\win32\\lib', 'C:\\PROGRA~2\\QGISVA~1\\apps\\Python27\\lib\\site-packages
\\Pythonwin', 'C:\\PROGRA~2\\QGISVA~1\\apps\\Python27\\lib\\site-packages\\Shapely-
1.2.18-py2.7-win32.egg', 'C:\\PROGRA~2\\QGISVA~1\\apps\\Python27\\lib\\site-packages
\\wx-2.8-msw-unicode', 'C:\\PROGRA~2\\QGISVA~1\\apps\\qgis\\python\\plugins\\fTools
\\tools']


texts: Format tilelayer.wms URL properties such that I add an overlay iteratively to my base map in Leaflet

texts:
My base map was created using hard-coded template variables. I am using the leaflet tilelayer.wms constructor and including the URL like so:

url: "http://localhost:8080/geoserver/sf/wmsservice=WMS&version=1.1.0&request=GetMap&layers=sf:archsites&styles=&bbox=589851.4376666048,4914490.882968263,608346.4603107043,4926501.8980334345&width=512&height=332&srs=EPSG:26713"

So I'm grabbing a WMS layer from GeoServer and I'm trying to make it display in my Leaflet base map, but I don't think I have to actually go about creating a feature group. I think the URL could possibly be simply:

http://localhost:8080/geoserver/sf/wmsservice=WMS

but then that would mean I need to structure the layer property which I am also unsure how to do so correctly. I think essentially I want to satisfy the GetCapabilities, GetMap, and GetFeatureInfo requests using Leaflet's API. And I think this is possible to do by utilizing only the layer group methods and the tilelayer methods. In addition I don't want to have to utilize onclick to retrieve the GetFeatureInfo because essentially I am adding in data programmatically.

I've seen examples utilizing the onClick eventlistener to do this, but I'd like it to be a toggle under the layer controls for visibility as it still is an overlay. So essentially I'd like to have a controls based listener as opposed to onclick.


texts: ESRI Shapefile not loading into SQL Server

texts:
I have been loading a number of shapefiles into SQL server 2012 using ogr2ogr.exe and the MSSQLSpatial driver. Up until this point, I have been very successful with doing so.

However, I have now run into difficulty with certain files and the errors thrown by various tools are not terribly helpful in identifying what exactly is wrong with the shapefile.

The files in questions are between 270MB and 290MB big and most have a single layer with more than 80,000 features in that layer. However, size and feature count is not itself an issue, I believe, as a file that is part of this data set is 280MB big and has 111,000(ish) features in it's particular layer.

The problem files all load correctly and are viewable in QGIS, including any files/features thrown out by ogr2ogr.exe as causing an issue. They also load into MapInfo successfully and can be edited there just fine. Re-exporting the layer from QGIS as a shapefile does not fix the issue with loading via ogr2ogr.exe

I've also attempted to load it using Shape2SQL and that program just hangs outright after it reaches a certain point with no error messages thrown at all.

The specific message thrown by ogr2ogr.exe is:

  ERROR 1: INSERT command for new feature failed. [Microsoft][ODBC SQL Server Driver][DBNETLIB]ConnectionWrite (send()).

  ERROR 1: Unable to write feature < feature no > from layer < layer name >.

  ERROR 1: Terminating translation prematurely after failed
  translation of layer < layer name > (use -skipfailures to skip errors)

The command line statement for ogr2ogr.exe I'm using is:

  ogr2ogr.exe -overwrite -f MSSQLSpatial "MSSQL:server=...;database=...;trusted_connection=yes" "< filename >.shp" -t_srs "EPSG:4326" -lco "GEOM_TYPE=geography" -lco "SCHEMA=ETL" -lco "OVERWRITE=YES" -progress

The shapefiles are in WGS84 projection.

Does anyone know what tool I might be able to use from the GDAL toolkit, QGIS or MapInfo to fix the errors in these shapefiles so that they can be loaded by ogr2ogr.exe into SQL Server 2012? Or alternately is it my ogr2ogr.exe arguments?

Update:

I've managed to run the Vertices Counter plugin and the first feature that causes an issue has 5,558,760 vertices.


texts: How to include a tolerance when using QGIS topology checker

texts: I'm using QGIS 3.14.3. The topology checker in this version doesn't seem to include a tolerance. I'm trying to relate points to lines by 'covered by' but keep getting errors.


texts: Display different layer based on zoom level with OpenLayers 3

texts:
I am new to OpenLayers and I was making a demo application which displays administrative areas based on zoom level. I am getting an error Cannot read property 'register' of undefined in map.events.register("zoomend", this, function(){});.

Here is my code:

var admin1 = new ol.layer.Tile({
title: 'Country Map with administrative level 1',
source: new ol.source.TileWMS({
  url: 'http://localhost:8080/geoserver/abc/wms',
  params: {
    LAYERS:'abc:abc',
    TILED:true,
    transparent: true,
    format:'image/png'
  }
})
});
var view = new ol.View({
   projection: 'EPSG:4326',
   center: [85.0583449970001, 28.006449365],
   zoom: 6
});
var map = new ol.Map({
    target: 'map',
    layers: [admin1],
    view: view
  });
Window.map = map;

map.events.register("zoomend", this, function(){
  var zLevel = map.getZoom();

    if( zLevel == 13 || zLevel == 14 || zLevel == 15  )
    {
        admin1.setVisibility(false);
    }
});


texts: pyproj, proj4 doesn't return same coordinates after inverse conversion

texts:
I'm trying to interpolate between coordinates, as a step in that process I convert latlong to UTM so I can work in metric.

I don't get the expected latlong when converting back, reproduced below:

>>> import pyproj
>>> pyproj.__version__
'1.9.2'
>>> p = pyproj.Proj("+proj=utm +zone=33 +north +ellps=WGS84 +datum=WGS84 +units=m +no_defs")
>>> p(58.4, 15.6)
(5579373.890337878, 2328887.997220567)
>>> p(*p(58.4, 15.6), inverse=True) # Should return (58.4, 15.6)
(58.284316956873305, 15.63501250163002)

Is there a good explanation for this? How do I do what I want without this issue?


texts: Should geopackage work in QGIS on Linux?

texts:
I'd been converting some shapefiles to GeoPackage with QGIS 2.18.5 on Windows. Now I check my project on Linux (Ubuntu 14.04 with QGIS 2.18.5), and the gpkg files are treated as attribute-only layers.

To test, I made a new project and added some sample vector data from http://www.geopackage.org/#sampledata (the Esri Sample Data, GeoPackage 1.0 and 1.1), with the same result. No geometry is recognized (and the gpkg_* metadata tables are instead also displayed in the wizard):

Any clues what is wrong? I couldn't find any similar reports, or extra requirements for QGIS on Linux.


texts: Create multiple "administrative" geographies in Tableau and filter

texts:
We have shapefiles for:

a) voting districts
b) Wards
c) Neighborhoods
D) Police Districts

which all overlap in non-uniform ways. We also have other data around population, income, etc that we would like to visualize.

How do we create one map in which each shapefiles above is a single layer (in Tableau, ideally). It sounds like this should be obvious, but I can not find the solution I'm looking for.

For example: One layer would be voting district, a 2nd would be ward, and a third would be neighborhood. We could then look at income in the same viz across all three and see if/how income changes as those administrative groupings change.

Tableau version: 10.5.

Suggestions?


texts: How to calculate tolerance for TopologyPreservingSimplifier in JTS?

texts:
I am trying to simplify some areal geometries (POLYGON\MULTIPOLYGON) using JTS's

TopologyPreservingSimplifier.simplify(Geometry geom, double distanceTolerance)

From the JavaDoc:

 void   setDistanceTolerance(double distanceTolerance)
          Sets the distance tolerance for the simplification.

I cannot understand what should be specified for the distanceTolerance paramter, I would like to have a tolerance of 0.1 Km (or miles).

Do I need to convert the input to decimal degrees ? If so, what conversion can that be used accurately for most of the globe ?


texts: How can I add the Spatial Analyst ArcToolBox?

texts:
I tried to add the Spatial Analyst ArcToolBox by right click and selecting 'Add Toolbox' but I cannot find in any folder the available toolboxes. They seem that they have disappered. I work with ArcGIS 9.3.
Any help please?


texts: Need help creating Europe reference grid map using EEA shape and spatialite files

texts:
I'm trying to generate a reference grid map of Europe with country and provincial boundaries using the EEA files available here. Here is an example of country and provincial boundaries with a raster layer overlay.

I just want the underlying shape map with boundaries, no colors/data/raster stuff overlaid. The linked EEA website has files for all of the individual countries, but doesn't seem to have a single folder for Europe as a whole. Do I have to download every single country to make the Europe map? The eea_reference_grid_v1.pdf at the site seems to suggest that there is a single folder you can download for the whole Europe map, but I can't find it anywhere. Additionally, any pointers on loading the files and producing the map in R are very welcome.


texts:
I have a bunch of layers and corresponding SLD's running a WMS service. I was exploring the possibility of shifting to a TileMill's tile based TMS service.

For this, I would like if it were possible to convert the SLD to CartoCSS, and then start the styling from that point on in tilemill.

Is there a SLD to Carto CSS?

I have seen the question Are there tools to directly import QGIS projects into Tilemill?, but the answers there really do not answer my question.

It's almost 2016, and Mapbox has renamed their old product as Mapbox Studio Classic.  However, I'm still looking for a way to convert SLDs to CSS, so that we can convert old SLDs to CSS for conversion projects.


texts: shapefile driver table name mismatch for ogr vs. shp2pgsql

texts:
I have a shapefile which if I convert to sql statements using shp2pgsql
 return table name all in lowercase, for e.g. 'liberia'.

When I read the same shapefile using ogr and use layer.GetName() it returns
capitalized name e.g. 'Liberia'.

Is the  mismatch because I am not looking at the right variables?


texts: QGIS: Unable to project raster layers correctly on a Toshiba Tecra 64# Win7

texts: I know how to project a raster, I can do it on other machines. On this particular machine the raster loads but are no where near the correct position as if the datum is on another planet.


texts: GeoServer / GeoWebCache serve precached TMS

texts:
I am pretty new to GeoServer and GeoWebCache

I have a set of precached tiles (TMS) that I want to server via GeoServer / GeoWebCache.

I have saved tiles in C:\Program Files (x86)\GeoServer 2.13.0\data_dir\gwc

and want to serve them in the following form (as per the TMS structure of the cached data):

http://localhost:8080/geoserver/gwc/service/tms/1.0.0/layername/z/x/y.formatExtension

How do I specify the 'layername' configuration, I think it is in the geowebcache.xml

If so, does anyone have a thorough example that I can follow?


texts: What's the difference between US Census DMA and MMA?

texts:
DMA and MMA are commonly used to denote a large, metropolitan municipal entity.

Is there a difference between the two? Is there a certified definition to these terms? Google is a bit inconclusive about the issue.


texts:
There is an annoying green background feature displayed in the Georgia, USA OpenStreetMap data that doesn't appear in other adjacent states. See:

http://www.openstreetmap.org/relation/161957#map=8/33.246/-83.790

The green background is rather overwhelming and obscures other critical map data for my needs. Has anyone encountered this and found a solution to turning off the "green layer"?


texts: Calculating distance from a central point in three dimensions and visualising the result

texts:
We will shortly be collecting data from our remote systems that track ADS-B signals from Aircraft.
This data will provide the lat/lon in epsg:4326 (wgs84) and a height in meters.
The station location is known in EPSG:4326 as a lat/lon.

We want to be able to visualize the effective range of the signals for planning additional stations. The 3D is required to visualize LOS (line of Sight) shadows due to mountains, man made obstacles, etc.

How do we do the following using PostGis:

Calculate the bearing from True North from the station to the
aircraft position
Calculate the distance over the ground between the two points
Calculate the distance to the point factoring in the height as well.

Assuming that we then store in a PostGis Database the furthermost distance for each radial, how do we do the following using PostGis:

Create a cool looking polygon that connects the points with curves in
2D for use in OpenLayers
Create a cool looking polygon that connects the points with curves
taking in the height (z) for output to a kml file for use in Google
Earth to visualize the range three dimensionally


texts: QGis 2.18 Actions window does not display 'Actions Properties'

texts:
I have installed QGis 2.18.19 64bits in Windows 7. My problem is when I open the layer properties actions window, I do not see the menu Action Properties.
Apparently it should be there according to the QGis documentation (e.g. QGis 2.18 training manual). What I am missing.

I also noticed that the button is named "Create default actions" while in the guide is named "Add default actions"

Here is my screen capture:

Here is the capture of the training manual pdf image:


texts: How to get user's geolocation with ESRI API for JS 3.16?

texts:
I want to add some button to find user's current geolocation. Can I do it with Esri API for JavaScript 3.16? I want something like Search widget.
Also I have seen Locate widget, but it is from 4.0.


texts: How to batch a model created with "model builder" in arcmap 10?

texts:
I created a simple model that selects iteratively elements from a point dataset and calculates a cost distance raster for each of these points.  It works fine when I run it for just one dataset, but when I try to batch it get the following error.

ERROR 000865: Input raster or feature source data: I_atlanta_OBJECTID does not exist.
ERROR 000581: Invalid parameters.
Failed to execute (Cost Distance).

Please any help will be greatly appreciate.

Here a link to a snapshot of the model
http://img88.imageshack.us/img88/7432/modelmultidist.png


texts:
Geoserver SLD:

  <sld:LineSymbolizer>
    <sld:Stroke>
      <sld:CssParameter name="stroke">
        <ogc:PropertyName>geo_color</ogc:PropertyName>
      </sld:CssParameter>
      <sld:CssParameter name="stroke-width">2</sld:CssParameter>
    </sld:Stroke>
  </sld:LineSymbolizer>
  <TextSymbolizer>
    <Label>
      <ogc:Literal>&#x227b;</ogc:Literal>
    </Label>
    <Font>
      <CssParameter name="font-family">Lucida Sans</CssParameter>
      <CssParameter name="font-size">20</CssParameter>
    </Font>
    <LabelPlacement>
    <LinePlacement>
    <PerpendicularOffset>0</PerpendicularOffset>
    </LinePlacement>
    </LabelPlacement>
    <Fill>
      <CssParameter name="fill">
        <ogc:PropertyName>geo_color</ogc:PropertyName>
      </CssParameter>
    </Fill>
    <VendorOption name="followLine">true</VendorOption>
    <VendorOption name="forceLeftToRight">false</VendorOption>
    <VendorOption name="conflictResolution">false</VendorOption>
  </TextSymbolizer>

I have the above style set for a WMS layer containing line strings. For most zoom levels the direction arrow aligns properly with its corresponding line. But on zooming out beyond a point, the arrows seem to be displaced away from their respective lines. Shown here:

I am fairly sure that this has to do something with VendorOption parameters but unable to put a finger on it. Appreciate any inputs to get this sorted out.


texts:
I have about 1000+ running routes which I can load into memory as polylines using latitude and longitude data captured from a running watch. From these polylines I would like to be able to identify routes that match a segment of that polyline as shown on the map (black).

Any advice on a method or mathematical method to be able to give me a confidence score that the marked segment is contained within the bigger segment. I ask for a confidence score, as there will always be GPS drift so I want to match segments that have say 95% confidence that the two lines are the same say buffered to 10 meters.

I am going to implement this using C# when I work out how but just need a point in the right direction.


texts: Getting lines across Legend symbol in ArcMap export to PDF?

texts:
I am encountering an issue with an item in the legend when I use a page size of 36x48 for exported PDFs. The item, a library symbol, appears to have lines going through it. It looks fine with smaller page sizes. I have tried changing colors and enlarging the symbol, as well as opening up a new mxd and adding just the library data, but the problem persists.

Any solutions?


texts: High quality public transport trip planner examples? Trip visualization examples?

texts:
Do you know high-quality online public transport route/trip planners? I'm looking for sites with recognized usability, most carrier-owned planners are bad or very bad from Human-Computer Interaction point of view.

Also: do you know any papers or works on trip planner GUI design patterns? Example: instead of expandable list with search results we can use the time graph with bars like in Nokia Transit - I just want to see more of such stuff. Thanks!


texts:
How to zoom in map in OpenLayerss 3 on click of button?

I want to zoom in and zoom out the map on custom buttons instead of + and - buttons which are on default map.

I tried with:

onclick="map.zoomIn();"

But it's not working, though it was working in OpenLayers 2.


texts: I have published a GIS service with a relation to a table view (in ArcMap I created a "relate" between the point shapefile and the table and the related records all show up here - then I did "share as service"). When I display it in the web map (viewer for flex) you are able to scroll through the correct number of related records, but no data is displayed for each - it says only "no related records".


texts: Projecting data without spatial reference in ArcGIS Desktop?

texts:
I have a raster file downloaded from here.  When I load it into ArcMap, it says the data is missing spatial reference, and that it can be drawn in ArcMap but not projected.  Nevertheless, using the Define Projection tool, I can assign a coordinate system that seems to work pretty well.

How does ArcGIS Desktop do this if there was no spatial reference to begin with?

Is this a valid approach to data missing spatial reference?


texts: Seeking bare streets image of specific neighborhood?

texts:
For an infographic project, I would need the bare layout of the streets of a specific neighborhood. What I mean is just the streets "lines", without any name, color or symbol, really just the streets, on a transparent background if possible.

Is this something QGIS could do for me?


texts: Is there an alternative for JTS TestBuilder on Linux?

texts: JTS Testbuilder, an application for testing JTS functions, is available on Windows and I use it quite regularly. Are there anything like this on Linux? Does GEOS have anything similar to this on linux?


texts: Getting QGIS Server GetCapabilities working on Windows?

texts:
I followed this tutorial: http://anitagraser.com/2012/04/06/qgis-server-on-windows7-step-by-step/ but nothing happens.

And when I try to open it this error message appears

So I tried to follow httpd_qgis.conf edition step from http://hub.qgis.org/projects/quantum-gis/wiki/QGIS_Server_Tutorial but it gives me the same error.

I think I didn't need that step because I can see "mod_fcgid" in my apache modules directory

What am I doing wrong?


texts: Making file path part variable using Python

texts:
I am creating a python script to help automate updates.
I am just wondering does anyone know how to make part of my file path a variable, as each month the python script will run via a different file path (only the date will change)
For example
Current workspace is
env.workspace = r"X:\Support\Data\Shapes\Dec20\SA"

Next month it will be
env.workspace = r"X:\Support\Data\Shapes\Jan21\SA"

Is there a way i can make the "Date" a variable?
I.e
DATE_IS = Jan21

env.workspace = r"X:\Support\Data\Shapes\"DATE_IS"\SA"


texts:
I am experiencing a strange behavior when symbolizing a point layer using any kind of Feature Blending.

I am currently using QGIS 2.18.16 but I have seen this in previous versions (edit: tried now in QGIS 3.0.0 and the problem is still there). What happens is when some points are in the same exact coordinate, the feature blending shows the bounding box of the circle representing the point.

Is there any way to avoid this? Has anyone detected this behavior?
I guess this is a general bug but maybe is something related to our video card.

Edit: Stranger things even stranger. If I set all symbol colors to transparent, I can get a map of just the bounding boxes of the overlapping points.

Taking a closer look at the screen, the bounding box appears always but you can notice it better when the points overlap since it is almost transparent.

The composer manager has a different behavior. It doesn't show the bounding boxes and when exporting the map the bounding box is not there.
We have tried this in three computers and the result is the same.


texts: Making heat map with joined layer using QGIS

texts:
I need to make a heat map with the number of cases by neighborhood.

I have:

a polygon shapefile with the neighborhoods, and
a CSV file containing the name of the neighborhood and the number of
cases.

Is there a way to do this in QGIS 3.4 LTS?


texts: CHM (Canopy height model) = DSM (Digital surface model) - DTM (Digital terrain model)

texts:
I want to create a canopy height model (CHM) which has the height of objects (trees, building, etc.) from DSM and DTM. Below is the step by step procedures which I did so far.

I use open source s2p software to generate DSM from stereo pleiades images. According to the offical webpage of s2p, the output DSM's height is ellipsoid height. And the spatial resolution of stereo images and DSM is 50 cm.
I use python elevation package to download DTM (the website only said it's DEM but not specifically said it's DSM or DTM, so I assume that it's DTM). According to the website of the elevation package, the DTMs are obtained from SRTM which in theory should be geoid. And the spatial resolution of DTM is 30 m.
Now I have both DSM and DTM. I subtracted DSM by DTM and got CHM (I fixed the problem of two models having different spatial resolutions by resampling the DTM).
According to the graph below, I expect that the pixel value of the street height should be 0 or closed to 0.

However, the result is not what I expected.
The pixel value 36 is the height of a point on the center blue triangle which is the ground (not trees or building), which is not 0 or closed to 0.

Here are my questions:

Is my idea of creating CHM correct?
If 1. is correct, which step could be wrong? I think there might be two potentional problem. First, the accuracy of two DSM and DTM are not clear. However, I don't think they will influence the result too much. Second, the DSM and DTM are not on the same reference system.
If 2. is correct, how do I convert ellipsoid height to geoid height and vice versa in a given area? (I don't think it's even possible though.)


texts: Land use land cover change post- classification in ArcGIS for cross- tabulation?

texts:
I have finished the image classification process in another software package but many papares recommend the ArcGIS package to process post-classification which means from which one land scape type has been changed to what.

How can I process post-classification in ArcGIS 9.3 to look at the matrix of change or "From-to" change?


texts: Arcpy turns off all fields in a MXD file that all of the values are Null

texts:
I am trying to turn off all fields in a MXD that have all null values.

I have been looking at a SearchCursor without much luck.


texts: Looking for HowTos or Examples of Presentation Platforms targeting Management/Clients

texts:
Any good ideas about using Google Earth to present self contained geospatial data&app to a client or senior management?

The broader issue being explored here is that many of our clients are not likely to maintain GIS systems themselves, or be competent enough to install a free app and then browse the data. And we would like the legacy we leave to a client to be independent of us, so we don't want them to be relying on our server or any one else's server. A single file container (yes I know you will say SpatiaLite but what app?) that can be USB'd or the like.

One feature I really like in GE, is to geolocate photos and have them appear on the map, to provide a record, or some interactivity etc. I can do the same in QGIS, but not without knowing a lot about the installation. In GE, I know that if they have got the KMZ, they will find the photos.

And what about senior management? Because you know they are luddites at best, won't have time to view it when you have time etc.

So my thought here is to put GE up as the app, and KMZ as the container, and ask the question in this forum.

What is the best example of GE (or your choice app) as a self contained presentation?


texts:
I have a polygon and raster with the same extent:

library(rgdal)
library(raster)
library(rgeos)
library(sp)

 # make the polygons
df <- structure(list(LAT = c(35.61226, 35.34986, 35.17794, 34.60425,34.40356, 33.94433, 33.41834, 16.89566, 16.89561, 16.89561),
                     LON = c(-9.604802, -9.803048, -9.921968, -10.30782, -10.44971,-10.76656, -11.13053, -24.99776, -24.99788, -24.99773)),
                .Names = c("LAT","LON"), class = "data.frame", row.names = c(1L, 2L, 3L, 4L, 5L,6L, 7L, 8L, 9L, 10L))

dfSp <- SpatialPointsDataFrame(matrix(c(df$LON, df$LAT), nrow = nrow(df)), data = df)

grid <- makegrid(dfSp, n = 20)
gridSp <- SpatialPointsDataFrame(grid, data = data.frame(id = rownames(grid)))
gridSpRas <- rasterFromXYZ(gridSp)
gridPoly <- rasterToPolygons(gridSpRas, dissolve = T)

# make the raster
grid2 <- makegrid(dfSp, n = 80)
grid2Sp <- SpatialPointsDataFrame(grid2, data = data.frame(id = rownames(grid2)))
gridSpRas2 <- rasterFromXYZ(grid2Sp)
rvals <- sample(seq(1,8,1),ncell(gridSpRas2),replace=T)
r <- setValues(gridSpRas2, rvals)
extent(r) <- extent(gridPoly)

I now need to clip the polygons based on values in the raster. I want to remove any polygon that has a cell with a certain value within it. In this case, anytime there is a raster cell with value of 1 within a polygon I want to remove that polygon.

For a simple example such as this I can do:

values <- extract(r, gridPoly)
values2 <- lapply(values, table)

x <- do.call(rbind, values)
y <- cbind(gridPoly@data,x)
collist <- c("1","2","3","4")
sel <- apply(y[,collist],1,function(row) 1 %in% row)
z<-gridPoly[!sel,]

But the real raster has hundreds or thousands of cells per polygon. Is there a simpler way to subset a polygon by a raster? So going from this:

plot(gridPoly)

to this:

plot(z)


texts: Regridding a rasterbrick in R using bilinear interpolation


texts: Google Earth Enterprise Portable Server replacement?

texts:
Is anyone still using Google Earth Enterprise Portable Server? Is there a better alternative now? "The Google Earth Enterprise Portable Server solution enables organizations to access portions of their custom globes in the field when no network access may be available or possible." This is exactly what I am looking for.

All of the references to it that I found were at least 3 years old.
https://support.google.com/earthenterprise/answer/4485225?hl=en#google_earth_api

Reference from 2011:
Using Google Earth imagery data for portable solution?


texts: Reclassify function in ArcMap 10.5 does not correctly reclassify values?

texts:
I reclassified a raster layer containing three classes (0-0.05, color = white; 0.05-0.5, color = dark salmon; and 0.5-1, color = red)

to a raster layer where theses ranges were replaced with 0 (color = brown), 1 (color = purple) or 2 (color = blue), respectively. While the reclassify function correctly converted values for most raster cells, it does not do so for one (there is a class 2 cell missing in the top left).

I have 9 such raster files that I reclassified. For some, the reclassify function correctly reclassifies all cells, but for others it does not. Despite that the steps to reclassify are the same for all. Any ideas what causes this discrepancy, or what am I overlooking?


texts: Error in the connection to Oracle Spatial from QGIS

texts:
I am not being able to stablish the connection to an Oracle Spatial server in QGIS.

Apparently I have the connection configured accordingly with QGIS documentation, but it stills gives me an error:

  ORA-12514:TNS_listener does not currently know service requested in connect descriptor

Here are my connection parameters as set in QGIS:

I have used the same parameters that I have set in the connection settings in SQL Developer (that works fine):

I have tried to use tnsping command to check the connection parameters and it is successful, regardless use:

tnsping [server_ip]:[port_number]/[SID]

or

tnsping [tnsname]

where [tnsname] is defined in the tnsnames.ora file with the same connection parameters indicated above:

  [tnsname] = (DESCRIPTION = (ADDRESS_LIST = (ADDRESS = (PROTOCOL = TCP) (HOST = XXX.XXX.XX4.95) (PORT = 1521)) ) (CONNECT_DATA = (SID = anza1) (SERVER = DEDICATED) )  )

I am really running out of options... Can anyone help?

NOTE: I'am using QGIS 2.0.1-Dufour


texts: PyQGIS cannot set binCount in provider.initHistogram

texts:
I'm using initHistogram in QgsRasterDataProvider, however I have not been able to set binCount. I am using the following commands:
provider = rlayer.dataProvider() # Where rlayer is a raster layer
provider.initHistogram(QgsRasterHistogram(),1,100) # 1 is band # and 100 is binCount
hist = provider.histogram(1) # 1 is band #
histlist = hist.histogramVector

No mater what value I put in for binCount (100, 10000...), I get a histogram with a binCount of 1000, when I use QGIS3.12, and binCount of 10000000 when I use QGIS3.14.
Can you specify the number of bins you want in a histrogram using binCount in initHistrogram?
Why does the the binCount change depending on which version of QGIS I use?


texts:
I have 2 layers, point, and line. I want all the points to snapped to the nearest lines.

Effectively, at the end, all these points should now lie in the lines retaining all of their attributes. Is this possible?