texts: How to convert coordinates from Gauss-Krueger to Longitude/Latitude?

texts:
I'm a newbie in the GIS world so please advise me :)

In my .NET/C# project I need to show some shapes/points on the Bing Map. The problem is I have the Gauss-Krueger coordinates and need to convert them to Longitude/Latitude.

Can anyone point me to the right direction please? Is there maybe a free managed library for this purpose? A link with a good explanation of conversion algorithm would be helpful too.

Thanks!

Anil


texts: How to fix invalid polygon with self-intersection - Python

texts:
I have a self-intersecting polygon inside a shapefile. Does there exist a Python library which spots the self-intersections and makes it possible to correct them?

I have found it is possible in QGIS but how can I achieve the same result in Python? Do I need a PostGIS DB or I can handle it with a Python library?


texts: Annotating a cartodb torque timeline-control

texts:
I'd like to know what options there are to customise a cartodb torque timeline control.
I've seen awesome visualisations made with CartoDB/Torque for FIFA '14 matches by twitterdata ( for instance http://cartodb.com/v/worldcup/match/?TC=x&vis=31ed2e2a-06ce-11e4-8c1d-0e230854a1cb&h=t&t=Brazil,FFCC00%7CGermany,B40903&m=7%2F8%2F2014%2017:00:00%20GMT,7%2F8%2F2014%2018:52:00GMT&g=108%7C11,23,24,26,29,87,97#/2/-17.6/-5.6/0 ).

From trying to read that code I inferred that twitter data created a modified Torque to do just that.

Before going in that direction, what are the best ways to customise Torque Timeline controls?

What I specifically want to achieve is have a graph displayed above the timeline that 'guides' the user as to when events happen.
In this example: https://emileaben.cartodb.com/viz/6a6d789c-8430-11e4-a483-0e4fddd5de28/public_map , there is an event at around 7:32 CEST (lots of red dots in the US, representing Internet connectivity going down), and having a graph with number of events over time, on the timeline axis would really help guide users to what is going on.


texts: How to open layer using OpenLayers 2 and GeoServer?

texts:
I know it's a very basic question but I've been searching for an answer for over 4 hours and I can't get it to work.

I have GeoServer running locally and I'm trying to open a layer in my browser through OpenLayers 2.

This is the code I have:

   <title>OpenLayers Example</title>

I copied everything exactly like the example in the OpenLayers website, only changed the URL for the layer and the layer workspace/name.

I tried changing the url to "localhost:8080/geoserver/nurc/wms" and every other possible combination, and I can't get it to work. All I can see is the "+" and "-" for zooming in and out of the layer, but the layer is blank, and sometimes it shows only pink tiles.

Firebug does not show any errors.

Any suggestions please?


texts:
I have used the code from Mapnik python tutorial. Now I want to use zoom_to_box to zoom to certain extent. I want to zoom to Box2d(-180.0, 0.0, 0.0, 90.0), however, after appied zoom_to_box, I saw the map envelope is Box2d(-180.0,-45.0,0.0,135.0), so what do I miss here?

The testing shapefile can be downloaded here: shapefile

Here below my code:

import mapnik
m = mapnik.Map(256,256)
m.background = mapnik.Color('steelblue')
s = mapnik.Style()
r = mapnik.Rule()
polygon_symbolizer = mapnik.PolygonSymbolizer(mapnik.Color('#f2eff9'))
r.symbols.append(polygon_symbolizer)
line_symbolizer = mapnik.LineSymbolizer(mapnik.Color('rgb(50%,50%,50%)'),0.1)
r.symbols.append(line_symbolizer)
s.rules.append(r)
m.append_style('My Style',s)
ds = mapnik.Shapefile(file='ne_110m_admin_0_countries.shp')
print ds.envelope()
layer = mapnik.Layer('world')
layer.datasource = ds
layer.styles.append('My Style')
m.layers.append(layer)
print m.envelope(), layer.envelope()
m.zoom_all()
print m.envelope(), m.scale()
bbox = mapnik.Box2d(-180.0, 0.0, 0.0, 90.0)  # 00
m.zoom_to_box(bbox)
print m.envelope(), m.scale()
mapnik.render_to_file(m,'world00.png', 'png#')
print "rendered image to 'world.png'"


texts: New fields added via ArcGIS Online not appearing on Collector?

texts:
This question has appeared in some random threads (https://community.esri.com/thread/217461-field-added-on-arcgis-online-does-not-appear-on-collector) with no clear solution.

Basic premise:

Simple collector setup with one KML file (two small polygons, hosted feature layer)
Initially 8 fields that required answering
I have full permissions over the groups/settings
Maps can be used offline and are working properly

The problem:

I needed to add 3 more fields to the questionnaire, which I successfully added via ArcGIS online. However, these fields are not showing up on the collector application. Moreover, the fields are not selected as 'on' when viewing the online map, and ticking them 'on' does not persist when closing and re-opening the online map.

Attempted solutions:

Removed and re-downloaded collector survey

Is there a setting/step I am missing?


texts: Calculating percentage of number of features in layer in QGIS

texts:
I'm using QGIS to get some "stats" on a current layer e.g. to determine the number of features of a certain class. That I do by using the Aggregate-Algorithm on the class-attribute with count on "ID".  I get a table that shows the number of features per class.
But how could I get another attribute giving me the percentage of the total number of features?
While this is easily done by hand (the number of total features in the layer can be easily determined, then it's just count / total_number), I don't know how to access the total number of features of the layer via an expression.


texts:
I currently have two file folders with a mxd in each of the folder. These two file structures and mxd name will be kept constant. Once a week or month, I will update these two mxd. Then, I will generate a couple of pmf (with data package) from these mxd, and save them into a another folder (fixed location).

Is there a way I can generate pmf with data package automatically (update interval set to once a week / month)?

Im on 9.3.1 Desktop, Server license, with Publisher. Sorry to say this but I do not have knowledge of scripts/python.


texts: One-to-many joins on a feature class to a table

texts:
I am using ArcMap 10.2.2 with an Advanced license. I have used ArcMap for other basic analyses, but am fairly new to the program and have only basic coding skills, with no SQL experience outside of very basic queries using the SQL code builders in ArcMap.

I have polygons of species' geographic ranges and I have a table with a list of pathogens/diseases carried by those species.  Certain species carry multiple pathogens.  In an ideal world, I'd join the polygons of the ranges to the table of the species+pathogens and keep the polygon/spatial data, but it doesn't seem like there is a GIS tool to do this (I am going to later count overlapping pathogens, so it's important that each species-pathogen combination has its own feature, rather than one feature that has the information for multiple pathogens).  These are what my data look like (#1 & #2, #3 is for later):

Attribute table for spatial features:

Table without any spatial data:

(Please note that these are not real data, and the pathogen-species combinations may not be real.  I actually have 115 polygons and 519 table rows.  Species carry anywhere between 1-40 pathogens and thus species have anywhere between 1-40 rows in my table)

As an added bonus, because of future analyses I need to do, I would ideally like to be able to do this one-to-many join such that I could keep an extra column of data from the table (see below for an example of how this might look).  However, this is less important because if I can figure out how to do the join, I can just redo the analysis.

There are a few solutions that I found to this problem, but they didn't work for me for reasons listed below:

Solution 1: David Aalbers's blog - I am too much of a novice to figure out how to edit the newer script and figure out how to import it into ArcMap.  I tried the older script and it just generated a ton of empty polygons titled test, test_1, test_1_1, test_1_1_1, etc.

Solution 2: One-To-Many Join - I tried to use the Make Query Table tool (and I have no SQL knowledge, although I used the expression builder) and I got ERROR 000383: Issue with a table, cannot find the workspace & Failed to execute (MakeQueryTable).  All my files are in the same folder, but they are not in a geodatabase.


texts: Boolean precondition in ModelBuilder not exporting to Python script

texts:
I don't know python so I often use ModelBuilder and export to python script. In this case, I have a model (that in itself is made of python scripts that have been exported from models) that works fine where if the precondition (ISG Template and VC Template) is false, it does not execute that branch and if true, it does.

When I export it to python script, the boolean does not exist and even if no parameters are input, the full script runs.

Below is the python script that runs and produces an output regardless of whether the boolean is checked or not.  How can I get the boolean precondition to work in Python? I'm using ArcGIS Desktop 10.6.  All the parameters are optional in both the model and the python script.

# -*- coding: utf-8 -*-
# ---------------------------------------------------------------------------
# Parent.py
# Created on: 2019-07-14 16:21:30.00000
#   (generated by ArcGIS/ModelBuilder)
# Usage: Parent <ISG_Template> <Output_ISG_DWG_File> <VC_Template> <Output_VC_DWG_File>
# Description:
# ---------------------------------------------------------------------------

# Import arcpy module
import arcpy
import os

# Load required toolboxes
arcpy.ImportToolbox("//cadd/users/steve/CADGIS/Version2019/GIStoCAD.tbx")

# Script arguments
ISG_Template = arcpy.GetParameterAsText(0)

Output_ISG_DWG_File = arcpy.GetParameterAsText(1)

VC_Template = arcpy.GetParameterAsText(2)

Output_VC_DWG_File = arcpy.GetParameterAsText(3)

# Local variables:

# Process: VC
arcpy.VC_GIStoCAD(Output_VC_DWG_File)

# Process: ISG
arcpy.ISG_GIStoCAD(Output_ISG_DWG_File)

Furthermore, below is the script with my modifications.  I still doesn't work.  The script runs regardless of whether boolean is checked but does not produce an output.  It runs in about a second so its really not doing anything.  It should take a minute.

# -*- coding: utf-8 -*-
# ---------------------------------------------------------------------------
# Parent.py
# Created on: 2019-07-14 18:31:23.00000
#   (generated by ArcGIS/ModelBuilder)
# Usage: Parent <ISG_Template> <Output_ISG_DWG_File> <VC_Template> <Output_VC_DWG_File>
# Description:
# ---------------------------------------------------------------------------

# Import arcpy module
import arcpy

# Load required toolboxes
arcpy.ImportToolbox("//cadd/users/steve/CADGIS/Version2019/GIStoCAD.tbx")

# Script arguments
ISG_Template = arcpy.GetParameterAsText(0)

if ISG_Template == True:
    arcpy.ISG_GIStoCAD(Output_ISG_DWG_File)

Output_ISG_DWG_File = arcpy.GetParameterAsText(1)

VC_Template = arcpy.GetParameterAsText(2)

if VC_Template == True:
    arcpy.VC_GIStoCAD(Output_VC_DWG_File)

Output_VC_DWG_File = arcpy.GetParameterAsText(3)


texts: Handling null inside .map in Google Earth Engine

texts:
I want to map over an ImageCollection and calculate percentiles of a region. Some images are empty for my area of interest, due to cloud masking.

How do I handle something like:

if variable ==null :
  continue with next image

I prepared some code with the exact problem:

GEE code

Can the following normalize function be altered to incorporate the occasional null value?

var normalize = function(Im){
  var time = Im.get("system:time_start")
  var percentiles = Im.reduceRegion({
    reducer:ee.Reducer.percentile([10,90],["10","90"]),
    crs:proj.crs,
    crsTransform:proj.transform,
    maxPixels:1000000000
  })
  var p10 = ee.Number(percentiles.get("band_10"))
  print(p10,"p10")
  var p90 = ee.Number(percentiles.get("band_90"))
  print(p90,"p90")
  return Im.subtract(p10).divide(p90.subtract(p10)).set({"system:time_start":time})
}


texts: How to check if points are within a circular sector using PostGIS?

texts:
I want to get the stars within the black region as shown above. To find this I have the centre of circle (lat, long), radius and the angle of region (45 degree).

I am new to postgres. I worked on code to check if points exist within the circle of a given radius.

How I can find points within this black region?


texts: What is the difference between CQL_FILTER and FILTER?

texts:
In a CSW request what is the difference between the parameters CQL_FILTER and FILTER ?

Examples:

localhost:8080/geonetwork/srv/eng/csw?REQUEST=GetRecords&CONSTRAINT_LANGUAGE_VERSION=1.1.0&RESULTTYPE=results&VERSION=2.0.2&TYPENAMES=gmd%3AMD_Metadata&CONSTRAINTLANGUAGE=CQL_TEXT&SERVICE=CSW&CONSTRAINT=Product

and

localhost:8080/geonetwork/srv/eng/csw?REQUEST=GetRecords&CONSTRAINT_LANGUAGE_VERSION=1.1.0&RESULTTYPE=results&VERSION=2.0.2&TYPENAMES=gmd%3AMD_Metadata&CONSTRAINTLANGUAGE=FILTER&SERVICE=CSW&CONSTRAINT=Product

Is there any CQL documentation?


texts: medium resolution Earth coastal outline in tabular form

texts:
I'm generating my own images in Python and Blender. What I need is data that I can read directly from a table/file to draw coastlines with 1 to 10 kilometer resolution. I'm not familiar with high level Geographic or Cartographic software, I just need dots to connect (lat/lon or xyz) or something similar.

While the package Basemap exists for python, I'd like to avoid installing it in each of my python versions. I just need the points.


texts: Designing circle label with many attributes in ArcGIS for Desktop?

texts:
Without using ArcPy, how can I create a circle label show many pieces of information from an attribute table in ArcGIS?

I have image about planning zone map, with specific label on top of objects, now I need digitizing map in ArcGIS.


texts:
I am trying to create a user, say dev, for cartodb running in local machine. When i run the script create_dev_user

sh script/create_dev_user dev

It produces this RuntimeError: ERR unknown command 'HMSET' error. I am on ubuntu 10.04.


texts: After upgrading to Qgis 2.14 Essen, cannot load layers from .zip

texts: I had a project in 2.12 Lyon that used many .zip files for access to .shp files.  Now, after upgrading, I cannot add any content from .zip.  This is Win10 Machine and I have the option checked to scan contents of .zip (FULL SCAN)  I cannot access my old zip files for my layers, nor can I access any new downloaded content.


texts: Building a Google Maps Directions URL with waypoints (that aren't stops)

texts:
I want to go from a place to another using a Google Maps Directions URL in the Google Maps Mobile App. This is from Lisboa to Porto, in Portugal :)

https://www.google.com/maps/dir/38.759741,-9.153736/41.1454953,-8.6108788/

As you can see, I only gave the origin and destination coordinates and it works.

In my final solution I would like to ask the Google Maps App for the same route, but forcing it to pass through a location (39.657293,-8.827500). This is not a stop so I I cannot treat it as a waypoint.

I saw several options in other posts but none is working for me.

Any thoughts?

Than you very much for you help


texts: Is there any way I can use a Key-Value store for geospatial data?

texts:
I have used many relational databases in the past, but I have also read about all NoSQL databases, and the Key-Value stores looks interresting.

When I store geometric object I mostly use five indexed columns ID, MIN_X, MAX_X, MIN_Y and MAX_Y (where X and Y are in a map projection). I don't need an index on my other data.

I need the X and Y values to lookup objects in a specified place (map rectangle), and I need the ID value if I want to update an specified object.

Is there any way that I can use a Key-Value store for this?


texts:
What does the GEOS LineString "normalize()" routine really do?

GEOS has virtual void geos::geom::LineString::normalize() with the documentation:

  A normalized linestring has the first point which is not equal to it's reflected point less than the reflected point.

Ignoring the "it's" typo, I still can't figure out what this really means.


texts: Using inline variable %Value% from iterate feature selection to specify feature class name for clip tool

texts:
What I am trying to achieve: I have a dataset of ~3000 health facilities in 52 districts. I want to draw Thiessen polygons around the health facilities which conform (are contained within) to the district boundaries. I am able to select all the health facilities in each district district using iterate feature selection and draw Thiessen polygons around them. I have also created a  feature class file for each district to use as the clipping feature for each set of thiessen polygons. The name of these files is identical to the output of the %Value% inline variable produced by the iterative function.

I would really like to automate the clipping of the thiessen polygons. I have tried various methods:

Using the inline variable value as a child in the select data tool (does not work as value is not a member of gdb)
using the inline variable to specify the name of the clipping feature in the clipping tool e.g E:\PhD\Objective1\Objective1.gdb\%Value%

Any ideas?


texts: What happens to pending tasks when GeoServer reboots?

texts:
I am running a GeoServer/GeoWebCache system on a Windows server. We do regular updates which require retiling and I have an automated system (using Task Scheduler) that sends requests to the GeoWebCache REST API. There are times when there can be thousands of requests generated and the retiling can take several hours or more.

My question is: What happens with the pending / currently executing tasks if GeoServer (or the main server) is rebooted?


texts: Replacing masked cloud areas with another satellite image?

texts:
I have various Landsat images of which many are cloudy. Using the cloud masking plugin available in QGIS, I was able to identify and remove the areas with cloud and cloud shadows from every image.

With this, I would like to replace the masked areas of one image with the data available from the others. How should I go about this?


texts: How to find namespace for WFS request 1.1.0?

texts:
I am trying to do a GetFeature request to this  WFS Server (GetCapabilities), but unable to find the argument for TYPENAME and thus  get this error:

msWFSGetFeature(): WFS server error. Incomplete WFS request: TYPENAME parameter missing

A DescribeFeatureType request also works, but attempting to do the GetFeature request with the <element name="Baustellen_aktuell> does not seem to work in this format:

http://geoweb1.digistadtdo.de/OWSServiceProxy/baustellen.ows?SERVICE=WFS&VERSION=1.1.0&REQUEST=GetFeature&TYPENAME=ms:Baustellen_aktuell&MAXFEATURES=5&OUTPUTFORMAT=text/javascript

How can I find the appropriate namespace for the TYPENAME argument from the XML responses of the GetCapabilities or DescribeFeatureType requests?


texts:
I want to convert a shapefile from polygon to point:

pfile = pathtofile
file = gpd.read_file(pfile) #it is a polygon

Searched this site but didn't find anything similar:
http://geopandas.org/geometric_manipulations.html


texts: How to achieve parallel Kriging in R to speed up the process?

texts:
I have tried the gstat package in R to perform bathymetry predictions based on Kriging interpolation methods. However, the process is very time consuming due the large amount of sampling data I'm using. Also, the whole process in R is running only in one core of the four cores my PC have.

I wonder if it is possible to parallelize the Kriging (krige) function of the gstat package or
Any R package to perform parallel Kriging


texts: Setting up a section using shapefile.Reader

texts:
I would like to open some shapefile with shapefile.Reader but I would like to see some part of it, for example 10000-20000 records. How to set up a section? I can't find that in documentation.
with shapefile.Reader(shapefile_path) as shp:


texts: Fixing python error in QGIS Freehand raster georeferencer plugin?

texts:
I am using QGIS 2.8.9 stable release in the 32bit form on a Windows 7 64 bit compatible PC. I want to use the Freehand raster georeferencer plugin. It installs without difficulty but when I try to open a raster in the plugin I get the error code 'An error has occurred while executing python code' and this readout:

raceback (most recent call last):
  File "C:/Users/hpethen/.qgis2/python/plugins\FreehandRasterGeoreferencer\freehandrastergeoreferencer.py", line 190, in addLayer
    self.createFreehandRasterGeoreferencerLayer()
  File "C:/Users/hpethen/.qgis2/python/plugins\FreehandRasterGeoreferencer\freehandrastergeoreferencer.py", line 201, in createFreehandRasterGeoreferencerLayer
    layer = FreehandRasterGeoreferencerLayer(self, imagePath , imageName, screenExtent)
  File "C:/Users/hpethen/.qgis2/python/plugins\FreehandRasterGeoreferencer\freehandrastergeoreferencer_layer.py", line 46, in init
    self.initializeLayer(screenExtent)
  File "C:/Users/hpethen/.qgis2/python/plugins\FreehandRasterGeoreferencer\freehandrastergeoreferencer_layer.py", line 157, in initializeLayer
    self.resetScale(sw, sh)
  File "C:/Users/hpethen/.qgis2/python/plugins\FreehandRasterGeoreferencer\freehandrastergeoreferencer_layer.py", line 166, in resetScale
    wratio = sw / iw
ZeroDivisionError: float division by zero

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

QGIS version:
2.8.9-Wien Wien, 26c3ece

Python path: ['C:/PROGRA~2/QGISWI~1/apps/qgis-ltr/./python/plugins\processing', 'C:/PROGRA~2/QGISWI~1/apps/qgis-ltr/./python', u'C:/Users/hpethen/.qgis2/python', u'C:/Users/hpethen/.qgis2/python/plugins', 'C:/PROGRA~2/QGISWI~1/apps/qgis-ltr/./python/plugins', 'C:\PROGRA~2\QGISWI~1\bin\python27.zip', 'C:\PROGRA~2\QGISWI~1\apps\Python27\DLLs', 'C:\PROGRA~2\QGISWI~1\apps\Python27\lib', 'C:\PROGRA~2\QGISWI~1\apps\Python27\lib\plat-win', 'C:\PROGRA~2\QGISWI~1\apps\Python27\lib\lib-tk', 'C:\PROGRA~2\QGISWI~1\bin', 'C:\PROGRA~2\QGISWI~1\apps\Python27', 'C:\PROGRA~2\QGISWI~1\apps\Python27\lib\site-packages', 'C:\PROGRA~2\QGISWI~1\apps\Python27\lib\site-packages\PIL', 'C:\PROGRA~2\QGISWI~1\apps\Python27\lib\site-packages\jinja2-2.7.2-py2.7.egg', 'C:\PROGRA~2\QGISWI~1\apps\Python27\lib\site-packages\markupsafe-0.23-py2.7-win32.egg', 'C:\PROGRA~2\QGISWI~1\apps\Python27\lib\site-packages\python_dateutil-2.2-py2.7.egg', 'C:\PROGRA~2\QGISWI~1\apps\Python27\lib\site-packages\pytz-2014.2-py2.7.egg', 'C:\PROGRA~2\QGISWI~1\apps\Python27\lib\site-packages\win32', 'C:\PROGRA~2\QGISWI~1\apps\Python27\lib\site-packages\win32\lib', 'C:\PROGRA~2\QGISWI~1\apps\Python27\lib\site-packages\Pythonwin', 'C:\PROGRA~2\QGISWI~1\apps\Python27\lib\site-packages\Shapely-1.2.18-py2.7-win32.egg', 'C:\PROGRA~2\QGISWI~1\apps\Python27\lib\site-packages\wx-2.8-msw-unicode', u'C:/Users/hpethen/.qgis2//python', 'C:\PROGRA~2\QGISWI~1\apps\qgis-ltr\python\plugins\fTools\tools', 'C:/Users/hpethen/Documents/Docs/Work/Teaching/EES_FieldSchool']

I have had previous releases of QGIS on here. Is it possible that there is some conflict? Is this a 32bit vs 64 bit issue or is there something else I'm missing.


texts: Hosting HTML files so they are accessible through AGOL

texts:
I recently made a web mapping application using the ArcGIS for JavaScript API, but I am having a difficult time making it accessible to users of my organizations AGOL account.

I naïvely thought I could easily host HTML files on my organizations ArcGIS Server through either ArcGIS Server Manager or ArcGIS for Developers, but this was unsuccessful.

I am currently hosting the files on GitHub and I am capable of setting up a web server with either Python or PHP, but I am confused as to how to eventually make everything available in AGOL.

I am starting to get the impression that this may not be possible and I should be hosting the HTML files/app on a non Arc web server and just send people the URL for the app instead of having them access the app via AGOL? Any guidance would be helpful here.


texts: How to divide my area into adjacent polygons containing 80 points each?

texts:
I have a couple of thousand points in my dataset. I need to create clusters that each consists of 80 points. The specific selection of which points to put into a cluster is not critical - only the proximity. From these 80-point clusters I need to create 1 polygon pr. cluster, so I end up with dividing my area into adjacent polygons containing 80 polygons each.

A solution usable in QGIS and other open source tools would be great!


texts: Why does R(gdal) mangle data when converting to geojson?

texts:
I've encountered issues with projections in .geojson files created by rgdal, but I never thought it would destroy my data, until I tried this (reproducible example):

downloader::download("http://census.edina.ac.uk/ukborders/easy_download/prebuilt/shape/England_lad_2011_gen_clipped.tar.gz", destfile = "cuas.zip")
untar("cuas.zip")
cuas_shp <- raster::shapefile("England_lad_2011_gen_clipped.shp")
tmap::qtm(cuas_shp)
rgdal::writeOGR(cuas_shp, "cuas-ogr", "OGRGeoJSON", "GeoJSON")
cuas_ogr <- rgdal::readOGR("cuas-ogr", "OGRGeoJSON")
tmap::qtm(cuas_ogr ) # fails - something wrong with the data
library(sp) # plot using normal methods
plot(cuas_shp)
plot(cuas_ogr)

There are more examples I can show. It seems the geojsonio package does not solve these issues either. Ask if you'd like more reproducible examples - want to check I'm not missing anything obvious first though.

(Also not the projection seems to be missing in the geojson example.)


texts: which v.clean option is required to eliminate overlapping segments of lines?

texts:
I want to produce a polygon layer out of lines like in this image with the command poligonize:

Source: QGIS create boundaries from background image

The problem is that some of the lines are coming from polygons with the command polygon to line tool in QGIS like in this image:

Source: How to merge overlapping and duplicate lines from vector layer to one single line?

I am struggling to understand which options I should use in v.clean to remove the overlapping segments of the line layer, or which other stepts or tools I could use to get the desired output in a straight way.

The main asumption is that lines match exactly.

My software QGIS 2.18 and GRASS 7.4


texts:
I have tables where the primary ID field is an integer that starts at 2 000 000. This number is allocated as the next available number - entered manually or calculated from the last available number. The progress of the count (above 2 000 000) is spread across multiple tables so they are not necessarily in sequence in any particular table.

In QGIS v3.0, when I look at that ID in the table view, the integer is truncated and presented the way Excel does when the number is too large. So 2048273 shows as 2.04827+06. I can double click in the field to see the actual ID, but that is a nuisance. I want to be able to see the ID without having to click.

How can I do that? Or is this similar to the other question about the American format and was fixed in a later version?


texts: Seeking workflow to build cloud optimized tif?

texts:
By its definition, a cloud optimized GeoTIFF is a regular GeoTIFF file, aimed at being hosted on a HTTP file server, whose internal organization is friendly for consumption by clients issuing.

Based on the gdal doc, I can create a cloud optimized tif by running

gdal_translate in.tif out.tif -co TILED=YES -co COPY_SRC_OVERVIEWS=YES -co COMPRESS=LZW

on a file with existing overviews.

Because of the size of my tif files, I would like to avoid making my tif file "cloud optimized" after building the pyramid in order not to duplicate the writing.

Is it possible to use the -co COPY_SRC_OVERVIEWS=YES when I build my pyramids with gdaladdo or is it necessary to run it in an additional step ?


texts: Documentation on "Analyzer errors" codes?

texts:
As noted here, when you receive an error of the form

001272 : Analyzer errors were encountered (codes = <value>).

this means that the attempt to stage and publish a service failed.

The tool provides a set of codes in the error message, but the documentation I linked to only tells you to fire up an ArcGIS Desktop tool to find out the problem. Instead, I would prefer to use the code to look up what the error is; this would be especially helpful to me since the publish is happening automatically in Python via arcpy.

Is there a resource that tells us what the codes mean?

I am using 10.2 at the moment and will soon be upgrading to 10.3.


texts: Is there a program to read all the data contained in a waypoint in a .gpx track file and assign it to a variable?

texts:
This is related to geodata.  I am a plant explorer and go into Ecuador and Peru looking for new species of plants.  I take close to 1000 pictures each time I go, and I use a Garmin Oregon GPS device to create a track each day.   I can use GeoSetter to read all the GPS Coordinate data and the Elevation and put it into the appropriate pictures which match the Track timestamp.  My Garmin also records ambient temperature data in the track file, but I know of no program that will read it and put it in the appropriate field of the EXIF photos.

Is anyone aware of a program that will do this?

If not, is there a program that will read each waypoint in a track file and assign the various values to a list variables, so that I can then write a program to match them against the Creation Time of each of the photos in my directory?

I can then use ExifTool to write the temperature data back into the EXIF data.


texts: Finding which polygons lie along a polyline in ArcMap 10.4

texts:
I work in ArcMap 10.4.

I'm looking for a tool/function which would return a list of polygons through which a polyline goes. To be more specific - I have a polygon layer of administrative units and polyline features (gaspipes etc). I made a simple sketch of what I mean (see pic below).
Is there a way to return the list of "touched" polygons the polyline "goes" through?


texts: Efficient intersect test with OpenLayers

texts:
In an OpenLayers map I have data of several hundred/thousand buildings. There are also like ~100 other objects on the map which can move. I need to check each 100 ms if a object can see another with a maximum viewrange of say 1000 m.

For the maxviewrange I implemented an OK solution (checking distance only if src or dst moved).

Now I need to check if a line from src to dst intersects with any building. I would create a boundingbox using src and dst, then query all buildings in the boundingbox and check for an intersection with the line src dst.

I think that this is a pretty suboptimal implementation. Maybe there are some fancy implementations for the problem... But I just could not find one in OpenLayers.

Has anyone an idea to improve that or a link to a function I just did not notice?


texts: Marker Line Symbol with Character Marker Symbols gets grayed when published to map service

texts:
After publishing as Map Service character marker symbols become grayed even if their colour is black in definition. This odd behaviour happens in ArcMap 10.2.2 and ArcMap 10.3.

Are there are any issues registered in ESRI KB regarding this problem? I could not find a similar one.

EDIT: Tried to play with different symbols and it seems that text colour is black only for "thicker" characters or fully filled ones. For example Arial letters turns into gray but Arial Black (thicker font) letters are black. This could be anti-aliasing related issue.

EDIT2: Played with Text Anti-Aliasing options when Publishing Parameters dialogue. Gray with options Force and Normal, black only with None. But other texts look ugly with option None... What workaround is possible here?..

The line symbol definition is shown below, I don't see any magic here.


texts: Performing per-pixel linear regression over multiple rasters with NoData?

texts:
I have a stack of 10 Landsat satellite images over the same area that I have calculated NDVI for. I want to calculate NDVI trend over time on a per-pixel basis. This means that for each pixel in the images I want a value for the slope of the regression NDVI vs. Time.

From a previous question asked here I came across the Curve Fit ArcGIS add on. I input my NDVI rasters into the tool, and asked for an output being linear regression slope and p-value. I got those results, but only for a small subset of my study area. I found that any pixel set as NoData (aka areas I had masked due to cloud cover) led to that entire pixel stack being out put as NoData. So even if 1 out of 10 values in a pixel stack was NoData the linear regression would fail and the output would be NoData.

My goal is to find a way to run linear regressions on all the pixels in my study area, even if one or more of the values in a pixel stack are masked due to cloud cover (meaning if 1 value is set to NoData then the linear regression will simply be out of 9 values out of 10 in a stack, rather than just being output as NoData).


texts: Unit of $length in QGIS Field Calculator

texts:
I'm calculating field of various features in my shapefile using QGIS.

I'm getting some sample length as 0.134914579766986.

I want to know in what unit (eg: meters, km etc)  is the above value ?


texts: Converting CAD Annotation to Shapefile with attributes using FME?

texts:
What transformer should I use to add AutoCAD Mtext to GIS attribute table using FME?

The reader file in CAD and I want to transform it to GIS Shapefile. In this case I want to expose the attribute from text (fme_text_string) to DC_LAND_PARCEL_BOUNDARY(polygon geometry).
So how to expose the attribute along with the value, so that in ArcGIS i only turn on the labeling options to show the annotation.


texts: Min and max zoom levels/scales for a QGIS project

texts:
We've built a qgis project for use in education and we want to set up a min. and max. zoom level (or scale) for the whole project.

In the project properties, we can specify some "project predefined scales", but this does not prevent from zooming over of these limits.

On each layer, we can set a min/max scale for rendering the layer (in Layer Properties > Rendering), but again, this only hides the layer when zooming over of these limits.

Is this feature available? Can we make it using some pyQGIS chunks of code?


texts:
I read in a raster and turned it into a binary image (land/water mask) of 0's (water) and 255's (land). I've named this array src_landCoverArray.
I want to output this array as a shapefile with water as polygons and land as empty space.
I've tried using rasterio.features.shapes() described here (https://stackoverflow.com/questions/37898113/python-boolean-array-to-polygon) but it produced a shapefile with the wrong extent. My code:
import shapely.geometry
myShapes = rasterio.features.shapes(src_landCoverArray)
polygons = [shapely.geometry.Polygon(shape[0]["coordinates"][0]) for shape in myShapes if shape[1] == 0]
crs = {'init': 'epsg:4326'}
my_gdf = gpd.GeoDataFrame(crs=crs, geometry=polygons)
my_gdf.to_file("Water_Poly2.shp", driver='ESRI Shapefile')

I've tried using a similar method with rasterio described here (How to polygonize raster to shapely polygons) but failed again with the location/extent being incorrect. My code:
mypoly = []

for vec in rasterio.features.shapes(src_landCoverArray):
    mypoly.append(shape(vec[0]))
crs = {'init': 'epsg:4326'}
my_gdf = gpd.GeoDataFrame(crs=crs, geometry=mypoly)
my_gdf.to_file("Water_Poly_mypoly.shp", driver='ESRI Shapefile')

I feel that I'm close to success if I could figure out the issue of extent/location being wrong. Has anyone tackled a problem similar to this?


texts:
Plotting arrows in a map to show shifts in coordinates or points using ggplot2
I have three points in my map (latitudes and longitudes of the points: a, b and c are provided below). I would like to draw arrows of different colors and later show the points and arrows on the legend of my map. The arrow should move from point a to point b (this arrow can be red in color and named first). Another arrow can be drawn from point b to point c (this arrow can be green in color and named second). The points in the map can be of different shapes ie rectangle, triangle etc Point a can be named current, point b named future_1 and point c named future_2
below is my code
load require libraries
library(maptools)#
library(rgdal)
library(rgeos)
library(sp)
library(ggplot2)

get the map of the world
data("wrld_simpl")
crs(wrld_simpl)
plot(wrld_simpl)
#plot africa
africa <- wrld_simpl[wrld_simpl$REGION==2,]
plot(africa)
##change into a format that can be used in ggplot2
AG <- fortify(africa)
crs(africa)
names(AG)
##m points
a <  -1.28333, 36.81667
b <- 3.327709,  6.457171
c <- -17.82772, 31.05337
##plot using ggplot2
ggplot()+ geom_polygon(data=AG, aes(long, lat, group = group), colour = alpha("darkred"), size = 0.8, fill = 'skyblue', alpha = .3)


texts: Seeing whether set of postcodes have typical distribution to whole postcode population of UK?

texts:
For my dissertation I have a dataset of postcodes of 9984 that I am converting to coordinates atm.

I looking to find whether this dataset has a 'typical or 'atypical' spatial distribution compared to the whole of the UK Postcodes spatial distribution.


texts: Selecting points from QGIS when button is greyed out?

texts:
I am new to QGIS.

I added mailing addresses from arcgis map server url https://gis.nccde.org/agsserver/rest/services/BaseMaps/Addresses/MapServer. I want to select some points here and move them to a new layer to work with fewer points. But the "select features by area or single click" button is greyed out.

In arcgis pro, I would use "select" tool to select some points. Then use "copy features" tool to move it to a new layer.

How do I get the select tool to be enabled in QGIS?


texts: Batchtransform multiple shapefiles with FME

texts:
I would like to batch-transform multiple shapefiles which are stores in several different folders.

I used the batch deploy with the option "recreate source folder structure" which works well for creating the folder structure.

What's messed up are the attributes of the shapefiles. The "feature types to read" seem to be bound to the shapefile I used to insert the shapefile-reader.

Is there something I can do to transform the shapefiles with the correct attributes (and without additional FME-specific attributes)?

( I used a workbench with shapefile reader-->csmap-reprojector --> shapefile writer but tested also a version with a creator, feature-reader, csmapreprojector and feature-writer).


texts: Prevent Overlapping Features in QGIS from Being Faded

texts:
I have multiple layers in QGIS that I've set as transparent so I can see their respective features. However, I would like one layer (showing proportionally sized points) to NOT fade when it overlaps with a block group layer.
In other words I'd like to move my points on top and have UNfaded color fills, (aka "Move to Front" as one does with MS products), but otherwise retain the visibility of the underlying block group layer everywhere else.


texts:
I want to display a map in offline mode (without Internet connectivity), so I have imported OSM data into my PostgreSQL database using osm2pgsql installer.
I went through various links and tutorials and now I am currently using this link:
https://sharpmap.codeplex.com/wikipage?title=Add%20a%20PostGIS%20layer

When I run this code using SharpMap, it gives the error of "geometry column not found" in the table.
I am stuck as i don't know what should I enter as "tablename" and "idColumn" in the above reference link code and particularly WHY?

Can anyone guide me?


texts: How to resolve the conflict between two layer labels?

texts:
Despite the fact that all labels placement methods and styles was tried ArcMap still fails to resolve the conflict between two layers labels?

In my case, I have two layers (the built-up area) and (well-known sites). And I am using maplex engine I tried all placement and position changing procedures, but unfortunately the labels still have clashes between them as shown in the screenshots below:


texts: Using points in polygons tool with QGIS and CRS?

texts:
I am working with 2 layers.

1 shows a polygon of a district and is in CRS EPSG 4326 and another layer (address layer) is in CRS EPSG 3082.

I right click the address layer and "save as" and specify the same crs as the district layer (epsg 4326).

When performing the "points in polygon function" shouldnt the result be a layer that when I open the attributes table it should show the number of addresses in the "pnt cnt" field?

Is there something I am missing about setting the CRS for all layers to be the same in order to perform Geoprocessing?

Is there another issue that I am missing?


texts: How to disable loading fonts from fast.fonts.net?

texts:
I have installed ArcGIS JS API (3.19) locally. But fonts were loaded via Internet.

If I have no Internet connection loading of my app is very slow.

What can I do?
Is it possible to disable loading via Internet?


texts: Adding multiple plugins to custom menu in QGIS

texts:
I am new to Python and QGIS and attempting to create some plugins to replicate in-house tools used in other GIS programs. I have added a custom menu using Python and added the plugin to the menu using the following code within the initGui method:

self.menu = QMenu("&Menu Name", self.iface.mainWindow().menuBar())
actions = self.iface.mainWindow().menuBar().actions()
lastAction = actions[-1]
self.iface.mainWindow().menuBar().insertMenu(lastAction,self.menu)
self.action = QAction(QIcon(":/plugins/Trial/icon.png"),"Trial plugin", self.iface.mainWindow())
self.action.triggered.connect(self.run)
self.menu.addAction(self.action)

However, if I try and add another plugin to the same menu I end up with duplicated menu item names with only one plugin in the menu, rather than all the plugins under the same menu. Is there a way of adding multiple plugins to an existing custom menu item?


texts:
CREATE VIEW new_line AS
SELECT
point1.gid,
ST_MakeLine(point1.geom, point2.geom)
FROM point1, point2
WHERE point1.gid = point2.gid;

It's possible to load the layer but QGIS says 'There is no entry in geometry_columns!'. The SRID is missing, too.


texts: How to design class breaks for a time series of choropleth maps?

texts:
I'm working with a time series of choropleth maps showing the same indicator for the same administrative boundaries, but spanning period of nine years.

What would be the most appropriate strategy of designing class breaks for such a time series?

At the moment, I think the quintiles division would serve the purpose the best, since it allows to scan quickly for 20% of areas with highest/lowest values (for instance).

That works for one map, however I''m not sure how to porperly choose a scheme for time series.

Options that I came up with so far would be:

use class intervals created separately for each map - advantage being they reflect the data structure for each year; disatvantage - harder to compare the maps between years.
use class intervals of quintiles from the summary map for the whole time period (without yearly split of the underlying data) - advantage being easier comparison of the maps for each year; disadvantage - breaks might not be ideal in certain years.
use class intervals of quintiles taken from the whole distribution of the data - all areas+all years - advantages disadvantages will most likely be simmilar to the above scenario.

What are your thoughts?


texts: Why does QGIS algorithm 'Join attributes table' change the types of the shapefile attributes?

texts:
In QGIS 2.18.9...

I have a shapefile with two attributes 'ID' and 'POP', each of which is of the integer type.

I have a CSV file loaded as a delimited, text-only table layer with two attributes 'ID' and 'EMP', each of which is of the integer type.

I use the standard QGIS algorithm 'Join attributes table' to join the table to the shapefile based on the common attribute 'ID'.

The join executes as desired, except that:

1) I can only get an output shapefile if I choose to create a temporary output layer instead of specifying a file name and location for it to save to

2) in this temporary output shapefile, the two attributes from the original shapefile are now of the text/string type instead of the integer type even though the two attributes from the joined table are still of the integer type


texts:
I have a set of points(lat, lon and other attributes) stored as an Asset in Google Earth Engine. I import this asset as a Feature Collection and visualize the points over Satellite Imagery in Google Earth Engine. I would like to generate a heat map from these points.

I searched extensively over the web, but could not find any resource elaborating on doing this in Google Earth Engine. I did come across this documentation https://developers.google.com/maps/documentation/javascript/examples/layer-heatmap of the Google Maps JS API that uses |google.maps.visualization.HeatmapLayer|.

Is there an equivalent of this in Google Earth Engine?


texts: GeoServer GetMap with geopackage very slow (with high zoom levels)

texts:
I have a GeoServer instance and have imported multiple high res tiff files as layers and created a layergroup from the layers. When I call the WMS endpoint with a specified bbox and the geopackage format, it can take multiple hours to get the geopackage from the server. It appears when my gridset includes zoom levels 19 or above is when I begin to see the time increase exponentially. A jpeg or PNG returns quickly.

Are there optimizations that can be made either to the query or to how we are organizing our layers that could help? Even if the bbox provided is within a single layer, the query against the layer group is extremely slow.

Example query: http://localhost:8080/geoserver/wms?service=WMS&version=1.1.1&request=GetMap&layers=TestWKSP:TestLGroup&bbox=627883.09857,8299208.89693,628180.15347,8299506.25009&width=500&height=500&srs=EPSG:32735&format=geopackage&format_options=min_zoom:8;max_zoom:21;tileset_name:neighborhood&transparent=true


texts: Keeping name of metrics as band names when using grid_metrics (lidR)

texts:
When executing grid_metrics for a normalized LAScatalog the final results are GTiff files where each band corresponds to a metric. The problem is that the band names are "band 01", "band 02", "band 03"... and I can not identify which band is which metric.
Is there a way to keep the name of the metric as the band name?
Here is the part of the script that I use:

ctg <- readLAScatalog("C:/JAIME/dataLaS_test/", progress = TRUE)

# Load an extra buffer around the processed file to properly process the catalog without edge artifacts
opt_chunk_buffer(ctg) <- 10

# Write files to disc
opt_output_files(ctg) <- "C:/JAIME/00_tests/test2/{ORIGINALFILENAME}_normalized"

#Normalize LAS files
las <- lasnormalize(ctg, tin())

opt_filter(las) <- "-drop_z_below 0"

#export metrics
opt_output_files(las) <- "C:/JAIME/00_tests/test2/{ORIGINALFILENAME}_metrics"

#Create functionn to get metrics
get_metrics <- grid_metrics(las, .stdmetrics_z, 15.8114, filter = ~ReturnNumber == 1)


texts: Setting map rotation in ArcGIS API for Javascript?

texts:
I've created a web map using the ArcGIS Javascript API, but I need my map to have a fixed rotation of 1.76.

The street grid of the map doesn't totally align with the cardinal directions, so it looks slightly off without the rotation.

I've searched through the API reference, but haven't had any luck so far.

Does anyone know how to do this?


texts:
I am trying to add ECW support into Geoserver 2.1.1 (also tried 1.7.7) but no luck so far.

Does anyone have a step by step guide to add the ECW support ? (the one

However, my understanding is that Geoserver 2.1.1 uses imageio-ext that uses gdal-1.4.5

I am now lost between gdal version, geoserver version, imageio-ext version etc etc

This is what my gdal_native dir looks like http://i.stack.imgur.com/yEFpp.png (where libecwj2.dll is copied from http://trac.osgeo.org/gdal/attachment/ticket/2200/libecwj2.zip)


texts: How do you delete a QGIS CRS from the USER database with Python?

texts:
The "Delete CRS in QGIS" question Delete CRS in QGIS is about modifying the recently used list, not actually deleting CRSs.
For Custom USER:xxxxx CRSs you can delete them manually with the '-' icon on the Settings/Custom Projections... page

Is there a python scriptable way to delete a user CRS as the '-' icon does?
I do not see a delete or remove function under https://qgis.org/pyqgis/3.10/core/QgsCoordinateReferenceSystem.html
In other words, how do you do the (hypothetical) user_crs.deleteAsUserCrs() in the sample code below?
crs = QgsCoordinateReferenceSystem("EPSG:27700")
if crs.isValid():
    print("CRS Description: {}".format(crs.description()))
    print("CRS PROJ text: {}".format(crs.toProj()))
else:
    print("Invalid CRS!")
crs.saveAsUserCrs("junk",  QgsCoordinateReferenceSystem.FormatWkt)
if crs.isValid():
    print("CRS Description: {}".format(crs.description()))
    print("CRS PROJ text: {}".format(crs.toProj()))
else:
    print("User CRS is an Invalid CRS!")

crs.deleteAsUserCrs() # <<<<<< what would do this?
if crs.isValid():
    print("CRS Description: {}".format(crs.description()))
    print("CRS PROJ text: {}".format(crs.toProj()))
else:
    print("User CRS is an Invalid CRS!")


texts: 'Number' input parameter can not be used by algorithms in Graphic Modeler

texts:
Does anybody knows why I use the Graphical Modeler in QGIS 3.0, the 'Number'  input parameter can not be used by  algorithms, I have tried the algorithms of 'Buffer', 'create grid' and others, but they can not use the 'Number' parameter I create as the below picture:
(I am not good at English, I hope I describe my question clearly)


texts: Migrating symbols from ArcGIS Server to CartoDB

texts:
We decided to migrate to CartoDB after many years of using ArcGIS Server.
We already handled data migration and now we are looking for an easy to export our symbol rules defined in a MXD to CartoCSS.
Is there any tool available to ease this task?


texts:
I am stuck working with QGIS on an older iMac running OS 10.8.5. I've installed QGIS 2.18.14 and it generally works very well, but I'm running into some issues that get in the way of my work flow. I wonder if I should install an older version of QGIS to get the functionality I need on the older OS... specifically, I am finding that:

Styling layers works strangely - as per this thread from a couple of years ago: QGIS just installed on MAC OSX - vector layer styles do not work properly
Exporting the Atlas as images does not give me the option to choose a file format, as per this 2-year-old discussion - https://issues.qgis.org/issues/12435. It lets me choose the directory to save to, but always uses JPEG. I would prefer to export to PNG to get uncompressed output. My back-up is to export PDFs and convert PNG using other software, but I would rather just use a version of QGIS that worked...

If I should install an older version of QGIS, what version should I choose? 2.8 seemed to help some people with the first issue - but that's quite a few versions back now. Maybe something more recent could work? Trial and error?


texts: Creating polygons from points with undefined order using QGIS

texts:
I have a point shapefile that contain a batch of features. I want to create polygons from these points. Each point should define a boundary of polygon.
There is an unique "ID" field that group points, but I don't have a order field to connect the points by order.
If I run "Points to line" tool I get some zig-zag lines and I don't want it connected in that way, I want each point to define outer boundary of each polygon. This should be done in QGIS.
I have tried Points2One Plugin but I don't have an order field to run this tool and get what I need.


texts: Using Digital Earth Africa collect_training_data tool on image collection data

texts:
On the Digital Earth Africa sandbox, I tried to apply this code to my study. In the code, the tool collect_training_data extracts the pixel values of an input image file, per polygon. I want to use this spectral response information to make spectral signatures for the different classes in my study area. However, the example code uses ga_s2_gm data, which is a single image, as far as I can see, instead of an image collection (all images for the selected year).
I would like to use this tool to get all pixel values per polygon (what it does now), but for all images that I supply the tool; in essence, use the product s2_l2a instead of ga_s2_gm.
If I simply change ga_s2_gm into s2_l2a I get the following error:
      0%|          | 0/2 [01:21<?, ?it/s]
---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<timed exec> in <module>

~/Scripts/deafrica_classificationtools.py in collect_training_data(gdf, products, dc_query, ncpus, return_coords, custom_func, field, calc_indices, reduce_func, drop, zonal_stats, clean, fail_threshold, fail_ratio, max_retries)
    824     # column names are appended during each iteration
    825     # but they are identical, grab only the first instance
--> 826     column_names = column_names[0]
    827
    828     # Stack the extracted training data for each feature into a single array

<string> in __getitem__(self, *args, **kwds)

/usr/lib/python3.6/multiprocessing/managers.py in _callmethod(self, methodname, args, kwds)
    770             dispatch(conn, None, 'decref', (token.id,))
    771             return proxy
--> 772         raise convert_to_error(kind, result)
    773
    774     def _getvalue(self):

/usr/lib/python3.6/multiprocessing/managers.py in serve_client()
    246
    247                 try:
--> 248                     res = function(*args, **kwds)
    249                 except Exception as e:
    250                     msg = ('#ERROR', e)

IndexError: list index out of range

How do I adjust the script to work for s2_l2a instead of ga_s2_gm?


texts: Hover on legend to highlight multiple areas on map

texts:
I would like to highlight different areas on a map when the mouse hovers over the legend using leaflet.

This jsfiddle http://jsfiddle.net/expedio/3w4j2u9z/ is probably the closest to what I am looking for, however, with my limited knowledge of JS and Leaflet I'm struggling to figure out how to link the hover to multiple ID's.

I've searched the internet far and wide and can't seem to find a solution.

Here is an example of the final outcome we would like to achieve - http://atlas.id.com.au/adelaide#MapNo=10173&SexKey=3&datatype=1&themtype=3&topicAlias=population-density&year=2011


texts: Creating a GeoServer user with admin rights to a single store / workspace

texts: Is it possible to create a GeoServer user with admin rights to manage a particular store or workspace? If so how? I've read the online documentation and played with the admin UI but can't work it out.


texts: Mapnik to reuse colors from geotiff when rendering png format

texts:
I am trying to make mapnik read this geotiff file and render it as a PNG file.

Unfortunately all the colors from the geotiff file are not retained. The resulting PNG is completely blank.

This is my Python code:

import mapnik

m = mapnik.Map(800, 600)
m.srs = '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'

s = mapnik.Style()
r = mapnik.Rule()
s.rules.append(r)
m.append_style('Default', s)

ds = mapnik.Gdal(file='K100_2017_1_20_612_62.tif', band=1)

layer = mapnik.Layer('layer')
layer.datasource = ds
layer.styles.append('Default')
m.layers.append(layer)

m.zoom_all()
mapnik.render_to_file(m, 'mapnik.png', 'png32:z=1')

I can get some visible pixels and colors into the PNG by adding a RasterColorizer like this:

rs.colorizer = mapnik.RasterColorizer(mapnik.COLORIZER_DISCRETE, mapnik.Color(0, 0, 0, 0))
rs.colorizer.add_stop(-417, mapnik.Color(0, 0, 0))
rs.colorizer.add_stop(68, mapnik.Color(255, 255, 255))
rs.colorizer.add_stop(234, mapnik.Color(255, 0, 0))
rs.colorizer.add_stop(461, mapnik.Color(0, 0, 255))
rs.colorizer.add_stop(719, mapnik.Color(0, 255, 0))

s = mapnik.Style()
r = mapnik.Rule()
r.symbols.append(rs)

Which renders this PNG:

But then I need to manually define all the colours in the PNG, which just doesn't make sense when all the colours are already defined in the TIFF file. And in the example the PNG-colors are plain wrong and useless.

What am I missing to "reuse" the colors from the TIF in the rendered PNG ?

SOLUTION

Based on input from user30184 this Python code renders correctly:

import mapnik

map = mapnik.Map(800, 600)

style=mapnik.Style()
rule=mapnik.Rule()
rule.symbols.append(mapnik.RasterSymbolizer())
style.rules.append(rule)
map.append_style('Raster Style',style)
lyr = mapnik.Layer('GDAL Layer from TIFF file')
lyr.datasource = mapnik.Gdal(file='K100_2017_1_20_612_62.tif')
lyr.styles.append('Raster Style')

map.srs = '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'
map.layers.append(lyr)
map.zoom_all()

mapnik.render_to_file(map, 'mapnik.png', 'png32:z=1')

It renders this PNG:

Time to study the details...


texts: How to efficiently extract time series from a very large multi-band GeoTIFF in R?

texts:
First some context: I want to run a time-series analysis function that is only available in R for a continental sized multi-band time-series raster, which is stored in disk in GeoTIFF format, that was created translating a VRT file to a GeoTIFF with GDAL using parameter INTERLEAVE='PIXEL' (That parameter should store all band values for each pixel continuously on the disk an make it easy and fast to read.)

The function only needs the time-series to work, so no spatial dependency. And the file has 380 bands, each of continental size, summing up 420GB on disk space.

The problem: I used the raster package to read the file as a stack, but when I try to get the values for a pixel, R just stucks and appears to be processing, but this operation should be fast. I tried using the [] notation, using the extract function and none worked. (Even if I try to get only 3 bands from a small RGB file, R seems slow.)

Using gdallocationinfo it is extremely fast, and with Python and GDAL module it's also incredibly fast.

So shouldn't it be fast in R too? Am I doing something wrong or missing a function designed for that? Doesn't the Raster package use rgdal internally?

But most importantly, is there a way to do it using rgdal then?


texts: Spatial analysis of points with different levels of precision (decimal points)

texts:
I've got a dataset of points with various level of precision - some have two decimal points, other three, other four. Like this:
lat       lon
23.11     120.161
45.6667  -131.12
13.441    108.567

I want to do two types of analysis on these points: (1) calculate distances between points, and (2) calculate centroids of groups of points. But I am not sure if I can do this since they have different levels of precision.
Do I need to round these points to the nearest common decimal (in this case two decimal points) in order to do these analyses on them? Or does it make sense to calculate distances and centroids of points as they are?


texts: PyQGIS function (version 3.4) to delete all fields in a vector Layer except the required one

texts:
I wrote a function that is supposed to delete all fields of a layer except a required one.

My idea was to iterate over all fields of the layer and save all fieldnames except the required one in a List.

NextI want to execute the tool qgis:deletecolumn and pass the list as fields to drop. I'm using QGIS 3.4 that's why I chose 'qgis:deletecolumn' and not "qgis:dropfields'.

This is the code:

def Channel2(self):
    channel1 = self.Channel1()
    fieldname = self.dlg.Channel_ID.currentText()
    fieldNameList = []
    for field in channel1.fields():
        if not (field.name == fieldname):
            fieldNameList.append(field.name)
    param = {'INPUT': channel1,
             'COLUMN': fieldNameList,
             'OUTPUT': 'memory:'}
    channel2 = processing.run('qgis:deletecolumn', param)
    vlayer = channel2['OUTPUT']
    QgsProject.instance().addMapLayer(vlayer)
    return vlayer

I don't get any error message and I do get a new layer as output but the new layer still has all the fields. Checking the documentation it seems that I don't pass a list to the Parameter 'COLUMN' but a "tablefield". But I don't know find what a "tablefield" is supposed to be and how to create it.
Or is the error somewhere else?


texts: "Item Description" from Catalog window in ArcMap won't open

texts:
I can't access metadata information from the catalog window in ArcMap.

I can click on it but nothing happens. I haven't used this option before, so I'm not sure how long I've had this problem. Already checked on ESRI's site to see if there are any patches out for this issue and on GeoNet to see if anyone else is having this issue.

Any ideas?


texts: Cesium | Any Alternatives to 3dtilesets to visualize point cloud data?

texts:
I'm in the early stages of developing a web-based visualization for multiple GIS data sources (imagery, vector data sets, terrain, point cloud...etc).  I am using the Cesium javascript library for this visualization.  And so far, so good...

When it comes to visualizing point cloud data, Cesium has off-the-shelf rendering capability for 3dtilesets generated for a point cloud.  I am hoping to possible sidestep this requirement.  I'm hoping I could just stream some point cloud data to the browser, and render.  I started messing around with PointPrimitive & PointPrimitiveCollection - but initial results don't seem very promising.

Can anyone point me to some examples of rendering point cloud data in Cesium not using 3dtilesets?  I can say upfront I don't intend to visualize a billion points...  I don't have any hard & cold numbers - but something like 1-2 million points would be a good start.


texts:
I am trying to make a RESTful WPS using ZOO Project. According to the following OGC Papers WPS 2.0 Interface Standard and Web Processing Service for WPS 2.0 and 1.0 respectively, it does seem possible to make a RESTful WPS. There is also WPS 2.0 REST API TAMIS which seems to be a RESTful implementation of WPS.

Is there a way to make a RESTful WPS using ZOO Project? (I did not find anything about this mentioned in their documentation)


texts: Opening different XLS spreadsheets as table with ArcPy?

texts:
I am trying to write a Python script that uses a xls file with different spreadsheets as input to build features and  domains  in a geodatabase.

This is the data in Sheet1:

and this the data in Sheet2

and I would like to use this as input for making  features in a geodatabase. Below my script is displayed, but it returns the following error:

  string indices must be integers

#---Import modules---#
    import pandas as pd
    import arcpy
    import xlrd

#---Set paths---#
    SheetFeatureclass = pd.read_excel(r"C:\xlmtofeature\sheets.xls", sheet_name='Sheet1')
    SheetDomain = pd.read_excel(r"C:\xlmtofeature\sheets.xls", sheet_name='Sheet2')
    projection = r'C:\xlmtofeature\RD New.prj'
    out_path = r'C:\xlmtofeature\test.gdb'

#---Colums in Sheet1---#
    FEATURENAME = SheetFeatureclass['FEATURENAME']
    GEO_TYPE = SheetFeatureclass['GEO_TYPE']

#---Colums in Sheet2---#
    DOMEIN = SheetDomain ['DOMEIN']
    FIELDTYPE = SheetDomain['FIELDTYPE']
    DOMEINTYPE = SheetDomain ['DOMEINTYPE']

#---Create features ---#
    for row in SheetFeatureclass:
        try:
            arcpy.CreateFeatureclass_management(out_path, row[FEATURENAME], row[GEO_TYPE], "", "", "",projectie,"","","", "")
            print str(row[FEATURENAME]+ " Toegevoegd" )
        except Exception as err:
            print(err.args[0])
    print "Done with adding features"

#--- create Doamins ---#
    for row in SheetDomain:
        try:
            arcpy.CreateDomain_management(out_path, row[DOMEIN],"",row[FIELDTYPE],row[DOMEINTYPE],"","",)
            print str("Add "+ row[DOMEIN])
        except Exception as err:
            print(err.args[0])
    print "Done with adding Domains"


texts: QGIS 3.10 - Point cluster symbology sum attribute

texts:
I have a point layer.

Each feature represents a case where a soil receiving facility has received a weighed amount of soil, this amount is for each feature written in the attribute "Faktisk mængde".
Geographically the points overlap since there is a limited number of soil receiving facilities.

Therefore I use the point cluster symbology.

Automatically the clusters are styled to represent the number of features in the cluster, but i would also like to label the clusters with the sum of the attribute "Faktisk mængde".

Is this possible?

I have tried using the clustering tools to make a new layer as proposed, in this question, but using that tool or DBSCAN does not seem to work even though the features are perfectly overlapping.


texts: Testing if shapefiles intersect from a list using intersect: ArcGIS

texts:
I am new to programming and scripting in general, so please bear with me.  I have generated a list of the location of all the shapefiles for my workplace.  The files are all placed in folders, no geodatabases are used.  I am trying to create a tool which queries which of the shapefiles intersect another shapefile.  I would like to use the list as one of the inputs for the intersect parameters.  I want to:

Calculate the area of the intersected features
perform stats
If area of intersect > 0, convert the original data to feature layer
clip the feature layer and save as new fc in target location.

I would really appreciate some python script examples as to go about this task.

I am using arcGIS 10 with an arcView license.

Thank you.


texts: Rice Field Levee Contours - .2 foot interval needed - LiDAR data?

texts:
I am a rice farmer in Arkansas and would like to create levees along contour lines in a field at a .2 foot interval. How could I utilize LiDAR data to map 0.2 foot contour intervals across my area of interest?

I got 2014 Lidar data from USGS Earth Explorer. 900 mb.


texts: Performing a 'far analysis' using ArcGIS for Desktop?

texts:
I would like to perform an analysis that is the inverse of Near_analysis for ArcGIS. Instead of my result being NEAR_DIST and NEAR_FID, I'd like each feature to have a FAR_DIST and FAR_FID populated with the furthest feature information from a 'far feature'. I can do this with some scripting, but it would be slow, and involve iterating through each feature in a feature class and performing a Near_analysis on each, and picking the furthest distance.

Any other ideas? I'm ArcGIS 10.1 and proficient in python/arcpy.


texts: Debugging ArcPy ImportToolBox AttributeError: Object: Tool or environment  not found?

texts:
I am trying to use the ConcatenateRowValues tool provided by ESRI. My overall code for my script is working fine except when it tries to use the ConcantenateRowValues Tool. I start by importing the proper toolbox and then attempt to get my script to access it. Below is my code to do this:

arcpy.ImportToolbox("C:/Users/'USERNAME'/Desktop/BaseLayer/ConcatenateRowValue101/ConcatenateRowValues101.tbx")
arcpy.ConcatenateRowValues_AA(hlint, "SHORT_DESC", "NAME", "Hydro_Line", "-")

The Error that's returned is:

Runtime error
Traceback (most recent call last):
  File "<string>", line 67, in <module>
  File "C:\Users\'USERNAME'\Desktop\BaseLayer\ConcatenateRowValue101\ConcatenateRowValues101.tbx", line 54, in ConcatenateRowValues
    /jVXV+BHwwJ5+Gfgv/wnbP8A+NV2Eo+ZTjHPQ9qlR60cIvojBSkupxw+Avwvz/yTTwX/AOE7Z/8A
AttributeError: Object: Tool or environment <ConcatenateRowValues_AA> not found

This is a different error from when an invalid function name is called on arcpy:

>>> arcpy.ConcatenateRowValuesASDF_AA(hlint, "SHORT_DESC", "NAME", "Hydro_Line", "-")
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
AttributeError: 'module' object has no attribute 'ConcatenateRowValuesASDF_AA'

That means ArcPy is seeing the model but not wanting to run it. Why not? What's the reason for the Tool or environment <ConcatenateRowValues_AA> not found error?


texts: Exporting inset map with same quality in ArcGIS Desktop?

texts:
I am having trouble exporting my map that has the same quality in the inset.  The lines come out blurry.

You can see a big difference between the quality.  Any ideas?


texts: How to cut polygons that overlay other polygons in ArcMap?

texts:
I created a layer of buffer polygons following some conditions (distance from parks, stations etc). Some of the buffer polygons are now overlaying public_parks polygon layer.

What tool can I use in order to create a new layer where the buffer polygons are not covering the parks polygons?

I thought about using NOT boolean, but from what I know it works only on raster data.

Thank you


texts: How to Receive an Email after a script has run?

texts:
I'm working my way through my first attempt at scheduled tasks with modelbuilder and python. So far I have:

Created a model that reconciles our SDE versions and then compresses the database.
Exported the model to a script so I can set it up in Windows Task Scheduler

Its a fairly straight forward model so I expect it to work, however, I'd like to set it up to email me when it is completed and include the geoprocessing messages.

I found this support article that includes a script that does just that.

My question is: How do I include this or reference this in my reconcile and compress model?


texts: Correct coordinate systems, but trustworthy distance between points?

texts:
I have limited GIS experience with projections and I have some questions about trying to accurately display point data from two different sources in ArcGIS 10.1 in order to do a proximity analysis.

The first data set, which already contains Latitude and Longitude coordinates, is from a state government agency and its metadata lists Horizontal Coordinate Scheme (projection): Not Applicable, Ellipsoid: GRS80, Horizontal Datum: NAD83, Horizontal Units: Feet, Distance Resolution: Not Applicable, and Geographic Coordinate Scheme: Decimal Degrees.

The second data set contains addresses that must be geocoded.  I used the 10.0 US Streets Geocode Service (ArcGIS Online) from the Geocoding toolbar to create a point layer based on this second dataset.

Everything appears to line up perfectly well, but since the goal is an accurate proximity analysis, I hoped to confirm that the distances among points seemed plausible before doing more work.

I have tried different approaches for projections, ranging from assuming the geocoded points were WGS84 and then projecting as NAD83 to trying to re-project both of the layers to a Projected Coordinate System such as a relevant State Plane NAD83 system.  Again, the points look fine and the "Source" coordinate system information matches across both layers and the data frame.  The problem is that, no matter what I do, the distance tool is giving me what seem to be extremely low distance values compared to an online tool that measures linear distance between points.  For example, the online tool reports a linear distance of 6,214 feet between 2 reference points, while different projection approaches with ArcGIS are giving me distances of 4,418.753176 feet and 4,384.41136 feet between the same 2 reference points.

Here are my questions:

1) Does the 10.0 US Streets Geocode Service automatically display points with WGS84?  If so, is the correct strategy to then define the Geographic Coordinate System for the points layer as WGS84 and, after that, use the Projections and Transformations->Feature->Project tool to display in NAD83?

2) Is it a mistake to think of a website that uses Google Maps to calculate linear distance between points as some sort of verification?  Are the reported distances on the website also flawed due to the "Web Mercator" projection used for Google Maps?  If it is a mistake, can anyone suggest a better resource?  Does an 'authoritative' resource even exist?

