texts: For each polygon, return yes/no value if it intersects another layer

texts:
I have a polygon grid overlaying several other polygon vectors. I wish to visit each polygon in the grid and ask yes/no if each of these layers fall within the grid cell or not. I can do this easily with a spatial query to select the appropriate cells, but I want to have this information returned as a new attribute column, where each of the grid cells has a yes/no (or 1/0) value accordingly. This seems simple but I cannot figure out how to do it.

Thank you very much for any guidance


texts: Is there any tutorial about creating a heat-map in ArcMap 10.1?


texts:
I’m very new to GIS, and I’ve been reading a lot of manuals and information, but I thought I’d ask a couple of questions here.

One of my major hobbies is drawing maps of a completely imaginary country; it’s part of a massive project of mine. I’m getting to the point where I’m reaching practical limitations in the way I’m drawing the maps using traditional graphics software.

However, I’ve already done quite a lot of work on these SVG maps in Inkscape, and I really would prefer not to start completely from scratch: the current file represents several months, if not a couple of years worth of work. They aren’t georeferenced, but they are drawn to scale (1:10000).

Is there a method to convert or import the SVG into a GIS application, such as QGIS? I’ve been playing around with it all week to no avail. Ideally I want to just import the file, layers and all, and pick up drawing where I left off, but that seems to be all-but-impossible. Even just an easy way to add georeferencing to the exported image would be a good start.

EDIT: Here’s one of the maps in question. It’s an export from Inkscape at 60 dpi; the original SVG file is ~ 3×3 metres(!) in size and tends to crash browsers. I’ve been documenting my progress with an album too. Although it’s unfinished, I figure if I’m planning to switch to QGIS, the earlier I do so the less re-mapping I’m probably going to have to do too. The bottom left corner is the point 0°E 40°N, and the map covers an area of roughly 1100 km² (33×33 km).


texts: Are there any geology symbol libraries for QGIS?

texts: Are there any structural-geology symbols libraries (and geology in general) for QGIS?


texts: Create a point x meters from origin using st_translate in postgis

texts:
How do I create a point 1000 meters away from a reference point?

The ST_Translate only supports geometry type.

I also tried using ST_GeomFromText but that returns a unit and not in meters.

Here is my sample query:

SELECT ST_AsText ( ST_Translate ( ST_GeogFromText('SRID=4326;POINT(-70.01 43.27)'), 1000, 0 ) );

Is there a way that I can convert the output of ST_Translate into meters?


texts: How to import shapefiles into MS SQL 2008 and then view that data using QGIS?

texts:
I'm a QGIS users but not a database person and I need some help on the DB side.
I'm currently in the need to import several shape files into MSSQL 2008 and I have no clue how to do that and my IT guy either.
Can anyone explain to me the components needed,and the step by step how to proceed, or point me towards an online tutorial for unexperienced DB person like me please.

Furthermore, lets assume that all my shape files are imported into MSSQL 2008 that is installed on an internal secured server, Do I need to install another component to have access to my shape files from outside my organization (ie: from a remote location via the web ) ? or will I be able to access my data easily using the QGIS Browser ?. If yes, what do I need to do (or tell my IT guy to do !!) ?

As you can tell, I'm a real rookie when it comes to DB servers...I just hope soemone can point me towards the right direction so I can help my IT guy in the process.

Thanks for all your help.
Yann


texts: Cleaning blocks pattern from urban LIDAR (blocks uplifted from streets)

texts:
We have a 1-meter LIDAR DEM from a city.

A small subset can be downloaded from this link:

This screenshot shows raw DEM with gray palette (darker belts are streets, and greyish and whitish rectangles are blocks):

This correspond to a place in Santo Domingo city, which can be seen in this Google screenshot:

On average, blocks are "uplifted" ca. 2 meters from streets, which is not true. We want to have a clean DEM to generate stream network and topographic wetness index (TWI). With the DEM supplied (we don't have the original bands from laser scanner), the hydrgraphic network seemed to follow a rectangular layout, and TWI resulted in blocks pattern. These pictures show the results:

This is the stream network result, generated with r.watershed in Grass GIS:

And this is TWI result, generated with SAGA:

We tried some procedures to solve this inaccuracy without success:

1) Denoising tool. We applied r.denoise tool in Grass GIS, but got some problems with the module installation. We ran it again with a shell in Windows, and received an insufficient memory message.

2) Filters. We ran different types of filters (low-pass, median, mean, etc.), with different windows sizes, and trying to put weight in the direction of the streets (Grass GIS, SAGA, QGIS).

3) Geostatistics. We generated points cloud strictly over streets (tried 1000 and 2000 points), generated a variogram model and then ran an ordinary krigging to fill blocks. Variogram modelling, and ordinary krigging, was done in R, using different packages. We got a linear variogram, so we wouldn't rely on kriging results.

4) Other tools. Installed ALDPAT tool, but couldn't get it to work because the program couldn't read the DEM.

In all cases, results in terms of the drainage network wasn't good, because we couldn't avoid the rectangular stream network; also, TWI still resulted in a blocks pattern.

In particular, with OK interpolated result, we got a point-like pattern DEM that affected the network result. However, the effect of the blocks pattern was diminished.

Also, we took a look at this question and answers...

Filtering out canopies and buildings from DSM to have a bare earth elevation

...which redirected us to Whitebox Geospatial Analysis Tools, but couldn't convert our DEM to LAS format. Also, we were not sure about the effectiveness of Bare-Earth DEM tool for us, because it's designed for removing semi-transparent objects, not blocks incorrectly "uplifted", which is our case.

We still want to generate a high quality DEM to make our hydrographic analysis, but don't know what else can we try.


texts: How to fill a symbolLayer with a custom texture?

texts:
I have this code :

pLayer is a QgsVectorLayer *

QgsFeatureRendererV2 * lrenderer = pLayer->rendererV2();
QList<QgsSymbolV2 *> lSymbols = lrenderer->symbols();
QgsSymbolLayerV2 * MySymbolLayer = lSymbols.at(0)->symbolLayer(0);
MySymbolLayer->setColor( QColor(90,90,90));
MySymbolLayer->setFillColor(QColor(90,90,90));

but how can I use the SetCustomTexture or setLineStyle from the 1.8 QGIS
Version , I want to fill the fill polygon with a specific texture or to define a type of line like DashLine

but I don't find the way to do it in C++ API  Qgis 2.8.2 for my application.


texts: How could I generate a CRS field inside the GeoJSON object created by Leaflet layer.toGeoJSON?

texts:
I'm using Leaflet to transform some drawn polygons to GeoJSON and send them to a backend. But the GeoJSON objects generated by layer.toGeoJSON are like this:

{
 "type":"Feature",
 "properties":{},
 "geometry":{
    "type":"Polygon",
    "coordinates":[[[18.69873046875,51.481382896100975],[18.69873046875,51.890053935216926],[20.456542968749996,51.890053935216926],[20.456542968749996,51.481382896100975],[18.69873046875,51.481382896100975]]]
  }
 }

When I receive that GeoJSON object in the backend of my application (Django based) and create a GEOS object from it, epsg:4326 is assumed. Even when epsg:3857 is used for the map.

I know GeoJSON spec supports CRS definition, but I don't know how to include it in the GeoJSON object generated by toGeoJSON.

Any suggestion?


texts: Combine multipolygons in feature collection while keeping boundaries in QGIS

texts:
This problem has been driving me absolutely crazy for a week now:

I have a map of a lot (1911) of districts. I need to add them all to a map, which is pretty heavy because each one requires a separate overlay, so I am trying to fit them all into one.

To do this, I was hoping to combine all the multipolygons in my feature collection (in a GeoJSON file). I've been trying to do all different kinds of merges, but they all end up erasing the boundaries, rather unsurprisingly.

Here is part of the map that I am working with:

Each of the districts above is a separate multi polygon.

This question asks something very similar to what I am asking, but it is about PostGIS, and I have no clue how to use PostGIS; from what I've read and attempted, there is a steep (and setup) learning curve.

How can I merge the features in a map while retaining their boundaries in QGIS?


texts: I would like to change the style of the component Flex Edition . The functionality is what I need but I find too much space it occupies such module for which I would like to make changes to style sheets


texts: Setting up rulebased labelling programmatically in QGIS?

texts:
In my plugin I do some labelling a vector layer with QgsPalLayerSettings. This works fine for me with QgsMapRenderer::setLabelingEngine.

But now I also want to implement the rulebased labelling. In API I found the class QgsRuleBasedLabeling::Rule which I thought works similar to setting up a QgsRuleBasedRendererV2::RuleList.

But I can't find the point where the QgsRuleBasedLabeling::Rule has to be set. Ther is no member in QgsMapRenderer. Does anybody has a hint?


texts: QGIS Processing GRASS7 output is not created in custom application on Mac | Linux

texts:
When trying to execute GRASS7s r.slope.aspect from QGIS using a custom application, no output is written on macOS. However, the algorithm runs fine and no error message is appearing. Doing the same on Windows runs just fine.

However, executing the algorithm in QGIS GUI works fine on macOS.

Input DEM: https://www.dropbox.com/s/mj7d42la95jlo3a/elevation.asc?dl=0

# export path variables
export DYLD_LIBRARY_PATH=/Applications/QGIS.app/Contents/MacOS/lib/:/Applications/QGIS.app/Contents/Frameworks/
export PYTHONPATH=/Applications/QGIS.app/Contents/Resources/python/
export QGIS_PREFIX_PATH=/Applications/QGIS.app/Contents/MacOS/
export PATH='/Applications/QGIS.app/Contents/MacOS/bin:$PATH'

# execute function
/usr/bin/python
import os
import sys
from qgis.core import *
from osgeo import ogr
from PyQt4.QtCore import *
from PyQt4.QtGui import *
from qgis.gui import *
QgsApplication.setPrefixPath(r'/Applications/QGIS.app/Contents', True)
app = QgsApplication([], True)
QgsApplication.initQgis()
sys.path.append(r'/Applications/QGIS.app/Contents/Resources/python/plugins')
from processing.core.Processing import Processing
Processing.initialize()
import processing
processing.runalg('grass7:r.slope.aspect', '/var/folders/5j/_1ts10x512sg_5q_3kzc_c1w0000gn/T/RtmpFnbVld/elevation.asc', '0', '0', 'True', '1.0', '0.0', '794599.107614635,798208.557614635,8931774.87460253,8935384.32460253', '0.0', '/var/folders/5j/_1ts10x512sg_5q_3kzc_c1w0000gn/T/RtmpFnbVld/slope.tif', 'None', 'None', 'None', 'None', 'None', 'None', 'None', 'None')


texts: Converting our Unicode to Persian in ArcGIS scripting by ArcPy

texts:
import arcpy
import os
from arcpy import env
env.workspace=r'I:\amozesh\samplegdb.gdb'
infeature=r'I:\amozesh\samplegdb.gdb\blocksclip'
cursor=arcpy.da.SearchCursor(infeature,['layer'])
cursor.next()

output
(u'\u0627\u0631\..........\u0631',)
How can I print it out using the Persian language?


texts: What is the smallest US geographic region with income data?

texts:
I'm looking for any region that is smaller than block groups. And it doesn't have to be for the whole United States. Does this exist? Does anyone have any ideas?

I feel like I saw somewhere that income data exists on the block level for some cities, but I haven't been able to find anything..


texts: Excluding pages without point within extent from Data Driven Pages?

texts:
ArcMap 10.2

I am creating a utility mapbook to show HVAC points on an installation. The data driven pages are named and sorted by a separate layer (area feature class). There are 70 named areas, however, not every page has HVAC units within the extent.