Any feedback or advice would be greatly appreciated.  Apologies if I'm incorrectly using some of the terminology or missing something very basic.  Thank you!


texts: Checking if points fall within polygon Shapefile

texts:
Zillow has a set of shapefiles for different neighborhoods of major US cities. I wanted to check if certain buildings were present in certain neighborhoods using R:

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

df <- data.frame(Latitude =c(47.591351, 47.62212,47.595152),
                 Longitude = c(-122.332271,-122.353985,-122.331639),
                 names = c("Safeco Field", "Key Arena", "Century Link"))
coordinates(df) <- ~ Latitude + Longitude

wa.map <- readOGR("ZillowNeighborhoods-WA.shp", layer="ZillowNeighborhoods-WA")

sodo <- wa.map[wa.map$CITY == "Seattle"  & wa.map$NAME == "Industrial District", ]

I can plot without any issues

plot(sodo)
points(df$Latitude ~ df$Longitude, col = "red", cex = 1)

I match the proj4 string from the shapefile to my data.frame

CRSobj <- CRS("+proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0 ")
df@proj4string <- CRSobj

over(df, sodo)

This just gives me a bunch of NA values. I have tried this answer

spp <- SpatialPoints(df)
spp@proj4string <- CRSobj
over(spp, sodo)

but still get only NA values. Any ideas what else I should try?