Is there a method to automatically exclude pages that do NOT include HVAC points within the viewing extent?

Notes:
- The data itself has to remain untouched so I can't add any additional fields.
- I'm familiar enough with python to run a script but I don't know what modules to call or run.


texts: Where can I find a Unesco World Heritage List as boundaries?

texts:
I was wondering if anyone knows where to find the UNESCO world heritage list with boundaries?  I can extract point data from the json feed for the map on the website http://whc.unesco.org/en/list via
 http://whc.unesco.org/?cid=31&l=en&mode=json but polygon or polyline data would be much more preferable.

Edit:  I was particularly interested in Australia, but it would be good to know where to go to find this sort of info worldwide, so if you have other countries, feel free to post them up as well.


texts: Storing related tables and toolboxes within feature dataset?

texts:
I am new to Enterprise geodatabases and am looking for some information regarding data structure and storage. Within my organization, we would like to store individual projects by feature dataset within a single Enterprise geodatabase. To my knowledge this is possible if the data does not contain any related tables or associated toolboxes.

Is it possible to store related tables and toolboxes within a feature dataset?

My goal is to not create an Enterprise geodatabase for each individual project.


texts: Setting correct options for displaying grid in QGIS composer?

texts:
I am using q. 2.5 4291337 in Windows 7 64 bits. I am working with raster and vector layers in EPSG::26714 (UTM NAD27 14N). In this version of QGIS, you are able to put a grid in a different CRS of the layers in Print composer. I tried to play with this option but I have some odd results.

When I set that I want a lat-lon Nad27 grid (EPSG::4267) and set the coordinate format to decimal, everything is OK. Longitudes are shown as negative numbers and latitude are positive.

When I set to show the labels to deg min sec with suffix the label shows wrong the hemisphere indicator. The suffix for my longitudes are labelled as "S" and for my latitudes are labelled "E".

I tried to replicate the problem with a project with layers in EPSG::4326 (Lat lon WGS84), and I do not have such a problem.

How do I set correctly the options to use a grid in a different CRS than my layers in print composer?

In the project with UTM layers, I chose the CRS for my grid to EPSG::4267. I set the spacing  as .01666666 for X and Y (which is equivalent to 1 minute).

Also I asked a friend of mine who uses Ubuntu and he does not have such a problem.


texts:
My colleague has several hundred xyz files (*.xyz) and wants to make a DEM to use it in ArcMap.
Files have three columns witohut first lane (variable names).

I suggested her:

to add first line with x, y and z column names (i know just manual way)
to convert xyz to dbf or similar
to add it in arcmap (as add x y data)
finally to use point to raster tool

I would probably know how to make a script for step 3 and 4, but i don't have any ideas how to solve step 1 and 2.

Or maybe there is an easier way (a tool,script) for direct conversion from xyz to dem?


texts: How to get all rivers from OSM using Overpass?

texts: i'm trying to get the files for Rivers, for most rivers i search for them in the Wizard using name=example and waterway=example and this works fine. For many rivers and canals overpass doesnt seem to be able to find the features and returns as map is left intentionally blank. I tried to use the selection box but then can't open in in QGIS all i need is the river does any one know of any way to get just the river file?


texts:
I've got some data that is in WGS84 UTM 29N and some that is in WGS84 Complex UTM 29N.