texts:
I have a RasterStack that I want to add two rasters to.  I am using this code

library(rgdal)
library(raster)

#original rasterstack
rasterfiles   <- list.files("C:\\path_to_rasters", "*.tif$", full.names = TRUE)

#stack the bands
stacked <- raster::stack(rasterfiles)

#read in new rasters
aspect <- raster("C:\\Aspect.tif")
slope <- raster("C:\\Slope.tif")

#create a shapefile extent of the stacked raster
e <- extent(stacked)

# coerce to a SpatialPolygons object
p <- as(e, 'SpatialPolygons')

#assign crs
projection(p) <- crs(stacked)

#function to mask the raster
mask_raster <- function(raster, shape){
      outr <- crop(raster, extent(shape))
      outr <- mask(outr, p)
      return(outr)
}

#apply function
slope <- mask_raster(raster= slope, shape = p)

aspect <- mask_raster(raster= aspect, shape = p)

#add new layers to stack
stacked <- addLayer(stacked, aspect)
stacked <- addLayer(stacked, slope)

This works completely fine for aspect but for slope when I try to add the layer I get the following error:

Error in compareRaster(c(x, rasters)) : different extent

The crs of all rasters are identical.

Before trying to mask this is what is in each object:

print(stacked) returns:

class       : RasterStack
dimensions  : 5629, 253, 1424137, 428  (nrow, ncol, ncell, nlayers)
resolution  : 1, 1  (x, y)
extent      : 480524, 480777, 5220075, 5225704  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=14 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
names       :   Band_1,  Band_10, Band_100, Band_101, Band_102, Band_103, Band_104, Band_105, Band_106, Band_107, Band_108, Band_109,  Band_11, Band_110, Band_111, ...
min values  :     0.00,    91.00,     0.00,     0.00,     0.00,     0.00,     0.00,     0.00,     0.00,     0.00,     0.00,     0.00,    80.00,     0.00,     0.00, ...
max values  :  1914.00,  6061.00,  8805.00,  8830.00,  8745.00,  8455.00,  8230.00,  8317.00,  8092.00,  7945.00,  8124.00,  8122.00,  6160.00,  7642.00,  7222.00, ...

print(aspect) returns:

class       : RasterLayer
dimensions  : 21152, 22837, 483048224  (nrow, ncol, ncell)
resolution  : 1, 1  (x, y)
extent      : 475140, 497977, 5212857, 5234009  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=14 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
data source : C:\NEON_2\lidar\Aspect.tif
names       : Aspect
values      : 0, 360  (min, max)

print(slope) returns:

class       : RasterLayer
dimensions  : 23239, 29315, 681251285  (nrow, ncol, ncell)
resolution  : 1, 1  (x, y)
extent      : 471887, 501202, 5211828, 5235067  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=14 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
data source : C:\NEON_2\lidar\Slope.tif
names       : Slope
values      : 0, 65.33745  (min, max)

After I mask each object contains this:

print(aspect) returns:

class       : RasterLayer
dimensions  : 5629, 253, 1424137  (nrow, ncol, ncell)
resolution  : 1, 1  (x, y)
extent      : 480524, 480777, 5220075, 5225704  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=14 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
data source : in memory
names       : Aspect
values      : 0, 360  (min, max)

print(slope) returns:

class       : RasterLayer
dimensions  : 5629, 253, 1424137  (nrow, ncol, ncell)
resolution  : 1, 1  (x, y)
extent      : 480524, 480777, 5220075, 5225704  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=14 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
data source : in memory
names       : Slope
values      : 0.0285432, 39.09481  (min, max)