The WKT for the complex projection looks like this:

  PROJCS["WGS_1984_Complex_UTM_Zone_29N",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator_Complex"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-9.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0],AUTHORITY["ESRI",102579]]

I reprojected both sets of data using FME to EPSG:3857, and the Complex UTM 29 data ended up 1500 KM southeast of where it should have been. My guess is that FME didn't know how to handle PROJECTION["Transverse_Mercator_Complex"] because all of the other parameters are identical to the regular UTM 29 data.

Now I'm curious as to the difference, so I started searching around.

I searched the EPSG database and couldn't find a CRS match searching on "complex."
No good results on the first few pages of a Google search
Nothing obvious in the couple of map projection books I own

What is the Complex Transverse Mercator projection?


texts: How to decrypt/decompile old toolbox tbx-file in ArcGIS 10

texts: I have inherited an arcgis 10 toolbox tbx file for a web application. Are there any way to open the fil in for example Model Builder to see its interiors?


texts: Displaying point cloud footprint from PostGIS using GeoServer and OpenLayers?

texts:
I'm using PDAL to import LAS files into PostGIS.

The tutorial provided by boundless is great, and shows how to display the data in patches, but it doesn't cover displaying the table 'footprint'.

I need to display the coverage area as a footprint/wireframe using Geoserver and OpenLayers.

Does anyone know how to do this?


texts: Cannot remove background fill/highlight from labels generated through Easycustomlabeling

texts: In QGIS, I cannot remove the background fill or highlight that appears as a rectangle behind the text from the easycustomlabeling labels. The background fill is white with black text, and I want the black text only. I have tried clicking and unclicking the "highlight pinned labels" button at the top but it does not change anything.


texts:
I have a list of rasters in a folder that I want to project, rescale values and clip. My current code looks like this

library(rgdal)
library(maptools)
library(raster

setwd("D:/working")

# define projections
utm33s <- CRS('+proj=utm +zone=33 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0')

# read in and stack files
r.list <- list.files(pattern="*.tif", full.names=TRUE)
r.stack <- stack(r.list)

# project stack to utm 33s
r.stack.prj = projectRaster(r.stack, crs = utm33s)

# function to calculate float values
fun <- fun <- function(x) {
  y = (x - 50) / 200
  return(y)
}

# apply function on the stack
calc(r.stack.prj, fun=fun)

# read in polygon file for use in clipping
aoi = readOGR(dsn=".", layer="aoi_shape")

# clip raster using polygon
r.stack.prj <- crop(r.stack.prj, aoi, snap="out")

# save as GeoTIFF
writeRaster(r.stack.prj, 'clean_stack.tif', overwrite=TRUE)

The script crashes at 'projectRaster' What I'm I doing doing wrong?


texts: Layer SRS WKT: unknown, shapefile created with pyshp and python

texts:
The problem is solved with a proj file.
I created a shapefile from a csv file. I used ogrinfo and appears Layer SRS WKT: unknown.

Technologies

Python 2.7.9
pyshp compatible with Python 2.7
proj 4.8.0, 6 March 2012
geos 3.4.2
gdal
django 1.11 and GeoDjango

Sample of the csv file

Name;Latitude;Longitude;
ALAVA SALIDAS;42.842003;-2.740091
TIPSA FUENLABRADA LLEGADAS;40.266063;-3.848882

Python code

import shapefile as shp
import csv, sys

class convertCsvToShp:
def __init__(self, csv_path):
    self.csv_path=csv_path

# funtion to generate a .prj file
def getWKT_PRJ (self, epsg_code):
    import urllib
    wkt = urllib.urlopen("http://spatialreference.org/ref/epsg/{0}/prettywkt/".format(epsg_code))
    remove_spaces = wkt.read().replace(" ","")
    output = remove_spaces.replace("\n", "")
    return output

def code(self):

    try:
        # Set shapefile writer

        w = shp.Writer()
        w.autoBalance = 1

        # Fields are the city, the latitude and the longitude.
        w.field("name", "C", 50) # C are ASCII characters.
        w.field("latitude", "F", 10, 6)
        w.field("longitude", "F", 10, 6)

        # Readeing data from csv file
        with open(self.csv_path + 'coord.csv', 'r') as csvfile:
            reader = csv.reader(csvfile, delimiter=';')
            next (reader, None) # Skip the header of the csv

            for row in reader:
                # Create the point of the geometry
                w.point(float(row[1]), float(row[2]))
                # Attributes of the data
                w.record(row[0], row[1], row[2])

        w.save(self.csv_path + 'agencies')

        # create a projection file
        prj = open(self.csv_route + 'agencies.prj', 'w')
        epsg = self.getWKT_PRJ("4326")
        prj.write(epsg)
        prj.close()
        return True
    except:
        print "Unexpected error:", sys.exc_info()[0]
        raise

ogrinfo

ogrinfo -so csvtoshp/static/files/agencies.shp agencies

Problem

The problem is that openlayers do not show the coordinates correctly. I do not know that this Layers SRS WKT unknown is involved.
It is required a proj file? Yes, it is required it.

New ogrinfo


texts: Error when trying to visualize geography data with QGIS

texts:
I use PostGis, and want to use QGIS to visualize data. But I get error message when trying to show the points in the table.

I create a table:

CREATE TABLE poi
(
  id text NOT NULL,
  coord geography(point),
  CONSTRAINT poi_pkey PRIMARY KEY (id)
);

CREATE INDEX poi_coord ON poi USING GIST ( coord );

I add data to the table:

INSERT INTO poi VALUES ('P01',ST_GeogFromText('POINT(19.034561 47.486079)'));
INSERT INTO poi VALUES ('P02',ST_GeogFromText('POINT(19.030486 47.490477)'));
INSERT INTO poi VALUES ('P03',ST_GeogFromText('POINT(19.030328 47.49059)'));

I see the data in the table:

"P01";"0101000020E610000091B75CFDD8083340DD2230D637BE4740"
"P02";"0101000020E610000061FC34EECD073340583849F3C7BE4740"
"P03";"0101000020E610000031276893C307334050DF32A7CBBE4740"

I try to add the layer in QGIS , but I don't see anything, the log says the following messages :

Query: DECLARE qgis_1 BINARY CURSOR FOR SELECT st_asbinary("coord"::geometry,'NDR'),"id"::text FROM "public"."poi" WHERE "coord" && st_makeenvelope(-180,-90,180,90,4326) returned 7 [ERROR:  Antipodal (180 degrees long) edge detected!]

1 cursor states lost.
SQL: DECLARE qgis_1 BINARY CURSOR FOR SELECT st_asbinary("coord"::geometry,'NDR'),"id"::text FROM "public"."poi" WHERE "coord" && st_makeenvelope(-180,-90,180,90,4326)
Result: 7 (ERROR:  Antipodal (180 degrees long) edge detected!)

Query: CLOSE qgis_1 returned 7 [ERROR:  cursor "qgis_1" does not exist]

I tried it with QGIS 2.4.0 and with two DB configs:

PostgreSQL 9.1.7 with PostGIS 1.5.3
PostgreSQL 9.3.4 with PostGIS 2.1.3

When I use Geometry instead of Geography, it works well.

What did I do wrong?


texts: ArcPy not seeing ArcSDE datasets from IDLE or PyScripter?

texts:
I tested the code below from ArcCatalog 10.5's python console and IDLE+PyScripter. While in ArcCatalog everything works correctly(code is supposed to list available datasets and feature layers), in other python enviroments(IDLE and PyScripter) Datasets and Feature Classes are not shown, only views are listed. No error message is produced, as if my Datasets don't exist.

# Import arcpy module
import arcpy
from arcpy import env

# Local variables:
databaseConnection = r"C:\Users\UserName\AppData\Roaming\ESRI\Desktop10.5\ArcCatalog\SDENAME.sde"
arcpy.env.workspace = databaseConnection

datasets = arcpy.ListDatasets()
for dataset in datasets:
    print(dataset)


texts: Is there a geo-fence map of Apple's or Google's reverse geocoding?

texts: Is there a map that shows how lat/long maps to locality/county/sublocality/subadmin and possibly covered area?


texts: How to replicate ArcGIS Intersect in PostGIS

texts:
I'm trying to replicate this ArcGIS process in PostGIS: http://blogs.esri.com/esri/arcgis/2012/11/13/spaghetti_and_meatballs/. This describes how to break buffered points down into polygons based on their intersections, counting the number of layers, and attributing that to the polygons in order to classify them. I'm using it to create a rough point density map with vectors, and the results were surprisingly nice for my data set in ArcGIS. However, I am struggling to come up with something workable in PostGIS where I need it for producing dynamic point density layers for a web map.

In ArcGIS, I simply ran the Intersect tool on my buffered points layer to create the shapes I needed.

In PostGIS, I ran this query:

CREATE TABLE buffer_table AS
SELECT a.gid AS gid, ST_Buffer(a.geo,.003) AS geo
FROM public.pointTable a;

CREATE TABLE intersections AS
SELECT a.gid AS gid_a, b.gid AS gid_b, ST_Intersection(a.geo,b.geo) AS geo
FROM public.pointTable a, public.pointTable b
WHERE ST_Intersects(a.geo, b.geo) AND a.gid < b.gid;

DELETE FROM intersections WHERE id_a = id_b;

The output looks pretty much identical to the ArcGIS output, except that it is not breaking the polygons down to the same extent that is required for a meaningful density map. Here are screenshots of what I mean:

ArcGIS is on the left, and PostGIS is on the right. It is slightly difficult to tell, but the ArcGIS image shows the 'interior' polygon created where all 3 buffers intersect. The PostGIS output, on the other hand, does not create that interior polygon and instead it keeps its components intact. This makes it impossible to provide a classification for just that interior area with 3 layers on top of each other compared to just 1 for the outer parts.

Does anyone know of any PostGIS function to break the polygon down to the extent I need? Alternatively, does anyone know of a better way to produce a point density map with vectors in PostGIS?


texts:
We want to indicate a point on a map and like the idea of using a teardrop but it seems that google has the teardrop idea patented (http://www.google.com/patents/USD645052 http://www.google.com/patents/USD661704 http://www.google.com/patents/USD620950)  Does this mean that even if we create our own, unique teardrop shape that we would be violating the patents that google has on the tear drop shape?  Obviously it wouldn't look like googles, but it would look similiar since it would be a teardrop shape.  Any other ideas on a good point marker?  We are using Esri's tools, not google maps.

Thanks!


texts:
I am trying to set all values in my raster that are not equal to -9999 to Null using Arcpy.
The problem is that using != operator gets me:
arcgisscripting.ExecuteError: ERROR 010416: Error in setting raster table filter for F:\Data\test.tif.
However, other operators like = or < works okay.
Any ideas?

I solved the problem using <> instead of !=, but I still do not understand why ArcGIS Desktop and ArcPy use different operators for the same functions.
# Import arcpy module
import arcpy

# Set the current workspace
arcpy.env.workspace = "F:/Data"

# check out any necessary licenses
arcpy.CheckOutExtension("Spatial")

inRaster = arcpy.Raster('test.tif')
outSetNull = arcpy.sa.SetNull(inRaster, in_false_raster_or_constant=1, where_clause='Value != -9999')


texts:
Is there any way to add parallel offset for those lines which oneway road are False?
All line here are straight (segment). So We just have two points for each line (start and end points).
# There is 43 857 edges with 14 561 oneway roads and 29 296 bi-
  directionals roads

G = nx.from_pandas_edgelist(edges,'origin', 'destination',
                              ['key','id_edge', 'name', 'lanes','length', 'speed',
                              'capacity', 'function', 'geometry'],
                              create_using=nx.MultiDiGraph)

# Add nodes to the graph G
  data = nodes.set_index('id_node').to_dict('index').items()
  G.add_nodes_from(data)

# Add oneway column (it will be a boolean)

  for u,v,d in G.edges(keys=False, data=True):
      if G.has_edge(v,u):
          G.edges[u,v,0]['oneway'] = False
      else:
          G.edges[u,v,0]['oneway'] = True

m = folium.Map(location=[48.86088, 2.33765], zom_start=2)
folium.features.GeoJson(edges).add_to(m)

To have parallel side, I tried this :
points = []
for u,v,data in G.edges(keys=False, data=True):
    if 'geometry' in data:
        xs,ys = data['geometry'].xy
        point = list(zip(xs,ys))

        if not data['oneway']:
            points.append(point)

gdf = gpd.GeoDataFrame({'geometry': [LineString(l) for l in points] })
gdf.set_geometry('geometry', inplace=True) # set the geometry
gdf.set_crs('EPSG:4326', inplace=True) # set the actual CRS

gdf1 = gdf['geometry'].apply(lambda x:x.parallel_offset(0.6, 'left', 2)

m = folium.Map(location=[48.86088, 2.33765], zom_start=2)
folium.features.GeoJson(edges).add_to(m)
folium.Choropleth(gdf1,line_weight=1,line_color='red').add_to(m)

But I got strange lines everywhere


texts: Export Image in .CSV on Google Earth Engine

texts:
This is what I'm dealing with:
I have three filtered image collections concatenated in one (Sentinel-1, Sentinel-2 and SRTM), and I would like to export this in CSV for external analysis.
The SRTM dataset generates 2 bands (aspect and slope);
The Sentinel-1 filtered collection has 158 images and 2 bands;
The Sentinel-1 filtered collection has 300 images and 10 bands;
Its a big data cube. I've been trying image.toArray(0).toArray(1) and then arrayFlatten, but its not working.
Any ideas?


texts: Using What3Words without tight-coupling to "computer tech"?

texts:
This question is a result of seeing a headline about Mercedes Benz adding What2Words to its in-dash navigation system. I think I understand what What3Words is about but I can't help but feel I'm missing something about how & why I'd want to use it.

In a nutshell, it places a roughly 3m x 3m grid over the planet. Each "square" is given a three word identifier. I give people the three words as my location. Great! Now that we have a basic understanding, let's meet for pizza and beer at whiskey.tango.foxtrot around 6pm.

Without using a computer, how do I find whiskey.tango.foxtrot? A street address (123 main street, anytown, USA) doesn't require any electronic tech and it even works when the power grid goes down. If I'm across the city but on the same street, I can look at the numbers and just walk. I might be 50 yards from cat.dog.chicken but how does that relate to my current position of donkey.cat.moose? I know that 250 main street, anytown, USA is fairly close to 123 main street.

A latitude / longitude coordinate pair gives me absolute location and I can estimate distance between multiple pairs. Think maps, sextants, stars, etc - you can locate yourself on the globe. Think USGS topo quad sheets. Latitude increase as you move away from the equator. Longitude increase as you move away from the Prime Meridian.

W3W seems to use a proprietary algorithm and requires using its service to decode W3W. You don't know what word is adjacent to another word and unlike lat/lon, you cannot determine direction from two words.

I understand W3W is easier to remember than addresses or coordinates but, in my opinion, the wheels quickly fall off when the technology (internet connectivity, cell service, electricity, a computer) is taken away.

I've considered using W3W on a project but determined the ease of remembering it -seemingly a huge selling point - was outweighed by actually being able to use the values.


texts: How to style a Vector layer in QGIS using a numeral of an attribute

texts:
I have a simple vector layer in QGIS that I would like to style with numerals based on the ID attribute.  For instance, point with ID 1 will show up as the numeral 1.  Is this possible with default libraries?

wkt_geom             ID
POINT(23.85  37.99)  1
POINT(23.8   37.9)   2
POINT(23.855 37)     3
POINT(23     37.9)   4


texts:
Solved.

In QGIS I turned off CRS. It then read the data of the image correctly.

I have brought a tiff photo of a side of a building into qGIS. It is about 300 images that were put together in Pix4d to an orthomosaic. When I took the photo I marked off 5m segments on the building with tape.

When I use the "Measure" tool in QGIS between the tape marks I get about 553,000 meters.

Is there an easy way to calibrate the measure tool to a known distance in the photo? I'd like to calibrate it and then use the measure tool to measure other things.

The reason I'm doing this is that I need to share the photo with someone, and I'd like them to be able to do measurements on the photo. QGIS seemed an ok way to do this.


texts: How to perform shapefile/raster integrated functions

texts:
Description:

I'm having trouble with some Raster/Shapefile manipulation. My problem is that I have a shapefile that I would like to use to update my rastervalues with. The shapefile has multiple polygons (Countries), that I would like to to loop through and change the value of the cells that are masked by that shapefile. Most other questions are about clipping a raster with shapefile, but I just want the cells on the original raster to be masked.
Here is an example of what I'm looking for:

for shapes in shapefile:
    masked_country=maskfunction(shapes,originalraster)
    masked_country=masked_country + arbitraryfunction(shape['layer'])

Do you have anything that works like that?  I work with python and gdal.

Here are some of the other options I have tried:

Trying to rasterize shapefile using python gdal.RasterizeLayer:

I only get a black tiff using this method. I have even implemented the target_dst=None at the end, and I'm still getting a 0 value raster. Here is my code sample for this:

from osgeo import gdal #geo library
from osgeo import ogr

pathtofile1="rPath\to\my\shapefile"
pathtofile2="rPath\to\my\raster"

shapefile = ogr.Open(pathtofile1)
rasterfile=gdal.Open(pathtofile2)

mb_l = shapefile.GetLayer()
geo_transform = rasterfile.GetGeoTransform()
source_layer = rasterfile.GetLayer()

x_min = geo_transform[0]
y_max = geo_transform[3]
x_max = x_min + geo_transform[1] * rasterfile.RasterXSize
y_min = y_max + geo_transform[5] * rasterfile.RasterYSize
x_res = lur250.RasterXSize
y_res = lur250.RasterYSize
pixel_width = geo_transform[1]

target_ds = gdal.GetDriverByName('GTiff').Create('my.tiff', x_res, y_res, 1, gdal.GDT_Float32)
target_ds.SetGeoTransform((x_min, pixel_width, 0, y_max, 0, 0-pixel_width))

band = target_ds.GetRasterBand(1)
NoData_value = -999
band.SetNoDataValue(NoData_value)
band.FlushCache()

gdal.RasterizeLayer(target_ds, [1], mb_l,options=["ATTRIBUTE=OBJECTID"])

target_ds=None

There is a possibility that the "OBJECTID" is not Int and it's string, but I don't know how to check that.

Also I have played with the arguments in SetGeoTransform, but to no use (Setting Ymax or Ymin for different raaster definition).

Trying to burn the values of the shapefile into a raster with similar projection and resolution using gdal_rasterize from bash:

I need to create an empty raster with the same coordinate as my rasterfile (which can be achieved by copying that raster), but I don't know how to set the value for the entire raster to "0" or NoDataValue. All other solutions for these kind of problems are using "gdal_calc.py" and our HPC gdal installation is without such a tool.

Extra information

Here is the gdal_info results for my raster:

Driver: GTiff/GeoTIFF
Files: final4.tif
Size is 141969, 60829
Coordinate System is:
PROJCS["World_Mollweide",
    GEOGCS["GCS_WGS_1984",
        DATUM["D_WGS_1984",
            SPHEROID["WGS_1984",6378137.0,298.257223563]],
        PRIMEM["Greenwich",0.0],
        UNIT["Degree",0.017453292519943295]],
    PROJECTION["Mollweide"],
    PARAMETER["False_Easting",0.0],
    PARAMETER["False_Northing",0.0],
    PARAMETER["Central_Meridian",0.0],
    UNIT["Meter",1.0]]
Origin = (-17619594.546999998390675,8722279.461999999359250)
Pixel Size = (250.000000000000000,-250.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-17619594.547, 8722279.462) ( 30d 5'16.03"E, 83d 8'14.83"N)
Lower Left  (-17619594.547,-6484970.538) (107d 3'58.73"E, 55d59' 2.97"S)
Upper Right (17872655.453, 8722279.462) ( 20d10'44.20"W, 83d 8'14.83"N)
Lower Right (17872655.453,-6484970.538) (103d26' 0.82"W, 55d59' 2.97"S)
Center      (  126530.453, 1118654.462) (  1d16'20.33"E,  9d 3'42.49"N)
Band 1 Block=141969x1 Type=Byte, ColorInterp=Gray
  NoData Value=-128


texts:
I want to merge 10 las files into one, and I am having problems with lasmerge of lastools. It isn't able to deal with such a high number of points (average of points per file almost 4 million of points).

Does anybody know other way to merge multiple las files into one? Other softwares or perhaps some python code?

I realised that there was not an error, simply the lasmerge put the new merged file into the bin folder, not in the specified folder.


texts: I have an NDVI for a very large region, so I clipped it so it would fit my study area. The problem I am encountering however is that when I apply the clipping function the pixels come out differently than before. I think this is due to the min and max value arcGIS sets but I have no idea how to change it. The image show what I mean, the red outline is study area.In the first image it shows the entire NDVI but when i clip it the values change greatly.


texts:
i'm using qgis 2.0 and i want to calculate a viewshed using r.los from grass tool box.
after typing in the python console and a few seconds, i don't get an error message, but there is no file.

the r.los gui works fine.

processing.runalg('grass:r.los',"Elevation",'642229,5598076',None,75,5000,False,None,None,'C:\\r_LOS')
{'output': 'C:\\r_LOS.tif'}

someone got an idea?!

thanks - niels...


texts: Export to shapefile to csv from spatial join code

texts:
I have created some code to perform a a one to many spatial keeping all rows, I would like to export the resulting shapefile to CSV as that is the required format for the program the the data will be used in?

Spatial Join code:

import arcpy
from arcpy import env

env.workspace = r"U:\GIO\Geospatial_Services\ANALYTICS\Tools\Parliamentary"
env.overwrite = True

Parcel = r"U:\GIO\Geospatial_Services\ANALYTICS\Tools\Parliamentary\MP_Cut_Many.shp"
int_shp = r"U:\GIO\Geospatial_Services\ANALYTICS\Tools\Parliamentary\MP.shp"

output_speed_fc = r"U:\GIO\Geospatial_Services\ANALYTICS\Tools\Parliamentary\MP_Join_2.shp"

arcpy.SpatialJoin_analysis(Parcel, int_shp, output_speed_fc, join_operation="JOIN_ONE_TO_MANY", join_type="KEEP_COMMON", match_option="INTERSECT")

print 'worked'

The CSV conversions I have tried before haven't worked, got to convert to DBF, I am using Arcgis10.1


texts:
I'm doing a Mann-Whitney U test to compare two very large samples (sample size 1 = 13250; samlple size 2 = 38871) originating from a shapefile with 52121 features. I know t-tests are not recommended to compare rasters, because since rasters have so many values, they will almost surely detect a significant difference, no matter how small that difference might be (see this post in gis.stackexchange.com and point 3 of the first answer to this post).

I ran the test on SPSS and got the following results:

Test Statistics
Mann-Whitney U    4094520,000
Wilcoxon W    759591276,000
Z -210,227
Asymp. Sig. (2-tailed)    ,000

While I did expect to find a difference between the groups, I don't know what to make of such large values.

There is a paper which did the same thing I did. See the last sentence of the Material and Methods section. However, they obtained much smaller values http://icesjms.oxfordjournals.org/content/69/1/75/T5.expansion.html (they don't make clear if these are Mann-Whitney Us or z-scores). Granted they used a smaller sample (combined sample size = 5629), but the difference in the magnitude of the values still seem strange to me.

So, are my results simply the result of a very large sample, but still valid? Or should I use another test?


texts: Why is raster boundary rectangular, not buffer?

texts: First, I took very few points and did IDW and i got a circular boundary raster. This is because of the search radius i applied. If that is the case, for large number of points, the raster must look like a buffer, with the buffer width of search radius. Instead, the edges are rectangular. Why is it so? I didn't find an answer in any previous post.


texts: Dissolving contiguous clusters of polygons and identifying them separately in attribute table of ArcGIS for Desktop?

texts:
I have a layer comprised of about 500,000+ simplified polygons that resulted from converting raster cells to polygons with raster cell values assigned to them. When I select the resultant polygons by attributes (say, those that have values above a certain threshold), I end up with a bunch of contiguous spatial clusters. My goal is to create an outer shell boundary circumscribing those clusters, which can be achieved using the Dissolve Tool in ArcMap.

The problematic twist here is that I want to be able to identify these clusters of dozens of small polygons individually after dissolution into "island" shapes. How would I go about it?


texts: How to update a raster layer's statistics (min/max) in PyQGIS after changing its data?

texts:
I want to alter a raster layer with Python code and update the style afterwards. The style part is out of reach because I encountered the following problem:
I added a Float64 raster layer to my project using the "Create constant raster layer" processing tool. The layer is in my project's table of contents, named "Constant".
Then I run this script in the Python console:
import struct

layer = QgsProject.instance().mapLayersByName('Constant')[0]
data_provider = layer.constDataProvider()
block = data_provider.block(1, layer.extent(), layer.width(), layer.height())

new_data = struct.pack("dd", *[-99, 99])  # two double values
block.setData(new_data)

data_provider.setEditable(enabled=True)
data_provider.writeBlock(block, 1)
data_provider.setEditable(enabled=False)

The TIF file is rewritten now. If I open it manually, I see the new values in the first two pixels in the top left of the image.
But the "Constant" layer that was already in my ToC and on which I operated, still shows the same old constant value as its min and max values. I tried setting the "Min / Max Value Settings"' "Accuracy" to "Actual (slower)" with no success. I also tried data_provider.reload() and layer.triggerRepaint() without success.

I do not want to set the min and max values manually, but have QGIS execute whatever it does when I set the "Min / Max Value Settings"' "Accuracy" to "Actual (slower)".
How do I make QGIS recognize the new data in the layer and update its statistics?


texts: OpenLayers - is it possible for a feature to have both a single click and double click event?

texts:
I'm have an OpenLayers map with a single vector layer containing many features. I've been using SelectFeature to pop up some info when a feature is clicked, and that's working fine. However, I now want to add some functionality when the feature is double-clicked - specifically I just want to zoom in, and not have the popup appear.

I gave this a shot: http://blog.slashpoundbang.com/post/2342888467/openlayers-zoom-when-a-feature-in-a-vector-layer-is
But it removed my single-click functionality.

Essentially this is the same problem as SelectFeature with Double-Click

Any ideas?

UPDATE:

So I'm trying a workaround using SelectFeature to handle the single click and registering an event listener on the vector layer to handle the double click. All I want to do is zoom in on the double clicked feature. This is what I have:

VEClayer.events.register('dblclick', VEClayer, function(e){
    map.zoomToExtent(e.features.geometry.getBounds(), closest=true);
    }
);

However now I have another problem, and I'm not sure if it's being caused by having both SelectFeature and the listener active or what. The new issue is that try as I might I can't access any of the event attributes. So when I try to get e.features.geometry.getBounds(), Firebug gives me the error "e.features is undefined". I've tried to get this multiple other ways (e.xy, e.longlat, getMousePosition) and it's always the same result.


texts: CQL to XML with string substitution using geotools?

texts: I have some cql attr > PLACEHOLDER that I want to convert to valid XML using geotools. I have heard of terms like filter functions that can somehow replace the placeholder with a function, but not totally sure the best way to go about this.


texts: I have stacked layer (from ERDAS) of Landsat 1-5 and 7, bands. I want each of the bands in a separate layer. How can I do that?


texts:
I'm learning the Manual Photogrammetry process, Stereoscopy, to Process UAV Imagery.
I'd like to know about others software, especially open source or low cost.

ERDAS Imagine, Erdas Extensions for ArcGIS.
Menci StereoTools.
Digi3d.
E-foto (I tried this but without success, If you like to help me Photogrametry using e-Foto (Stereo Plotter, Stereo Imagery, Stereoscopy))
EnsoMOSAIC.


texts: Is 64-bit background geoprocessing native to ArcGIS 10.2?

texts: I cannot seem to find any documentation as to whether or not 64-bit background geoprocessing is native to ArcGIS 10.2.  Or, is a separate download and install required as in 10.1?  There were rumors prior to the 10.2 release that 64-bit background geoprocessing comes pre-installed with the upgrade/download.  Is there an authoritative reference on this subject?  No relevant system information is available on my Arc Administrator configuration information.


texts: Is it possible to host MapServer in the cloud, specifically in Azure?

texts:
Can MapServer be hosted in the cloud?

I presume as a virtual machine, or is it possible to have it as an Azure Web App, and avoid the hassle of VM maintenance

And by cloud, I mean Microsoft ... but if I put Azure, my question gets flagged as a duplicate of my question about hosting GeoServer in Azure


texts:
I have a series of LineString that I want to have the length in meters. I use the shapely.geometry module, and importing LineString.

I've seen this post, explaining how to do it, but it seems using Django : Get LineString length in meters (Python, GEODjango)

# Create a line string
>>> line = LineString([(0, 0), (360, 0)])

# Specify srid using [WGS84][1]
>>> line.srid = 4326

# Transform into projected coordinate system (using web mercator)
>>> line.transform(3857)

# Line length. This is zero within numerical precision, because the start
# point (0, 0) is the same as the end point (0, 360)) on a map.
>>> line.length
>>> 1.1329847282581795e-08

# An example forcing the line through the other side of the world,
# and specifying srid on definition.
>>> line = LineString([(0, 0), (180, 0), (360, 0)], srid=4326)
>>> line.tansform(3857)
>>> line.length
>>> 40075016.6855784

I try to use this post, but it makes me the following error :

'LineString' object has no attribute 'lenght'

IS there any module to import?


texts: Including join features shape geometry in output feature class using Spatial Join in ArcPy?

texts:
I am trying to join two feature classes together and also keep the spatial geometry of the join features. I want to keep the target and output geometry the same but also have a field like JoinSHAPE where I can access the geometry of the feature that it was joined with.

I have tried to add a field mapping for the join feature but an error occurs.

RuntimeError: FieldMap: Error in adding input field to field map

when I run the code below.

import arcpy

target = arcpy.GetParameter(0)
join_features = arcpy.GetParameter(1)
output =  r"S:\GIS TOOLS\Workbench.gdb\asdf"

fieldmappings = arcpy.FieldMappings()
fieldmappings.addTable(target)
fieldmappings.addTable(join_features)

shape = arcpy.FieldMap()
shape.addInputField(join_features, "SHAPE@")
fieldmappings.addFieldMap(shape)

arcpy.SpatialJoin_analysis(target, join_features, output, "#", "#", fieldmappings)


texts:
I am working with an excel/csv file which contain Easting Northing & Elevation value (Its all about point feature). But when i am converting this file in dwg by GIS for use in auto cad then the elevation value is not showing in cad only point feature is showing.

How would I resolve the problem?


texts:
I have many points that were picked by  handheld GPS and I want to fit them to their relevant cadastral plots automatically.

The attribute information (plot number) for these points are exactly similar to the attribute information (plot number) for the cadastral plots.

However, I am unable to automate the process of moving the many points which are outside their relevant plots. Instead I am forced to do them manually.

I have tried spatial join but the points can't be moved to their plots.

Any idea on how to do this in ArcGIS 9.3?


texts: Measure collected polygons in Geometry Generator, avoiding extra texts from multipart polygons

texts:
I combined a set of polygons on the first instance of Geometry generator using something like this:

buffer(buffer(collect( $geometry, concat("MZNA","REGLABEL")), 0.1), -0.1)

Buffers are meant to join poligons with a little separation. Thing is, after that I defined marker lines to simulate measurement, but in the case of multipart polygons, I get random measures floating between the different 'islands'. So I want to separate them before measuring, but I dont know how...

EDIT:

Im using the approach shownin QGIS How to label individual sides of polygon by @she_weeds. My problem is, as shown in the pic, after collecting and buffering, I get some multipolygons and those ones gets an extra text floating between the polygon parts. Is there a way to avoid it without dissolving the generated geometry?

Im not experienced enough to know if its possible (started using QGis 2 weeks ago), but i thought some options:

Dissolve polygons instead of collecting them, using a python function. Didnt know how to call dissolve from the function editor, didn't even know if it was possible.
Apply a filter to the collect statement.
When generating lines from points, or text from those lines, somehow test the points to condition the generation of lines.
When generating points from nodes, test multipart polygons and generate different sets of points. Managed to successfully test that but didnt manage to create more than one set of points.


texts: Is it possible to embed a QGIS map with all the information I added, in my web site?

texts:
I am not a cartographer or geography student. I work in a network organization for developmental policy and I need a good software which enables me to record the cities in which our members work on a worldwide map. The problem is, that I must be able to embed that map in our website.

Do you know if that is possible with a map from QGIS? And if it's not possible, do you know a software which I could use for that purpose?


texts: std:exception error in QGIS when importing .dwg

texts: I have a .dwg file that when I try and import it into QGIS I get a std::exception error. I have tried converting the file to a .shp file and then loading it but when I do this it doesn't seem to load in the correct location. How can I either georeference the converted .shp file or import the .dwg file into a map?


texts: How to set the old .qgs file as the default project format in QGIS 3?

texts:
I see this discussion that one can change the default QGIS project format back to the text based .qgs. But there is no detail. Quoting the discussion:

  There is an option to start QGIS with another default setting file

  --globalsettingsfile

  https://docs.qgis.org/testing/en/docs/user_manual/introduction/qgis_configuration.html#running-qgis-with-advanced-settings

  You can have a look at how the french ministry of ecology handles this
  issue here:
  http://www.geoinformations.developpement-durable.gouv.fr/qgis-package-3-4-2-pour-test-a3662.html

But the linked page does not mention .qgs.

Does any one have a complete answer on how to change QGIS 3 default to use the old .qgs format?

(as of QGIS 3.12, Ubuntu 18.04).


texts: Extracting attributes from WMS layers in OpenLayers

texts:
I am serving up a WMS layer from GeoServer in OpenLayers and am trying to extract a specific attribute from an object in this layer when I click on the map.

I am defining my WMS layer like this:

       oas = new OpenLayers.Layer.WMS(
       "oas",
       "http://www.bwgeoserver.com/geoserver/UK_Admin/wms",
        {
          "LAYERS": 'wealthmap_bdy',
          transparent: 'true',
          extractAttributes: 'true'
        },
        {
          buffer: 0,
          opacity: 1,
          displayOutsideMaxExtent: true,
          isBaseLayer: false
        }
        );

map.addLayer(oas);

The layer is displaying OK on the map. I am using a proxy server and there isn't a cross-domain problem as I can select and display all of the attributes in a popup.
"OA11Code" is definitely a field on the WMS layer, and it displays OK with the other attributes in the popup.

But when I try to extract this attribute from the selected object using this code:

oaPicker = new OpenLayers.Control.WMSGetFeatureInfo({
             url: 'http://www.bwgeoserver.com/geoserver/UK_Admin/wms',
             title: 'identify features on click',
             layers: [oas],
             queryVisible: true
         });

   oaPicker.events.register('getfeatureinfo', this, selectid);
   map.addControl(oaPicker);
   oaPicker.activate();

   function selectid(e) {
     var val = e.features[0].attributes.OA11Code;
   }

I get an error in the Firefox debugger saying: "TypeError: e.features[0] is undefined".

All I want is to be able to click on the map, extract the value of an attribute from the WMS layer, assign it to a JavaScript variable and do something with it.


texts: gdal_grid unable to get nodata outside of original data area

texts:
I have a vector based csv file that I want to convert into a grid using gdal_grid, but I want to force null data values strictly around the area where there is data.  I can get gdal_grid to work but it is always just smoothing the whole grid rather than giving me null data outside of the shape. Can someone tell me what settings to use to make sure that I get null data inside the grid where there is no data within the original csv.

Many thanks


texts: Saving custom labelling (template) so it can applied to other feature classes?

texts:
I sometimes work with large NTDB datasets, and was just wondering the best way to go about saving all of my custom labelling.

Is there a way I can save all of my custom label classes so I can apply them to all of my NTDB datasets so all maps made with this dataset are standardized?

I can do it the long way by changing the classes every time, but just wanted to know if I could save it so I can just apply it i.e. like a layer file?


texts: How do I have logic in ESRI's popup configs for the Flex api/viewer?

texts:
The popup window config functionality I'm referencing is this... http://help.arcgis.com/en/webapps/flexviewer/help/index.html#//01m30000002q000000

This is my config...

<?xml version="1.0" ?>
<configuration>
    <title>Wind</title>
    <description>
        <![CDATA[
            Forecast is valid for {DATA_DATE} UTC
            Wind Direction: {WIND_DIRECTION} degrees
            Wind Speed: {WIND_SPEED} knots
            Wind Chill: {WIND_CHILL} F
        ]]>
    </description>
    <fields>
        <field name="DATA_DATE">
            <format useutc="true" />
        </field>
        <field name="WIND_DIRECTION"/>
        <field name="WIND_SPEED"/>
        <field name="WIND_CHILL"/>
    </fields>
</configuration>

The problem is that sometimes I want to display N/A for the Wind Chill value. This only occurs when wind chill is -999. Is it possible to display something different than what is in the layer data in this config driven popup? If so, how?

I'd rather not have the column for wind chill be a varchar and have the string "N/A" in it.


texts: Bulk inserting with Python and cx_Oracle into st_geometry enabled table?

texts:
I'm trying to make a bulk insert of a CSV file with Python and cx_Oracle into an st_geometry enabled table of an Oracle database with ArcSDE.

I could do it with the method execute() of cx_Oracle but for big files is not the faster approach. That's why I'm trying to make the method executemany() work. With executemany() I'm able to insert registers with strings, integers,... but not st_geometry fields. With st_geometry fields I get the error:
     cx_Oracle.DatabaseError ORA-01036: illegal variable name/number

With that code "print cursor.bindnames()" returns ['1','2'], what it means that :3 y :4 are not recognized.
Just for the sake of the argument, if I modify st_geometry('point(:3 :4)',0) to st_geometry(point(:3 :4),0) "print cursor.bindnames()" returns ['1','2','3'], what it means that :4 is not recognized.

This is the sample code I'm using right now:

Workaround found, see next answer

# Bulk insert
#
# Sample test st_geometry table
# create table test (col1 varchar2(30 byte), col2 varchar2(30 byte), shape sde.st_geometry);
#
# Sample sql insert to insert data into the table
# insert into test (col1, col2, shape) values ('val01', 'val02', sde.st_geometry('point(12.67 2.32)', 0));
#
# Sample file.csv
# val11,val12,34.22,23.22
# val21,val22,12.33,14.22
# val31,val32,13.21,18.32
#

import cx_Oracle

dsn_tns = cx_Oracle.makedsn('server', 1521, 'sid')
dbCon = cx_Oracle.connect('user', 'password', dsn_tns)
cursor = dbCon.cursor()

list = []
f = open('file.csv','r')

# srid must be a valid code in sde.st_spatial_references
cursor.prepare("""INSERT INTO test (col1,col2,shape) VALUES (:1,:2,sde.st_geometry('point(:3 :4)',0))""")
for line in f:
    pair = line.split(',')
    list.append((pair[0],pair[1],float(pair[2]),float(pair[3])))

print cursor.bindnames()
cursor.executemany(None, list)
dbCon.commit()
f.close()
cursor.close()
dbCon.close()


texts: How to write CRS info from QGSVectorLayer

texts:
I'm writing a plugin which will create a new virtual layer containing points. I know the CRS of the new layer, as it is the same as the layer I am using for input. I would like to specify the CRS in code but the CRS selector dialog pops up when I create the layer. I would like to suppress that dialog so that something like this will work.

def createVirtualLayer(self):
    self.vl = QgsVectorLayer("Point","block_labels", "memory")
    self.vl.setCrs(self.layer.crs())


texts: Uploading TAB data to MSSQL using ogr2ogr

texts:
I have been trying various methods to import a shapefile into MSSQL but keep getting

FAILURE: no source datasource provided.

ATTEMPTS

ogr2ogr -f "MSSQL:Server=xxx;Database=_GIS_Test;UID=yyy;PWD=zzz"  -nln “D:\QGIS_Development\mssql\bridges.shp” “MGM_Bridges_test” a_srs EPSG:4326 -t_srs EPSG:28354 -nln “MGM_Bridges_test”
CD to folder and

ogr2ogr -overwrite -f “MSSQLSpatial” “MSSQL:Server=z;Database=_GIS_Test;UID=xxx;PWD=yyy"  bridges.shp a_srs EPSG:4326 -t_srs EPSG:28354 -nln “MGM_Bridges_test”

3.

ogr2ogr -overwrite -f “MSSQLSpatial” “MSSQL:Server=z;Database=_GIS_Test;UID=xxx;PWD=yyy" a_srs EPSG:4326 -t_srs EPSG:28354 -nln bridges.shp “MGM_Bridges_test”

ogrinfo on the input dataset opens it.

I have checked many posts here but no one reports this issue. Tried bridges.shp in quotes. I have removed the "“MGM_Bridges_test”" and I get the error No Target Datasource provided.

---- UPDATE

I tried typing the following but still get the same error...

ANOTHER attempt -same issue


texts: Points contained in envelopes that cross date-line - java

texts:
I'm using the ESRI Geometry Java API https://github.com/Esri/geometry-api-java

I want to construct an envelope and see if it contains certain points. It's possible the envelope crosses the International date line.

Here is some sample code I've tried:

    double lon = Double.parseDouble("179.99");
    double lat = Double.parseDouble("0");
    Point pt = new Point(lon, lat);

    double minlon = Double.parseDouble("170");
    double minlat = Double.parseDouble("-20");

    double maxlon = Double.parseDouble("-170");
    double maxlat = Double.parseDouble("20");

    Envelope ev = new Envelope(minlon, minlat, maxlon, maxlat);

    SpatialReference sr = SpatialReference.create(4326); // decimal degrees

    boolean test = GeometryEngine.contains(ev, pt, sr);

    System.out.println(GeometryEngine.geometryToJson(sr, ev));
    assertTrue("envelope containts point true", test);

The test fails, and if I check the print statement I get a json that has swapped my x-coordinates

{"xmin":-170,"ymin":-20,"xmax":170,"ymax":20,"spatialReference":{"wkid":4326}}

Is there a different/better way to do this?


texts: Creating a "Stream Order" of canal system using ArcMap

texts:
I'm interested in creating a "Stream Order" for our canal system. I am the GIS Specialist for an irrigation district. I would like to have the ability to have a layer for when I click on a section of canal, it traces/highlights upstream all the way to our source of surface water. All of our water comes from the same source to the east, and disperses among hundreds of miles of canal to cover our district.

What tool would I use to do that in ArcMap?


texts: Median image of merged image collections

texts:
Following this example I have merged 2 image collections to create a synthetic twice daily Aqua and Terra composite. Now I want to generate the median image of this combined collection. For a single collection this already works with:

var collection = ee.ImageCollection('MODIS/006/MOD09GQ');
var filtered_date = collection.filterDate('2010-10-01', '2010-12-10');
var medianImage = filtered_date.median();
Map.addLayer(medianImage, {bands:['sur_refl_b02'], min:100, max:4500});

However, the merge seems to interfere with the following script.

var collection_1 = ee.ImageCollection('MODIS/006/MOD09GQ');
var collection_2 = ee.ImageCollection('MODIS/006/MYD09GQ');
var collection_date_1 = collection_1.filterDate('2010-10-01', '2010-12-10');
var collection_date_2 = collection_2.filterDate('2010-10-01', '2010-12-10');
print('Terra collection', collection_date_1);
print('Aqua collection', collection_date_2);

// merge collections
var collection_merge = collection_date_1.merge(collection_date_2);
print('Merged collection',collection_merge);

// sort merged collection
var collection_merge_sorted = collection_merge.sort("system:time_start");
print('Sorted merged collection', collection_merge_sorted)

var medianImage_10 = collection_merge_sorted.median();
Map.addLayer(medianImage_10, {bands:['sur_refl_b02'], min:100, max:4000});

The script returns that 'collection_merge_sorted.median is not a function', which seems to be caused by the merge itself. How to avoid this?


texts: Routing problem finding shortest route to line

texts:
inspiration on how to solve this problem.
(Ascii art, I don't have points to insert picture)

I need to find Z

Z = Shortest route from X Using Road (no point exists here)

O = The points that makes the line

X = Start point

Y = Closest point on line (Can't be used)

The problem right now is:
The nearest point (Y) from X will give a long route.

Route from X to Points (O) on line also returns a long route.

Any idea on how I can find Z ??

    O(Point)
    | ||
Line| ||       O2----------O2 (closer line but no road to this line so this can't be used)
    | ||
    | ||
    |Y||         X (Start point)
    | ||         ||
    | ||    ||   ||Road
    | ||    ||   ||
    | ||    ||   ||
    |Z||----||----Road
    | ||    ||
    | ||    ||------------- Road
    | ||
    | ||Road
    | ||
    | ||
    O (Point)

Is there anyway I can take the line and add more points to it on every 5 meters? So the Line would look like O-O-O-O-O-O instead of O---------O

Setup:
I have a postgres 9.3 with postgis 2
The database has all X (points)
All Lines (linestring)
And all roads (Export from openstreet)

(can't get pgrouting to work right now)
And a local osrm server running.


texts: How to get imported data from excel in the shapefile properties?

texts:
I imported some data from excel (as csv) to ArcMap. The information from the "excel table" has been joined to my shapes and is shown at the attribute table of every shape/layer. (As you can see below).

However, when I export it to GIS online I see the information isn't there. With some research I found out the information is visible at the attribute tables, but it isn't at the fields in the shapefile properties.

Is there a way to add the 2006,2009,2010,Total fields from the attribute table to the shapefile properties?


texts: Maintaining coordinates when converting SpatialPolygonsDataFrame to simple feature in R

texts:
I am ultimately trying to determine the area of a home range, but I need to remove areas that overlap land. To accomplish this, I have calculated and plotted my home range, and want to use one polygon (land) to clip the other (home range), leaving only the area over water.

For the clipping code to work, both the HR polygon and the land polygon need to be in the same format (sf) and projection. After converting my HR polygon to sf format, the coordinates change (by a lot) and no longer overlap my land area. I've put together an example below.

library(sp)
library(sf)
library(ggmap)
library(adehabitatHR)
library(rnaturalearth)
library(lwgeom)

data<-data.frame(x=c(-50.3, -49.9, -50.0, -50.6, -55.3, -55.4, -55.5, -55.3, -54.9, -54.4, -51.5, -51.2, -50.8, -50.3),y=c(50.3, 48.8, 48.1, 47.4, 48.2, 47.4, 50.1, 48.1, 47.5, 50.7,50.4, 50.7, 50.5, 48.3))
data$id<-as.factor("a")

#create a SpatialPointsDataFrame by defining coordinates
coordinates(data)<-c("x","y")

#set CRS
proj4string(data)<-CRS("+proj=utm +zone=21 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")

# run code to determine minimum convex polygon
mcp<-mcp(data, percent = 95, unin = c("m"),
           unout = c("km2"))

#view resulting polygons by plotting them
plot(data, col=as.factor(data@data$id), pch=16)
plot(mcp, col=alpha(1:7, 0.5),add=TRUE)

#now plot polygons on map
#import map
canada <- ne_countries(country="canada",scale = "large", returnclass = "sf")

ggplot(data = canada) +
  geom_sf(data=canada, fill="gray80") +
  coord_sf(xlim=c(-58,-47), ylim=c(46, 52)) + # map coordinates
  labs(x = "Longitude", y = "Latitude", colour="") + #labels
  geom_polygon(data=mcp, fill="goldenrod",alpha=0.7, aes(x=long, y=lat, group=group))

# since these MCPs overlap land a lot, I want to clip away the parts that are on land. I will do this by clipping the MCP polygon with the map polygon.

#first, see what class of data the HR polygon and map polygon are
class(canada) # sf, so this is good.
class(mcp) # sp, so needs to be converted to sf

# coerce sp object to sf
mcp.poly <- st_as_sf(mcp,"+proj=utm +zone=21 +ellps=WGS84 +datum=WGS84 +units=m +no_defs",agr="constant",coords=c("x","y"))

class(mcp.poly) # now this is sf too

#next check that the CRS matches for both the HR polygon and land
#check utm format using:
st_crs(canada)
st_crs(mcp.poly)

#Since they differ, convert both to match
nl<-st_transform_proj(canada,"+proj=utm +zone=21 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")
mcp.poly<-st_transform_proj(mcp.poly,"+proj=utm +zone=21 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")

#plot these to check that they still make sense before clipping

ggplot() +
  geom_sf(data = nl, colour = "light gray", fill = "gray80")+
  geom_sf(data = mcp.poly,aes(group=id),fill = "gray")

#Coordinates are wrong
# map plots correctly alone:
ggplot() +
  geom_sf(data = nl, colour = "light gray", fill = "gray80")

# but the coordinates on the polygon have changed. Why?
ggplot() +
  geom_sf(data = mcp.poly,aes(group=id),fill = "gray")

## from here, I want to use the below code to clip the area of the polygon over land, and then determine the area of the polygon that covers water in km2

# difference between world polygons and the mcp
difference <- st_difference(mcp.poly, nl)

# coerce back to sp
difference <- as(difference, 'Spatial')

#determine area
poly.area<-area(difference)


texts: Convert part of a raster file to .asc file

texts:
I'm a bit new at this. I have a DEM map (.tif) of an area in QGIS. I know that I can transform it to a text file (.asc) with raster -> translate.

However, what I need is only a certain cross-section of the area. So I want to choose a straight line, with specific coordinates, and get an .asc file containing the data of only the points along this line.


texts:
I was wondering if PgQuery for QGIS is still available.  I am trying to install it on some new machines with QGIS 1.8.0 and I can not find it in the Python Plugin Installer.  Does anyone know of a place where I can download it directly?  Do I need to add a specific repository?  I currently only have the Official Repository loaded.

Thanks...


texts: What is causing invalid spatial reference error using ogr2ogr converting PostGIS to FileGDB in a Feature Dataset?

texts:
I am trying to use ogr2ogr to convert a layer from PostGIS to a FileGDB.  I want the layer to be created inside a Feature Dataset.

The spatial reference of the source layer is EPSG:3500.  The spatial reference of the Feature Dataset within the FileGDB is as follows, copied from the properties page:

NAD_1983_StatePlane_California_VI_FIPS_0406_Feet
WKID: 2230 Authority: EPSG

Projection: Lambert_Conformal_Conic
False_Easting: 6561666.666666666
False_Northing: 1640416.666666667
Central_Meridian: -116.25
Standard_Parallel_1: 32.78333333333333
Standard_Parallel_2: 33.88333333333333
Latitude_Of_Origin: 32.16666666666666
Linear Unit: Foot_US (0.3048006096012192)

Geographic Coordinate System: GCS_North_American_1983
Angular Unit: Degree (0.0174532925199433)
Prime Meridian: Greenwich (0.0)
Datum: D_North_American_1983
  Spheroid: GRS_1980
    Semimajor Axis: 6378137.0
    Semiminor Axis: 6356752.314140356
    Inverse Flattening: 298.257222101

My assumption is that the WKID: 2230 with Authority: EPSG, means that this should be spatial reference EPSG:2230

I think the tricky part here is that the spatial reference is already determined with the Feature Dataset.  I've seen errors rise before when trying to designate a spatial reference where the output is a Feature Dataset.

Given the above, I thought the best step would be to simply try to convert to EPSG:2230 without specifying the output spatial reference.  This resulted in an error.

Here is a sample command string and result:

ogr2ogr -overwrite -of "FileGDB" c:\fgdbData.gdb PG:"host=server port=5432 dbname=ib_data active_schema=ibcommunitydev user=gis_admin password=xxxx" cdpo_zoning -t_srs "EPSG:2230" -lco FEATURE_DATASET=fdCommunityDev

ERROR 1: Error: Failed at creating table for \fdCommunityDev\cdpo_zoning (Invalid spatial reference.)

I tried a number of different combinations , including specifying the source spatial reference as EPSG:2230 and letting it do a straight transformation, but nothing working.  They all returned the same invalid spatial reference error.

I tried using the same source layer and exporting to the File GDB outside of the Feature Dataset using the following command string, and it worked fine.

This string has to contain the variable -nlt MULTIPOLYGON, otherwise the table is imported as a non-spatial layer.

ogr2ogr -overwrite -f "FileGDB" c:\fgdbData.gdb PG:"host=server port=5432 dbname=ib_data active_schema=ibcommunitydev user=gis_admin password=xxxx" cdpo_zoning -t_srs "EPSG:2230" -nlt MULTIPOLYGON

The fact that the transformation to the spatial reference worked, leads me to believe it is an issue with the spatial reference of the Feature Dataset.


texts: Display XY Data points not displaying correctly?

texts:
I added a csv file with lat and long coordinates, and am trying to display the XY data. The X field is set to longitude, and Y field to latitude, the coordinate system and projection match the data frame, and I have taken out spaces in the csv file. However, when the points load, they project incorrectly and off and away from the rest of the map.

How can I fix this?

[


texts: Store reprojected coordinates in attribute table

texts:
I want to transform coordinates from gauss kruger to wgs84.
transforming ain't the problem as you can just set the layers reference to another but is there a way to print the new coordinates format into the attribute table?


texts: How to install PostGIS 1.5 on Mac OS 10.5?

texts:
I am running PostgreSQL 9.0 on Mac OS 10.5. I've installed postgis-pg90-1.5.3-1-osx through "Stack Builder" but still don't see postgis in my list of databases when I open PgAdmin. When I open my Qgis 1.7, I can connect to my postgres database and display all tables with no spatial geometry. However, when I try to query a non-spatial table, I get the following message "Your database has no working PostGIS support". I appreciate your ideas to get postgis to integrate into postgresql. Thanks!

Thanks for your ideas! Unfortunately, the first link (http://www.kyngchaos.com/software/postgres) is limited to "...Lion and Snow Leopard, 64bit only." I'm using Leopard (Mac OS 10.5). The second link (http://trac.osgeo.org/postgis/wiki/UsersWikiPackages) provides two options. Both of these (http://yum.postgresql.org/repopackages.php) & (http://www.postgresonline.com/journal/archives/204-postgis15-install-yum.html) point to Linux installations. So, I'm still stuck. I've downloaded PostGIS 1.5 using Application Stack Builder 3.0.0 from PostgreSQL 9.0. It installs, but does not run. I appreciate any other
suggestions.


texts: Network analyst Solve function getting stuck (never finishes)

texts:
I'm trying to create a series of Service Area polygons within a for loop. It performs OK creating the Service Area layer and loading Facilities, but for some features it begins running the arcpy.na.Solve function but never completes, just gets stuck running (no error, no crash). When done manually within Arc, the process is on the slow side but still works OK.

Strangely, on other features it works fine, there doesn't seem to be any difference  between the features which do and don't work, and it isn't related to size (sometimes has been left running for hours with only 8 features, but on others works fine when there's 100+).

Is there a way to improve it?

#Import system modules
import arcpy
from os import path
arcpy.env.overwriteOutput = True
arcpy.CheckOutExtension("Network")

# When run from within Arc...
streets = arcpy.GetParameter(0)
rls = arcpy.GetParameterAsText(1)
selection_field = arcpy.GetParameterAsText(2)
outGDB = arcpy.GetParameterAsText(3)

#Set environment settings
arcpy.env.workspace = outGDB

# Create list of selection values
#   This is required to process in batches (or else Arc falls over)
selection_values = []
with arcpy.da.SearchCursor(rls,"{}".format(selection_field)) as cursor:
  for row in cursor:
    selection_values.append(row[0])

# Loop through each selection value
for selection_value in selection_values:

    outNALayerName = "ServiceAreas_small"
    impedanceAttribute = "minutes"

    #Create a new SERVICE AREA layer.
    outNALayer = arcpy.na.MakeServiceAreaLayer(streets, outNALayerName,
                                  impedanceAttribute, "TRAVEL_FROM", "10",
                                  "DETAILED_POLYS", "NO_MERGE", "RINGS",
                                  polygon_trim = "TRIM_POLYS", poly_trim_value = 100)

    # Get the layer object from the result object. The service layer can now bereferenced using the layer object.
    outNALayer = outNALayer.getOutput(0)

    selection_code = selection_value.split(" ")[1]
    outFC = "ServiceArea_smaller_{}".format(selection_code)

    # Select by relevant selection value
    arcpy.SelectLayerByAttribute_management(rls, "NEW_SELECTION", """{0} = '{1}'""".format(selection_field,selection_value))

    # Get the names of the SUBLAYERS within the service area layer.
    subLayerNames = arcpy.na.GetNAClassNames(outNALayer)
    facilitiesLayerName = subLayerNames["Facilities"]
    polygonsLayerName = subLayerNames["SAPolygons"]

    # LOAD THE FACILITIES using default field mappings and default search tolerance
    arcpy.na.AddLocations(outNALayer, facilitiesLayerName, rls, "", "")

    # SOLVE the service area layer
    arcpy.na.Solve(outNALayer)

    # Get an object of the POLYGON layer
    PolygonsSubLayer = arcpy.mapping.ListLayers(outNALayer, polygonsLayerName)[0]

    # Create FEATURE CLASS
    arcpy.management.CopyFeatures(PolygonsSubLayer, outFC)

arcpy.SelectLayerByAttribute_management(rls, "CLEAR_SELECTION")


texts: How to create a lowest function for decimate_points() using lidR

texts:
How do I create a custom algorithm for decimate_points() in lidR analogous to the highest() algorithm but that selects the lowest points?
decimate_points(las, lowest(2))

I know that lidR allows to add custom algorithms as it was suggested by @JRR in this question: How -thin_with_grid works?


texts: Use Lat and Long to Zoom and Crop A Color GEOTIFF Based Map In R

texts:
I have a strong background in R but have only done limited GIS work.

I have a color Geotiff file that I downloaded from the USGS Earthexplorer. It is a "LandsatLook Image with Geographic Reference". https://earthexplorer.usgs.gov/

Given:

download.file("https://cdn.filesend.jp/private/4DxFm0MU9rS3nl-VdLqX1CKUGUBhasI0wQL4wg_solTUoS3QizjrJRQvub0IM8dl/LC08_L1TP_011031_20190731_20190819_01_T1.tif", destfile = "myraster.tif")
LL_Limits<-list(x=c(-70.4,-69.9), y=c(42.1,41.5)) # Lat/Long of Desired Plotting Location
POI<-c(-70.2,42.05) #Point Of Interest I Want To Mark

So far in my simple example I can load raster and look at it in color with,

library(raster)
x<-brick("myraster.tif")
plotRGB(x)

But I am unsure how to proceed with utilizing CROP and EXTENT to get the output I need. I want to crop to my desired lat and long, have the image be at the full resolution of the file, and put a marker on my point of interest. For reference the output of GDALinfo is:

library(rgdal)
GDALinfo("myraster.tif")

rows        7931
columns     7801
bands       3
lower left origin.x        319785
lower left origin.y        4503885
res.x       30
res.y       30
ysign       -1
oblique.x   0
oblique.y   0
driver      GTiff
projection  +proj=utm +zone=19 +datum=WGS84 +units=m +no_defs
file        myraster.tif
apparent band summary:
  GDType hasNoDataValue NoDataValue blockSize1 blockSize2
1   Byte          FALSE           0        256        256
2   Byte          FALSE           0        256        256
3   Byte          FALSE           0        256        256
apparent band statistics:
  Bmin Bmax Bmean Bsd
1    0  255    NA  NA
2    0  255    NA  NA
3    0  255    NA  NA
Metadata:
AREA_OR_POINT=Area

Any suggestions?


texts: Combining Google Maps coordinates with an existing shapefile

texts:
I would like to add GPS coordinates from Google Maps (see image below) as point data in QGIS and combine it with an existing shapefile [dvg1_EPSG25832_Shape.zip] https://www.opengeodata.nrw.de/produkte/geobasis/vkg/dvg/dvg1/.

Unfortunately the dots don't fit to the shapefile (see image below).

I think the problem is due to different Coordinate Reference Systems, but I don't know how solve it.


texts: From the city point I want to create multiple buffer zone for 8 directions. For example I am attaching an image. What is the procedure to do that in QGIS or any suitable method like Google Earth Engine.


texts:
I have calculated Kernel Density for feature class "fc" using following code:

dens = KernelDensity(fc,"NONE", 1, SearchRadius, "SQUARE_METERS", "DENSITIES", "PLANAR")

The raster attribute table:

I have then saved the raster as lyr file.

I can apply "Unique Values" in the symbology tab.
However:

lyr.symbologyType
Out[30]: u'OTHER'

How can i modify my raster, so that symbologyType = "UniqueValues"?

I am using:
ArcGIS Desktop 10.3.1
Python 2.7, 32bit

EDIT

I found a workaround:
Convert raster to polygon, setting field "Value" as raster field.
In the new polygon the field "GRIDCODE" now contains the unique values


texts:
I am trying to update parcel data for 29 different counties. Typically I would truncate the old data and append the new data, but clicking both those commands 29 times is time consuming.
I have written a script to automate the process. Truncate works, but I am getting stuck on the Append command. I am getting:

ERROR 000840: The value is not a Table View.
ERROR 000840: The value
is not a Raster Layer.

I think that the code is not recognizing the input value as individual feature classes and I am stuck on how to fix that. I am using a test environment with copies of 3 different counties.
import arcpy, os
from arcpy import env
env.overwriteOutput = True
env.workspace = "H:\Parcel Update Testing\Tax Parcels Test.gdb"

FCList = arcpy.ListFeatureClasses()
inpt = "H:\Parcel Update Testing\2020.gdb" # do i need to spell out individual Feature classes?
outLocation = r"H:\Parcel Update Testing\Tax parcels Test.gdb"
schemaType = "NO_TEST"
fieldMappings = ""
subtype = ""

for fc in FCList:
    arcpy.TruncateTable_management(fc)
    print("Truncating Table")

for fc in FCList:
    arcpy.Append_management(inpt, outLocation, schemaType, fieldMappings, subtype)
    print("Appending Table")


texts:
i'm using the split features tool to create two polygons from a bigger one. The problem is, when i save the split i have done, half of the polygon disappears. But sometimes, after doing some more polygons they appear again, and the last ones in with i was working disappear.
 Do you know what is happening and what is the solution?

PS: This Shape was made by other people and i don't know how it was made, neither in which program has been  created.


texts: Align images with self-defined CRS in Google Earth Engine

texts:
I have uploaded a dataset to GEE and I want to align every image in the collection to have the same geographic extent as the first image of this collection.

My code so far :

var ref_image = imageCollection.first();
var limited_collection = imageCollection.limit(10)

var align = function(image) {
  var im = image.select('b1');
  return image.reproject({
    crs: im.projection().crs(),
    scale: 1000})};

var new_collection = limited_collection.map(align)

print(new_collection)
print(imageCollection)
Map.addLayer(new_collection, {min:0, max:25, bands:['b1']}, 'new collection')

https://code.earthengine.google.com/9107bd80f72128963e617fe260d54ee4

The error I get is:

ImageCollection (Error)
Error in map(ID=0901012350):
Projection: Parameter 'crs' is required.

If I print out the CRS of the first image I get the information:

crs:
PROJCS["Radolan projection",
  GEOGCS["Radolan Coordinate System",
    DATUM["Radolan_Kugel",
      SPHEROID["Erdkugel", 6370040.0, 0.0]],
    PRIMEM["Greenwich", 0.0],
    UNIT["degree", 0.017453292519943295],
    AXIS["Longitude", EAST],
    AXIS["Latitude", NORTH]],
  PROJECTION["Polar_Stereographic"],
  PARAMETER["central_meridian", 10.0],
  PARAMETER["latitude_of_origin", 90.0],
  PARAMETER["scale_factor", 0.9330127019],
  PARAMETER["false_easting", 0.0],
  PARAMETER["false_northing", 0.0],
  UNIT["m*1000.0", 1000.0],
  AXIS["X", EAST],
  AXIS["Y", NORTH]]

So the CRS information is stored in the image. Is it only possible to align images with official CRS? Or have I misunderstood the .reproject command?


texts: ArcGIS Online producing isochrones which are skewed North-South, is there a way to fix this?

texts:
I am trying to run an isochrone analysis in ArcGIS Online using Analysis > Use Proximity > Create Drive-Time Areas.

The resulting polygons are stretched North-South, I suspect this is due to the analysis being done in WGS 84 instead of a projected coordinate system.

Is there a way to fix this? Possibly by running the analysis using a projected coordinate system?

EDIT:

As per Vince's request:

I used the Create Drive-Time Areas analysis tool with the following parameters:

Walking Time (5 minutes)
Towards Facility
Areas from different points - Overlap
Use current map extent - Unchecked

There is no option available to specify performing geodesic calculations.

The study area is at approximately 53° North.


texts: How to create raster stack from names on R

texts:
I have a dataframe containing the best predictors by their contribution to a model:

bestpredictors
           var contrib
2  Pfallwinter  8.1141
15      Tofall 58.0000
23          ai  5.3573
26     radmean  7.9753
27   radseason  0.9027

Every variable name is also a raster object. So, how can I make a raster stack with those variables? Like this:

predictors<-stack(Pfallwinter,Tofall,ai,radmean,radseason)

I have tried to extracting the string character and using the do.call and eval functions but seems that isn't the correct way, because they can't call the object already saved.


texts: How to connect GeoServer to a non-spatial database or Web service?

texts:
We have a system built using an Oracle database with geographical informations (latitude/longitude) but without spatial extensions (Locator). We would like to connect a GeoServer to it in order to deliver WFS services (i.e. physical values measured at some coordinates).
Is it possible to adapt a GeoServer for such a database access?

Our system delivers also REST web services able to give this same latitude/longitude/value information.

Can GeoServer translate it to WFS service?


texts: Compiling and Using GeoExplorer in Windows

texts:
I have been using GeoServer on Apache Tomcat in Windows and I want to install GeoExplorer. It is an OpenGeo Suite tool, and its source code is found here.

I have downloaded the source code and tried to compile it using WinAnt to generate a .WAR file. However, this resulted in error messages and I am stuck. Is my way correct or what must I do to fix this problem?

The error occurred when trying to build using WinAnt. The source code was added to C:\Temp\GeoExplore, then I navigated to that directory in a Command window and issued the ant command. The resulting errors are below, and are also linked at http://notepad.cc/errorGeo

C:\Temp\GeoExplorer>ant
Buildfile: C:\Temp\GeoExplorer\build.xml

init:
     [echo] pulling in submodules

BUILD FAILED
C:\Temp\GeoExplorer\build.xml:26: Execute failed: java.io.IOException: Cannot run program "git.cmd": CreateProcess error=2, System cannot find the file(translated to english)
        at java.lang.ProcessBuilder.start(ProcessBuilder.java:1029)
        at java.lang.Runtime.exec(Runtime.java:615)
        at org.apache.tools.ant.taskdefs.Execute$Java13CommandLauncher.exec(Execute.java:862)
        at org.apache.tools.ant.taskdefs.Execute.launch(Execute.java:481)
        at org.apache.tools.ant.taskdefs.Execute.execute(Execute.java:495)
        at org.apache.tools.ant.taskdefs.ExecTask.runExecute(ExecTask.java:631)
        at org.apache.tools.ant.taskdefs.ExecTask.runExec(ExecTask.java:672)
        at org.apache.tools.ant.taskdefs.ExecTask.execute(ExecTask.java:498)
        at org.apache.tools.ant.UnknownElement.execute(UnknownElement.java:291)
        at sun.reflect.NativeMethodAccessorImpl.invoke0(Native Method)
        at sun.reflect.NativeMethodAccessorImpl.invoke(NativeMethodAccessorImpl.java:57)
        at sun.reflect.DelegatingMethodAccessorImpl.invoke(DelegatingMethodAccessorImpl.java:43)
        at java.lang.reflect.Method.invoke(Method.java:601)
        at org.apache.tools.ant.dispatch.DispatchUtils.execute(DispatchUtils.java:106)
        at org.apache.tools.ant.Task.perform(Task.java:348)
        at org.apache.tools.ant.Target.execute(Target.java:390)
        at org.apache.tools.ant.Target.performTasks(Target.java:411)
        at org.apache.tools.ant.Project.executeSortedTargets(Project.java:1399)
        at org.apache.tools.ant.Project.executeTarget(Project.java:1368)
        at org.apache.tools.ant.helper.DefaultExecutor.executeTargets(DefaultExecutor.java:41)
        at org.apache.tools.ant.Project.executeTargets(Project.java:1251)
        at org.apache.tools.ant.Main.runBuild(Main.java:809)
        at org.apache.tools.ant.Main.startAnt(Main.java:217)
        at org.apache.tools.ant.launch.Launcher.run(Launcher.java:280)
        at org.apache.tools.ant.launch.Launcher.main(Launcher.java:109)
Caused by: java.io.IOException: CreateProcess error=2, Sistem belirtilen dosyay?
 bulam?yor
        at java.lang.ProcessImpl.create(Native Method)
        at java.lang.ProcessImpl.(ProcessImpl.java:189)
        at java.lang.ProcessImpl.start(ProcessImpl.java:133)
        at java.lang.ProcessBuilder.start(ProcessBuilder.java:1021)
        ... 24 more

Total time: 0 seconds


texts: I want latitude longitude information as an ASCII file for a location of my choice in 1 km resolution. The picture shows data for each grid point. Now I want the lat long info for each grid, also as an Ascii file. How do I get it?


texts:
I am running

r.neighbors --overwrite --verbose input=bf03pits@PERMANENT output=bf03pitsdenoise method=mode size=5

in the GRASS GUI using the r.neighbors dialog.

bf03pits has 3955848 all with values either 0 or 1.

When I run r.neighbors the result has 1 cell with a range of -nan to -nan.

Yet the output shows no errors.

Why might r.neighbors be failing? How can I go about doing what I want?


texts: Making a buffer cover a specific area (need to calculate overlap)

texts:
I work in precision agriculture and use QGIS. Today's technology makes it possible to adjust the spreading width of, for example, sprayers and manure spreaders according to where you have driven before, to reduce overlap. Such implements can have a spreading width of over 30 m. The overlap, when driving towards the headland at an angle <> 90 degrees, will be significant when considering investing in such equipment.
So, to the point. I want to make a tool in QGIS Graphic Modeler with the field boundary and the width of the implement and direction of travel as input. QGIS would then calculate the overlap. I have reached a point where I have made a headland (buffer on the inside of the border with the width corresponding to the tool). After specifying a direction of travel, I have made parallel lines indicating the center point of the track. These are then "clipped" to the headland.

A buffer on each side of the track (with a flat end cap) marks the area worked. The result is then an overlap on one side of the track and a "negative" overlap on the other side of the track (see figure) when I enter the headland on an angel <> 90 degrees.
As I want to automate the process, I need help to figure out how to extend the worked area a little longer on each end to just cover the headland as the figure shows.
Will it be possible to calculate the maximum extent of the worked area by figuring out point A and B, then extend the rectangle to that pint and then repeat this process for all the worked rectangles?
https://i.stack.imgur.com/IwcaY.png


texts: Network Analysis Layer in ArcGIS Mobile Map Cache

texts: Does anyone know how to extract a Network Analysis Layer into an ArcGIS Mobile Map Cache?


texts: Using newly created features in next part of script using ArcPy

texts:
I am trying to write a python script which buffers a large number of features from a geodatabase (but not all) and then dissolves the newly created buffer shapefiles together. Currently I have only been able to achieve this by listing the specific locations of the newly created buffer shapefiles (which is OK for a couple of features but a nightmare with more). I am new to Python and unsure how to references too (or ideally make a list) of the newly produced buffer shapefiles without individually listing there locations. This is also further complicated by producing 3 different size buffers for each feature class. Here is the first bit of my script before I carry out the dissolve.
# imports arcpy and os
import arcpy, os
from arcpy import env

# sets the workspace to my gdb
arcpy.env.workspace = "C:/GIS Home/Mining Features (MATTHEW).gdb"

# Output workspace location is defined
output_workspace = "C:/GIS Home/project_1"

# list of features to buffer in gdb
fcList = ["Shaft_Metalliferous", "Quarry_Metalliferous"]

for featureClass in fcList:

arcpy.Buffer_analysis(featureClass, output_workspace + "/" + featureClass[:-4] + "_b250.shp", 250, "", "", "ALL")
arcpy.Buffer_analysis(featureClass, output_workspace + "/" + featureClass[:-4] + "_b500.shp", 500, "", "", "ALL")
arcpy.Buffer_analysis(featureClass, output_workspace + "/" + featureClass[:-4] + "_b750.shp", 750, "", "", "ALL")


texts: R Raster: extract weighted mean within circle of specific radius

texts:
I have a raster object from which I want to extract area weighted mean values from a circle with a specific radius in m. The raster has latitude and longitude wrt to the WGS84 datum.

I can extract values within a circle using extract and specifying a buffer, but how can I get a weighted mean?

The raster layer can be downloaded here.

WOA <- readRDS('WOA.RDS')
crs(WOA)
CRS arguments:
  +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0

xy <- data.frame(x = -40, y = 60)

# extract values from cells from circle with 100 km radius
extract(WOA, xy, buffer = 1e+5)

but to make use of the weights in extract I think need to convert xy to a SpatialPoint

xy_sp <- SpatialPoints(xy, proj4string = CRS('+proj=longlat +datum=WGS84'))

and use gBuffer from the rgeos package to expand the point to a circle with the desired radius (i.e. a SpatialPolygon), but how do I take care of the units since xt_sp is in degrees? Can I use spTransform (rgdal) and what (metric) coordinate system would be appropriate for my global dataset?


texts: How to overlay my interpolation within the boundaries of Quebec?

texts:
I am trying to create interpolated maps within geographical borders.

I have the coordinates and values to the points I am trying to interpolate. I also have a map of the geographical land mass in .shp file, within which I am trying to confine the interpolation.

Using GRASS in QGIS, I have managed to interpolate my points with v.surf.idw. But I seem unable to confine them to the map of my geographical region (Northern Quebec, bordering James and Hudson Bay), meaning that my gradient goes into the bay!

How do overlay my interpolation within the boundaries of Quebec?


texts: Seeking tools to convert ways of osm data to format including coordinates?

texts:
I have asked Getting coordinates of ways in *.pbf file? for getting coordinates of way data in osm. In detail, the ways including some kinds of amenities and tourist attractions of whole Japan. So the volume is very large (That is one of the reasons why I do not choose to use Overpass API). If using the code

osmosis.bat --read-pbf JP.pbf --tf accept-ways tourism=* --tf reject-relations --used-node  --write-xml de-pipeline.osm

The output did include the coordinates but it is like this:

<node id="1" lat="" lon=""/>
<node id="2" lat="" lon=""/>
<way id="1">
 <nd ref="1"/>
 <nd ref="2"/>
</way>

And if I want to get the coordinates of a way I need to
traverse the entire document which I think will take a lot of time if the file is large. So I am wondering whether there are some tools to directly get each way as an element(object) that includes the node coordinates like this (or any other formats that do not need to traverse the file again):

<way id="1">
 <node id="1" lat="" lon=""/>
 <node id="2" lat="" lon=""/>
</way>


texts: Focal statistics means for very large area of NDVI raster to be used to extract value to points

texts:
Environment: ArcPro running on server, NDVI from Living Atlas, 30,000 points covering the entire eastern half of Colorado.
Task: get the mean NDVI mean value around each point for 270 meters
So far my software has crashed twice running focal statistics I need for the value extraction. I split my area in half to see if that helps and it's now been processing for 40 hours (at least it hasn't crashed yet). Is there another more efficient method to derive the mean values around each point? I still need to do the other half and really don't have 4 more days to get this one task done. I am also open to R or Python code-based solutions.


texts:
I want to show an image in a menu in an ArcMap AddIn toolbar. Here is my code in the  config.esriaddinx:

<Menus>
       <Menu id="GridViewMenu" caption="Grid View" image="Images\GridView.png">
          <Items>
            <Button refID="Button1" />
            <Button refID="Button2" />
            <Button refID="Button3" />
            <Button refID="Button4" separator="true"/>
            <Button refID="Button5" />
            <Button refID="Button6" />
          </Items>
        </Menu>
</Menus>

Why does it not show an image in menu?


texts: How to open RS2 AWiFS imagery on SNAP software?

texts:
I am working on an Indian satellite, ResourseSat2 AWiFS imagery. The input imagery is in erdas imagine '.img' format.  I would like to use the tools available on ESA SNAP software.
File--> Open--> image_directory is not able to open the imagery
I am getting an error "java.lang.NullPointerException. Failed opening a dataset...
I am able to successfully open the imagery on QGIS, but not on SNAP software.
Is there any way to make SNAP software to read the RS2 AWiFS imagery?


texts: PostGIS spatial reference ID (SRID) for regular cartesian coordinate system?

texts:
I would like to store a polygon as a set of 2D points represented as x, y coordinates of a Cartesian coordinate system.

Which spatial reference can I use? Looking at my spatial_ref_sys table, all the spatial references seem to be geography related.

The points of the polygon would represent satellite measurements which I would invoke ST_ConvexHull on to get the shape of the satellite footprint.