texts: Adding the details of the map maker into qgis map

texts:
I am making maps as part of my academic work.
It is important to have the details of the map maker on these kind of maps, and I am guessing it is also true for maps that are made for clients.

I was wondering is there an easy way to add my details, so they will we be in the corner of the map and up to scale?

Something like the way to credit OSM maps in QGIS.


texts: Converting Heatmap into Grid of Points to regenerate less accurate Heatmap?

texts:
Question: How do you take a heatmap generated using a GIS software (arcgis, qgis...) and convert it into a uniform grid array of points which can be used to regenerate the heatmap (although with lesser accuracy)?

Initial to final product:

(Initial)Millions of Datapoints in MySQL -> (... some processing, maybe GIS software) -> (Final)Hundreds of Datapoints in a grid array -> Render density map in browser

How do I get from Millions of Datapoints in MySQL to Hundreds of Datapoints in a grid array?

Purpose: Transferring millions of points from database to browser where the heatmap is created will cause heavy database load and network traffic and strain the browser. So instead of transferring millions of points to the browser, it is more optimal to simplify the millions of points to 100s of points in a grid arrangement which will be easier for the browser to handle.


texts: QGIS dissolves all different shapes into one feature instead of multiple features

texts:
I'm pretty new to QGIS and geoprocessing. I have a layer of points representing traffic lights - where each point is a different feature:

I wish to create a layer with traffic light junctions based on this layer. My approach is to create buffer around each point - which I've done using QGIS:

But then, when I try yo dissolve them into single shapes, QGIS creates one feature for all the shapes instead of a feature per junction:

How can I tell QGIS to create a feature per junction?


texts: OpenLayers: my proxy truncating response?

texts:
I'm using OpenLayers on IIS 7.5 and I want to get my proxy (OpenLayers.ProxyHost) working on my localhost machine. So i copy/pasted the source for this example and have it running on localhost. However my WFS layer is not showing up.

I took a proxy code from here i am using:
http://code.google.com/p/iisproxy/issues/detail?id=8

Now everything seems to work fine except my response seems to be truncated. Using firebug i can see my request goes through perfect and my RESPONSE is truncated for some reason.

What's going on and HOW CAN I FIX THIS?

Here is my proxy code source: (I've commented out credential stuff as i'm using forms authentication).

<%@ WebHandler Language="C#" Class="proxy" %>

using System;
using System.Configuration;
using System.IO;
using System.Net;
using System.Web;

/// <summary>
/// adapted from: http://code.google.com/p/iisproxy
/// </summary>
public class proxy : IHttpHandler
{
    public void ProcessRequest(HttpContext context)
    {
        var request = context.Request;
        var remoteUrl = request.QueryString["url"];

        var req = (HttpWebRequest)WebRequest.Create(remoteUrl);
        req.AllowAutoRedirect = false;
        req.Method = request.HttpMethod;
        req.ContentType = request.ContentType;
        req.UserAgent = request.UserAgent;
        //var basicPwd = ConfigurationSettings.AppSettings.Get("basicPwd");
        //req.Credentials = basicPwd == null ?
        //    CredentialCache.DefaultCredentials :
        //    new NetworkCredential(HttpContext.Current.User.Identity.Name, basicPwd);
        req.PreAuthenticate = true;
        req.Headers["Remote-User"] = HttpContext.Current.User.Identity.Name;
        foreach (string each in request.Headers)
            if (!WebHeaderCollection.IsRestricted(each) && each != "Remote-User")
                req.Headers.Add(each, request.Headers.Get(each));
        if (request.HttpMethod == "POST")
        {
            var outputStream = req.GetRequestStream();
            CopyStream(request.InputStream, outputStream);
            outputStream.Close();
        }

        HttpWebResponse response;
        try
        {
            response = (HttpWebResponse)req.GetResponse();
        }
        catch (WebException we)
        {
            response = (HttpWebResponse)we.Response;
            if (response == null)
            {
                context.Response.StatusCode = 13;
                context.Response.Write("Could not contact back-end site");
                context.Response.End();
                return;
            }
        }

        context.Response.StatusCode = (int)response.StatusCode;
        context.Response.StatusDescription = response.StatusDescription;
        context.Response.ContentType = response.ContentType;
        //if (response.Headers.Get("Location") != null)
        //{
        //    var urlSuffix = response.Headers.Get("Location");
        //    if (urlSuffix.ToLower().StartsWith(ConfigurationSettings.AppSettings["ProxyUrl"].ToLower()))
        //        urlSuffix = urlSuffix.Substring(ConfigurationSettings.AppSettings["ProxyUrl"].Length);
        //    context.Response.AddHeader("Location", request.Url.GetLeftPart(UriPartial.Authority) + urlSuffix);
        //}
        foreach (string each in response.Headers)
            if (each != "Location" && !WebHeaderCollection.IsRestricted(each))
                context.Response.AddHeader(each, response.Headers.Get(each));

        CopyStream(response.GetResponseStream(), context.Response.OutputStream);
        response.Close();
        context.Response.End();
    }

    static public void CopyStream(Stream input, Stream output)
    {
        var buffer = new byte[1024];
        int bytes;
        while ((bytes = input.Read(buffer, 0, 1024)) > 0)
            output.Write(buffer, 0, bytes);
    }

    public bool IsReusable
    {
        get { return false; }
    }
}


texts: What constitutes a valid spatial layer when connecting QGIS to MySQL?

texts:
I'm trying to work with a MySQL 5.5 table using QGIS 1.8.  I can do "Add Vector Layer", choose "Database", and successfully connect to the MySQL connection I previously set up.  QGIS shows me a "Select vector layers to add" dialog that lists the tables in my database, including the one I want to use which shows a "Geometry type" of "Point" (which is correct, there's a POINT column in that table).  When I highlight that table from the list and click OK I get an error message:

The layer MySQL:[redacted connection info] is not a valid layer and can not be added to the map

Is there something that needs to be done to make a MySQL table a "valid layer" other than making sure that it has a geometry field (like a point)?


texts:
I have a raster file containing roads (left image). I polygonize the raster file with GDAL (see the script below). In the end I like to have vector lines. However GDAL only gives me something like this back (right image). Obviously this is correct, since polygonize creates polygons. Is there a way to "linize" a raster?

Here is my code:

import gdal,ogr,os

# open raster file
raster = gdal.Open('test.tif')
band = raster.GetRasterBand(1)

#create new shp file
newSHPfn = 'test.shp'
shpDriver = ogr.GetDriverByName("ESRI Shapefile")
if os.path.exists(newSHPfn):
    shpDriver.DeleteDataSource(newSHPfn)
outDataSource = shpDriver.CreateDataSource(newSHPfn)
outLayer = outDataSource.CreateLayer(newSHPfn, geom_type=ogr.wkbLineString )

# polygonize
gdal.Polygonize(band, None, outLayer, 1)


texts: Maximum Likelihhod Classification - Error 010108 and Error 010067

texts: I'm trying to run a supervised classification on a Landsat8 multi-band raster in ArcMap. Every time I run the classification it returns Error 010108 and Error 010067. The 010108 error says "The number of bands in the signature file is less than the number of bands in the raster", so I tried using the Composite Bands tool to compress my raster to a single band, but it still returns a raster with the same number of bands as the original. Is compressing my multi-band raster to a one band raster the right solution? If so, how do I do so without extracting bands individually? I can't find any other suggested tool than Composite Bands and it's still just producing an identical multi-band raster.


texts: Variation of algorithms to find point from another point with bearing and distance

texts:
There are several algorithms to calculate distance between two latitude/longitude points.

Haversine
Hubeny
Lambert-Andoyer
Vincenty's formulae for inverse problem

But there are not so many algorithms (actually I can't finally find it) which can calculate latitude/longitude of point B from another latitude/longitude point A with bearing and distance.

Vincenty's formulae for direct problem

I'm glad to know if there are any algorithms which is lesser accurate, but faster than Vincenty's formulae.


texts:
How can I auto clip text and append dots on a label if the current text is too long in GeoServer?
Like in GoogleMaps:

I only used AutoWrap option.


texts:
I have the following table defined, which uses PostGIS geography.

CREATE TABLE test_geog (
    id SERIAL PRIMARY KEY,
    boundary GEOGRAPHY(Polygon)
);

I added the following test polygon to the table:

INSERT INTO test_geog VALUES (ST_GeographyFromText('SRID=4326;POLYGON((0 0, 0 5, 5 5, 5 0, 0 0))'));

I am trying to determine whether a point lies within any of the polygons in this table.  I have this query:

SELECT ST_ContainsProperly(ST_GeographyFromText('Point(2 2)'), area)
FROM (SELECT boundary FROM test_geog) AS area;

This yields the following error:

ERROR:  function st_containsproperly(geography, record) does not exist
LINE 1: SELECT ST_ContainsProperly(ST_GeographyFromText('Point(2 2)'...
            ^
HINT:  No function matches the given name and argument types. You might need to add explicit type casts.

How do I convert this "record" into a POLYGON?  I'm confused because it seems like the column has already been declared to only hold POLYGON types, but for some reason that is not what I am pulling out of the database.

I attempted to cast the record to a POLYGON like this:

SELECT ST_ContainsProperly(ST_GeographyFromText('Point(2 2)'), CAST (nitf_area AS POLYGON))
FROM (SELECT boundary FROM source_imagery) AS nitf_area;

But that gives me this error:

ERROR:  cannot cast type record to polygon
LINE 1: ...tainsProperly(ST_GeographyFromText('Point(2 2)'), CAST (nitf...

What am I not understanding here?


texts: What do you think is happening to my QGIS Georeferencer

texts:
I'm pretty new with QGis, but have used tools as Openlayers and Geoserver for some years now without never being asked to diplay a georeferenced raster image. Which gets me to use QGis in order to build a Geotiff from a simple scanned PNG or JPEG.

I've followed many tutorials, but I'm somehow still struggling to import any image file in my Georeferencer plugin window.

Here's what I get when I import a jpeg image:

It looks like a badly read jpeg. Any idea ?

Thanks for your help !

EDIT

I've tried with a PNG by converting the original jpeg with an image tool (like IrfanView) and the image is now diplayed.
It's better, but doesn't answer exactly why the jpegs aren't working ...


texts: See highlighting when an overlapped point is selected?

texts:
I have two points with the same location (with the same geometry point). It seems that I'm dealing with a overlapped selection bug because having those two points when I choose the second (the one on top) through the attribute table will be shown a highlight selection which is OK but if I choose the first one (the one below) it simply doesn't highlight. Which it is very confusing since that it means that the first point is a point without geometry, which is not true. When you have overlapped points QGIS will select only the top ones (in my case the last record on my table).

Is this a bug or in fact is a property that I can change in some place on QGIS settings? If so, where can I change that to show always the highlight selection for overlapped points?

Note: My QGIS version is 2.18.16


texts: How to use ST_DelaunayTriangles to construct a Voronoi diagram?

texts:
(edit 2019) ST_VoronoiPolygons available since PostGIS v2.3!

With PostGIS 2.1+ we can use ST_DelaunayTriangles() to generate a Delaunay triangulation, that is a dual graph of its Voronoi diagram, and, in theory, they have an exact and reversible conversion.

Does any safe  SQL-standard script with an optimized algorithm exist for this PostGIS2  Delaunay-to-Voronoi conversion?

Other refs: 1, 2


texts: GDAL C# Error: "Failed to import coordinate system"

texts:
I've been struggling with building C# code that will do the equivalent of gdalwarp, specifically reproject a GeoTIFF file and save to an ESRI .bil file (EHdr)

To that end, I thought the following code might do the trick:

var vrt = Gdal.AutoCreateWarpedVRT(srcDs, null, epsg5070, ResampleAlg.GRA_Average, 0.0);
var bilfile = bilDriver.CreateCopy(dstPath, vrt, 1, null, null, null);

where epsg5070 is the WKT for EPSG:5070.  The source GeoTIFF is EPSG:4269.

The AutoCreateWarpedVRT statement errors with:

Failed to import coordinate system `GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS 1980",6378137,298.2572221010042,AUTHORITY["EPSG","7019"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6269"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4269"],EXTENSION["CENTER_LONG",-96.5]]'.

The WKT in the error is 4269, the projection of the source file.

I am using the NuGet packages (GDAL / GDAL.Native, v2.4.4), have called the GdalConfiguration methods, the GDAL files appear to be copied correctly to the project output (including the gdal/data and gdal/share folders), and it appears the paths are being set correctly.  But my impression is that GDAL is not finding its data.

Where should I look next?  Is there a trick I don't know about?  Or am I going about this the wrong way?
 I'd rather not shell out to the command-line utilities.


texts: Where do I find information where I can download the legislative mapping data

texts:
Possible Duplicate:
  Where I can download legislative mapping data?

It is a collection of files downloaded from aristotle to map Zip to District data.
Among the files are
1.  Census Code (ILCENSUS) which has County Census Code, State Census Code and State County Name
2.  Legislator Identification (IL105ZD) which has district address information, legislator number and name.
3.  Legislative district mapping(ILLEGYY) which has County Code, State Census Code and  Zip.


texts:
I'm trying to build internal overviews on a big GeoTIFF using gdaladdo and it is way too slow.

From what I could find out it is a known problem (#1442) that has been solved in May last year

However in my case it is still excruciating slow. I'm wondering what I could do wrong (I use GDAL 3.0.2 released on 28/10/2019).

I tried to generate 1 level of external overview with the following command. It lasted 68 min which is fine

P:\>gdaladdo -ro -r average --config BIGTIFF_OVERVIEW YES --config INTERLEAVE_OVERVIEW PIXEL --config GDAL_TIFF_INTERNAL_MASK YES --config PHOTOMETRIC_OVERVIEW YCBCR --config COMPRESS_OVERVIEW JPEG --config JPEG_QUALITY_OVERVIEW 85 "C:\superorto_ext.tif" 2

With the exact same command on the same file but with internal overviews (without the -ro switch)  it lasted 6 hours and 50 min which is not so fine

P:\>gdaladdo -r average --config BIGTIFF_OVERVIEW YES --config INTERLEAVE_OVERVIEW PIXEL --config GDAL_TIFF_INTERNAL_MASK YES --config PHOTOMETRIC_OVERVIEW YCBCR --config COMPRESS_OVERVIEW JPEG --config JPEG_QUALITY_OVERVIEW 85 "C:\superorto_ext.tif" 2


texts: Creating many-to-many table from intersection between two polygon layers in QGIS?

texts:
I have created a buffer around several points. I want to intersect this layer with a Census Tract layer. My end goal is to create a table that would list all the census tracts that intersect with the buffer area of each point. Something like this:

The result must show that some buffer entities may be related to many census tracts, and the other way round. I think ArcGIS' "Relate" function is suited for that purpose, but I don't know how to do this in QGIS.

Do you have any idea how I could do that in QGIS 2.6?

UPDATE:
The way I am planning on doing this if I don't get a simpler answer is to

Create a new attribute in the Census Tract shapefile for each of the 15 buffer polygons
For each of the 15 buffer polygons: run a spatial request to select all Census Tracts that intersect with the chosen buffer polygon
Use the Field Calculator to write the name of the Census Tract in the variable I created in Step 1
Export the Census Tract shapefile to Excel and play around..

There has to be a simpler way around, but I don't know it.


texts: Get all city information (roads, buildings, ..etc) in OpenStreetMap via Overpass API

texts:
I want to draw a city (Avignon, France) with JavaFX and OpenStreetMap (via Overpass API).
What is the QL request that I have to send to Overpass API to get this information?


texts: Dynamic connection between Mapinfo and a database

texts:
I'm trying to create a dynamic connection between a database and a Mapinfo table. My database countain only one table, with a list of bus stops, characterized with their ID, name, and X Y coordinates. I already have created a Mapinfo table that updates with the database, using a DBMS connection, and I have made it mapable so the stops are placed on a map at the right coordinates.

But I can't make it work the other way; The table created with the DBMS connection is read only. I'm trying to make the Mapinfo Table editable, and to apply the changes of the table on the database.

Can anybody help me out?

(I use Mapinfo Professional 12.5 and PostgreSQL)


texts: "Object based Image Segmentation" and some options are absent in SAGA

texts: Although I have the latest version of SAGA, "Object based Image Segmentation" and some options do not appear when I click on "Geoprocessing>Imagery> Segmentation>? " Why are they absent?


texts: Solar radiation flux at 5 minute interval using GRASS GIS r.sun

texts:
I have been trying to estimate instantaneous total solar irradiance(w/m2) at  five minute interval for a day (eg. 21 sept). I can do this in SAGA GIS by using "moment, mode 1" and running the solar module for 12*24 time.

Therefore, for automation purpose I want to stick to GRASS. I have looked into r.sun and r.sun.hourly however I am unable to use these function to produce output at 5 minute interval.

Here is the sample code from r.sun.hourly

r.sun.hourly -t elevation=example_data@PERMANENT start_time=6 end_time=18
             time_step=.08 day=277 year=2011 beam_rad_basename=beam
             diff_rad_basename=dif glob_rad_basename=tot nprocs=6


texts: QGIS problem with exporting from print composer

texts:
I am using QGIS (version 2.10.1-Pisa) on a Mac (Yosemite). I am trying to export a map in the print composer as jpg or pdf and whenever I do try, the page in the print composer goes blank. when trying to save as jpg, a file is saved but the image is black. it also saves a pdf but the file is empty.
I am basically trying to print some simple vector layers on top of a google hybrid background.


texts: How Do You Autopopulate Attribute Based on Spatial Location of Underlying Grid?? ArcGIS Mobile

texts:
User drops a point, then Attribute Table Autopopulates based on Underlying (Polygon)Grids Spatial Location.

I want to do this for field editors using ArcGIS Mobile.

Does anyone have an idea on how to do this?


texts: Circle radius is not being applied in ArcGIS API for JavaScript?

texts:
I am working with the ArcGIS API for JavaScript 4.3

Everything is going well except this one little issue.

I am trying to allow a user to add a location to our db.
When they put in the address the [x,y] coords are created from the address to make a point

I then take the point, make it the center of a circle and then set a radius of that circle and add it to the map.

The problem:
No matter what I do I only get a circle on the map without the radius.

Here is my code.

NearbyLocationsMap.prototype.GetPointsNearby = function (x, y, cb) {
    var PointsNearbyPromise = new Promise(function (resolve, reject) {
        resolve(1);
    });
    PointsNearbyPromise.then(function () {
        return point = nl.MakePointFromXY(x, y);
    })
    .then(nl.MakeCircleFromPoint)
    .then(nl.MakeCircleLayer)
    .then(nl.GoToCircle);
}

NearbyLocationsMap.prototype.MakePointFromXY = function (x, y) {
    return new nl.Point({
        x: x,
        y: y
    })
};

NearbyLocationsMap.prototype.MakeCircleFromPoint = function (point) {
    var radius = 25;
    return new nl.Circle(point,{
        radius: radius,
        radiusUnit: "miles",
        geodesic: true
    })
};

NearbyLocationsMap.prototype.MakeCircleLayer = function (circle) {
    nl.CircleLayer;
    nl.CircleLayer = new nl.GraphicsLayer({id: "circles"});
    nl.map.add(nl.CircleLayer);
    nl.N_graphic = new nl.Graphic({
        geometry: circle,
        symbol: nl.getRendererCircle()
    });
    nl.CIRCLE = circle;
    nl.CircleLayer.add(nl.N_graphic);
}

NearbyLocationsMap.prototype.GoToCircle = function () {
    var c_layer_extent = nl.CIRCLE.extent;
    nl.c_layer_extent = c_layer_extent;
    return setTimeout(function () {
        nl.view.goTo({
            target: nl.c_layer_extent,
            zoom: 14
        }, {
            duration: 1000,
            easing: "linear"
        }, 500)
          .then(nl.findLocationsInCircle)
          .then(LocationModalLoading);
    });
};

NearbyLocationsMap.prototype.getRendererCircle = function () {
    return new nl.SimpleMarkerSymbol({
        style: "solid",
        color: "rgba(255, 255, 255, 0.5)",
        outline: {
            width: 2,
            color: "rgba(24, 69, 59, 1)",
        }
    });

    //return new nl.SimpleLineSymbol({
    //    color: "rgba(24, 69, 59, 1)",
    //    width: "2px",
    //    style: "short-dot"
    //});
}