texts: Listing datasets doesn't bring back the correct result

texts:
If I do a ListDatasets() in python, then I don't get the correct result in the list except for datasets belong to other schema owners, so why can I still see the datasets in ArcCatalog.

The code is staightforward:

workspace = r'Database Connections\\enidb.sde'
arcpy.env.workspace = workspace
dslist = arcpy.ListDatasets()

The dslist is empty bar the CULTURAL.ESRI_DATASETS - I should be able to see all the ENIDB ones as I am logged in as that schema owner:

Any clues folks?


texts: Changing Background Color of Selection Rows using MapBasic?

texts:
In MapInfo Pro 12.5, in the Browser window, selected rows (selected either by clicking rows with the mouse or via a query) are highlighted dark blue but running the same queries in a MapBasic program selection rows in the browser window are highlighted light grey.

In MapBasic is there a system setting to change the selected row highlight colour in the browser window (either to a deeper grey or dark blue). I have checked the MapBasic Reference manual, MapInfo System settings, googled etc. but have not been able to find anything.

Where might I be able find this in the Reference or User manuals?


texts: Using "Merge Divided Roads" tool with lesser parallel lines

texts:
When using the "Merge Divided Roads" tool, it merges roads that "trend generally parallel to one another". This means that from figure one.. It will merge the beginning, where the roads are close to each other (so the result is like figure 2), but not the entire length like I want it to (figure 3). I believe this is because the lines are not parallel enough (the angle between some of the lines are too big). Because even if I extend the merging distance to 300m (which is more than how far the roads are apart), I still don't get the result from figure 3.

Anybody know how to decrease how parallel the lines have to be to get merged?


texts: QGIS and SQLite wrongly changing layer type on import

texts:
I have been struggling for a few days now with importing Linestring
layers into Spatialite.

Every time I save the layers into a sqlite file, they are automatically
converted to MultiLinestring, even though they are perfectly normal
Linestring files in all other formats I tested (SHP and TAB).

Has anybody come across this before? Is this a bug?  Or am I doing it wrong?


texts: Does ArcPy leak memory when importing CSV files?

texts:
I'm trying to import a great number of csv files of about 250MB each. I wrote a python script for that as follows:

import gc
import glob
tables = glob.glob("C:\\Users\\emportella\\Documents\\ArcGIS\\2517")
for table in tables:
    arcpy.TableToGeodatabase_conversion(table, "C:\\Users\\emportella\\Documents\\ArcGIS\\2517\\Tables.gdb")
    gc.collect()

It works for the first 4 to 6 files and without giving any warning is closes. i've payed attention to memory usage in Windows Task Manager and the process cosumes more and more memory until it craches.

I've already tried using some different aproches like using:

env.workspace = "C:\\Users\\emportella\\Documents\\ArcGIS\\2517\\*.csv"
tables = arcpy.ListTables()
arcpy.TableToGeodatabase_conversion(tables, "C:\\Users\\emportella\\Documents\\ArcGIS\\2517\\Tables.gdb")

The TableToGeodatabase_conversion accepts list of tables but it craches the same way.

So there is the question is there a way to manage the memory used in arcpy framework?


texts: Difference between two Z-values extracted using "Extract -Z values" tool in QGIS 3

texts:
When I run Extract-Z tool (vector geometry) and select Q3, I get two values for each object:

In many cases, the two values are almost identical, but in some, they differ quite a bit. There is no explanation in the docs for these two values.

Can anybody explain what the difference is?


texts: Update Multiple Records with Unique Values

texts:
I'm looking for a way to programmatically update multiple records in attribute tables.  The information to be updated for each record is unique.

History of Problem: Datafile names were truncated at the time of shapefile creation.  The attribute tables for the features have been exported to a spreadsheet and the full datafile names have been identified.  Basically, I have a spreadsheet with OBJECT IDs and the corrected full datafile name that needs to overwrite what is currently in the attribute table for the record.

It makes sense to me that it would be some sort of loop, but I'm not a savvy programmer and the fact that the replacement text is unique for most of the records makes it harder for me to wrap my head around.

I'm using ArcMap.


texts:
Using ggmap, I want to include municipality borders (polygon) from a shapefile on a map containing some location points. This script does everything except plotting the polygon:

library(rgdal)
library(ggmap)

# Get shapefile with Drammen municipality borders
tmpzip<-tempfile()
tmpdir<-tempfile()
dir.create(tmpdir)
download.file("http://www.kartverket.no/Documents/Kart/N50-N5000%20Kartdata/33_N5000_shape.zip",tmpzip)
unzip(tmpzip, exdir=tmpdir)
kommune <- readOGR(dsn=tmpdir, layer="NO_AdminOmrader_pol")
kommune<-kommune[kommune$NAVN=="Drammen",]
kommune<-spTransform(kommune, CRS("+init=epsg:4326"))

# Get location point data
subscr<-data.frame(lon=c(10.1237,10.2161,10.2993),lat=c(59.7567,59.7527,59.6863), pop=c(58,12,150))
coordinates(subscr)<-~lon+lat
proj4string(subscr)<-CRS("+init=epsg:4326")

lon <- c(10.0937,10.3293)
lat <- c(59.7916,59.6563)
map <- get_map(location = c(lon[1], lat[2], lon[2], lat[1]),
               maptype = "roadmap", source = "osm", zoom = 11)
p <- ggmap(map) +
  geom_point(data = as.data.frame(subscr), aes(x = lon, y = lat, size=pop),
             colour = "darkgreen") +
  theme_bw()
print(p)

How can I plot the polygon from the shapefile? I have tried substituting the second last line with the following:

p <- ggmap(map) +
  geom_point(data = as.data.frame(subscr), aes(x = lon, y = lat, size=pop),
             colour = "darkgreen") +
  geom_polygon(data = as.data.frame(kommune)) +
  theme_bw()

But then I get the following error:

Error: Aesthetics must be either length 1 or the same as the data (1): x, y


texts:
I want to copy a raster file into SQLite database. Previously, I did the same with vector data like this:

def store_vector_output(self, file_name, identifier):
    from osgeo import ogr
    drv = ogr.GetDriverByName("SQLite")
    dsc_out = drv.CreateDataSource(self.dblocation)
    layer = dsc_out.CopyLayer(dsc_in.GetLayer(), identifier,
                              ['OVERWRITE=YES'])

    dsc_out.Destroy()
    dsc_in.Destroy()

With raster data, I expected it to be similar:

def store_raster_output(self, file_name, identifier):
    from osgeo import gdal

    dsc_in = gdal.Open(file_name)
    drv = gdal.GetDriverByName("SQLite")
    dsc_out = drv.Create(self.dblocation, xsize=30, ysize=30)
    layer = dsc_out.CopyLayer(dsc_in.GetLayer(), identifier,
                              ['OVERWRITE=YES'])

    dsc_out.Destroy()
    dsc_in.Destroy()

However, in store_raster_output, dsc_in.GetLayer() is None. There are no layers in dsc_in (its datatype is <class 'osgeo.gdal.Dataset'>)

So, I tried to use dsc_in.CreateLayer("some_name") but that resulted in this:

ERROR 6: CreateLayer() not supported by this dataset.

Therer is also a method ReadRaster() but since I am using CopyLayer the argument must be of datatype layer.

How do I write a raster file to SQLite?

Perhaps using a different command than CopyLayer?

I didn't find any such command here but maybe I just didn't know where to look.


texts: How to randomly reclassify raster pixel values in ArcGIS?

texts:
I have a raster dataset (USDA Cropscape dataset) where each pixel is given a crop or landcover type (eg, corn, rice, urban area, forest). Some of the pixel values are classified as more than 1 crop/landcover type. For example, there are 40,000 pixels that are classified as “Double Crop Barley/Corn” as the class name. This is because certain acres of land engage in crop rotation practices. I have attached a screenshot of the attribute table.

I would like to reclassify the pixels (shown in 'Count' column of attached screenshot) as either Barley OR Corn (but not both), so that 20,000 pixels are classified as barley and 20,000 are corn. Now for the important part: I would like to randomly distribute the reclassification so that there is not a clustering of one crop type (ie, corn or barley) but an even and random distribution of each reclass across the entire extent. I don't want the first 20,000 pixels to be corn and the next to be barley, but every other pixel being classified as either one would be nice.

How would you go about this? General outline or step by step guidance, if possible. The important part is to get half of the pixels reclassified in randomized distribution


texts: Adding photos as attachments to geotagged photo points in ArcMap?

texts:
I'm using ArcGis 10.2, I have photos in .jpeg that are geotagged with lat/long. I'm using Data Management Tools > Photos > GeoTagged Photos to Points tool in the ArcToolbox. So far I have been able to get it to successfully import the photo locations as points onto my map. However I would like to include the photos as attachments so that when I click on the individual point I can see the photo it represents.

On the tool pop up, there is the option to check "Add Photos As Attachments" but I can't select it. It's grayed out and when I click on it the message on the side of the popup reads: "the output feature class must be in a version 10 or higher geodatabase for photo files to be added as attachments".

I'm not sure what it means by a version 10 or higher geodatabase. The output feature class sections fills out automatically to the ArcGIS folder on my E drive. I've tried changing it so that it's the same folder as the input folder but I still can't select the "Add Photos" option. Has anyone else had this issue and/or does someone know what it means by my choice of geodatabase?

Update: I was able to find a solution that produced different but similar results. I enabled hyperlinks in the properties of the point shapefile for the photos. If I use the hyperlink tool (little lightning bolt) and click on an individual point the actual picture opens in a photo viewer on my computer. So the links to the network are working and I can see which photo corresponds with which point (ultimately what I needed) but I would still like to see if I can find a way to add photos as attachments as mentioned above.


texts: Add layer to GeoServer workspace which is not default

texts:
Adding a layer to GeoServer in the default workspace, I do it this way:

Layers> "Add a new layer"
In the "Add layers from" menu, select the source
Click the "Configure new SQL view" link
Type view name and SQL statement
If the query is correct, the layer attributes will be displayed
Check the Identifier checkbox next to the Id line
In the geometry line, select the geometry type (eg polygon)
Enter 3765 in the SRID column
Save
The editing form opens
In the "Bounding Boxes" section, click the "Compute from SRS bounds" and "Compute from native bounds" links
12.save

How can I add a layer to workspace that is not default?
There is some checkbox or I need to change my default workspace?


texts: QGIS Python error on OS X while trying to simplify shapefile geometries

texts:
I am trying to reduce the file size of the shapefile by simplifying geometry in QGiS on OS X.

But I get the following error...

An error has occured while executing Python code:

Traceback (most recent call last):
  File "/Applications/QGIS.app/Contents/Resources/python/plugins/fTools/tools/doSimplify.py", line 255, in run
    self.runSimplify()
  File "/Applications/QGIS.app/Contents/Resources/python/plugins/fTools/tools/doSimplify.py", line 280, in runSimplify
    shapeFileWriter = QgsVectorFileWriter( self.outputFileName, self.outputEncoding, shapeFields, wkbType, crs )
TypeError: arguments did not match any overloaded call:
  QgsVectorFileWriter(QString, QString, dict-of-int-QgsField, QGis.WkbType, QgsCoordinateReferenceSystem, QString driverName="ESRI Shapefile", QStringList datasourceOptions=QStringList(), QStringList layerOptions=QStringList(), QString newFilename=None): argument 2 has unexpected type 'NoneType'
  QgsVectorFileWriter(QgsVectorFileWriter): argument 1 has unexpected type 'QString'

Python version:
2.7.1 (r271:86832, Jun 16 2011, 16:59:05)
[GCC 4.2.1 (Based on Apple Inc. build 5658) (LLVM build 2335.15.00)]

QGIS version:
1.8.0-Lisboa Lisboa,

Python path: ['/Applications/QGIS.app/Contents/MacOS/../Resources/python', '/Users/admin/.qgis//python', '/Users/admin/.qgis//python/plugins', '/Applications/QGIS.app/Contents/MacOS/../Resources/python/plugins', '/Library/Frameworks/SQLite3.framework/Versions/3/Python/2.7', '/Library/Python/2.7/site-packages/numpy-override', '/Library/Frameworks/GDAL.framework/Versions/1.9/Python/2.7/site-packages', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python27.zip', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/plat-darwin', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/plat-mac', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/plat-mac/lib-scriptpackages', '/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/lib-tk', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/lib-old', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/lib-dynload', '/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python/PyObjC', '/Library/Python/2.7/site-packages', '/Applications/QGIS.app/Contents/Resources/python/plugins/fTools/tools']

Any idea what the issue could be ?

Gaurav


texts:
I am using QGIS 2.4.0 and I have a vector layer which contains a remoteness indicator which is broken up into 5 categories. I would like to turn this into a layer that has these 5 categories collapsed into 2. On the included map I would like to join the green areas together into one group and the yellow and red areas into another group.

In the attribute table included I would like to join together the categories in the RA_Name11 field.


texts: How to optimize API calls for a large dataset using Python?

texts:
Objective: Send a list of addresses to an API and extract certain information(eg: a flag which indicates if an address is in a flood zone or not).

Solution: Working Python script for small data.

Problem: I want to optimize my current solution for large input. How to improve the performance of the API calls. If I have 100,000 addresses will my current solution fail? Will this slow down the HTTP calls? Will I get a request TIME out? Does the API resist the number of API calls being made?

Input: a list of addresses

Sample input

777 Brockton Avenue, Abington MA 2351

30 Memorial Drive, Avon MA 2322

Output: A data frame of all addresses and a flag which indicates if it is in a flood zone or not.
API: https://hazards.fema.gov/gis/nfhl/rest/services/public/NFHL/MapServer/28/query.

My current solution works well for a small dataset.

# Creating a function to get lat & long of the existing adress and then detecting the zone in fema
def zonedetect(addrs):
    global geolocate
    geocode_result = geocode(address=addrs, as_featureset=True)
    latitude = geocode_result.features[0].geometry.x
    longitude = geocode_result.features[0].geometry.y
    url = "https://hazards.fema.gov/gis/nfhl/rest/services/public/NFHL/MapServer/28/query?where=1%3D1&text=&objectIds=&time=&geometry="+str(latitude)+"%2C"+str(longitude)+"&geometryType=esriGeometryPoint&inSR=4326&spatialRel=esriSpatialRelIntersects&relationParam=&outFields=*&returnGeometry=true&returnTrueCurves=false&maxAllowableOffset=&geometryPrecision=&outSR=&returnIdsOnly=false&returnCountOnly=false&orderByFields=&groupByFieldsForStatistics=&outStatistics=&returnZ=false&returnM=false&gdbVersion=&returnDistinctValues=false&resultOffset=&resultRecordCount=&queryByDistance=&returnExtentsOnly=false&datumTransformation=&parameterValues=&rangeValues=&f=json"
    response = req.get(url)
    parsed_data = json.loads(response.text)
    formatted_data = json_normalize(parsed_data["features"])
    formatted_data["Address_1"] = addrs

    #Exception handling
    if response.status_code == 200:
        geolocate = geolocate.append(formatted_data, ignore_index = True)
    else:
        print("Request to {} failed".format(postcode))

# Reading every adress from existing dataframe
for i in range(len(df.index)):
    zonedetect(df["Address"][i])

Instead of using the for loop above is there an alternative. Can I process this logic in a batch?


texts: Which map API to draw polygon and export my data that is within boundary?

texts:
I have 5 point datasets that I want to display and be able to filter by fields. I will then be able to hand draw a polygon and export Excel of the data from each dataset within the polygon.

Which web map API should I use? I was looking at mapbox/turf and geojson, but this seems limited.


texts: Using label properties/expression in ArcMap to label several ranges of numbers within feature class red?

texts: I have about a thousand labels, some of which need to be labeled red. For example, 511-533 and 659-671 need to be labeled red, and the rest can be labeled in black. They are all part of the same field and feature class. I've already used python to filter out the even numbers (different part of the work flow) and need to move on to coloring said number ranges. Any ideas?


texts: Polygon creation (Clockwise rotation or not)

texts:
According to "ESRI Shapefile Technical Description An ESRI White Paper—July 1998", Vertices for a single, ringed polygon are, therefore,always in clockwise order. Otherwise,"Dirty" polygons occur when the rings that define holes in the polygon also go clockwise, which causes overlapping interiors.

In my understanding, you can create donuts polygon with "avoid intersection of new polygons" is checked on snapping option, even if both polygon(inner and outer) are created clockwise roated.

It has not necessary to care about clockwise rotation when you create polygon in QGIS.

Is this a right understanding?


texts: Joining points to points on the same road

texts:
I know there may be a really simple solution to this query but I have been looking for one for ages and have not been successful.

I have points that represent mammal sightings on roads, as well as points that represent traffic volume in the UK. I want to join the mammal points to the volume points so I know the traffic volume at each mammal point.

I have tried joining point to point, where the mammal points join to the closest volume point, which worked for some cases, but in most cases the mammal points joined to a volume point that was on a different road. I have also tried converting the volume points to a line but that doesn't seem to have worked either.

Does anyone know a way I can join 1 point to the 5 nearest points, or to join the mammal points to volume points that are on the same road as it?


texts: How to change black edge into transparent edge using GDAL?

texts:
I am trying to make seaice map in the Arctic Ocean.

I used a code like this.

gdal_translate -outsize 100% 100% -a_nodata value -co format=png -of KMLSUPEROVERLAY input.tif output.kmz

However when I upload kmz file into Google Earth, there is black edge.

How can I transform black edge into transparent edge?


texts:
I am trying to write a function to reproject a vector file , by taking the original shapefile as input and giving out the name of the reprojected layer . Though when I call this function I get an error saying The file size exceed 2GB, and I don't see why.

The file has only one polygon feature, though i wish to write a generalized code, i am specifying the type as

outlayer=outputdata.CreateLayer('',geom_type=ogr.wkbPolygon)

for simplicity so that if it works for this , would probably iteratively develop extra functionality

I am fairly new to programming using ogr.

This is what it looks like

#transformataion from one system to the other using ogr
from osgeo import ogr,osr
from epsg_ident import EpsgIdent
from os.path import basename,splitext,exists
def transformfromshapefile(file):
    driver=ogr.GetDriverByName("ESRI Shapefile")
    #to identify the epsg code of the input layer
    filebasename=splitext(basename(file))[0]
    fileprj=filebasename+'.prj'
    ident=EpsgIdent()
    ident.read_prj_from_file(fileprj)
    in_epsg_value=int(ident.get_epsg())
    src_epsg=osr.SpatialReference()
    target_epsg=osr.SpatialReference()
    src_epsg.ImportFromEPSG(in_epsg_value)
    target_epsg.ImportFromEPSG(4326)
    #define the transformation
    transformation=osr.CoordinateTransformation(src_epsg,target_epsg)
    #create an output layer if does not exist, or delete and create it if it already does
    outputfile="projected.shp"
    if(exists(outputfile)):
        driver.DeleteDataSource(outputfile)
    outputdata=driver.CreateDataSource(outputfile)
    #set the fields for the output layer using the fields of the input layer, which is inlayer, outputlayer is outlayer
    inputshapefile=driver.Open(file)
    inlayer=inputshapefile.GetLayer()
    outlayer=outputdata.CreateLayer('',geom_type=ogr.wkbPolygon)
    inlayerdefinition=inlayer.GetLayerDefn()
    for i in range(inlayerdefinition.GetFieldCount()):
        fieldDefn=inlayerdefinition.GetFieldDefn(i)
        outlayer.CreateField(fieldDefn)
    #Get the feature definition of the outputlayer
    outlayerdefn=outlayer.GetLayerDefn()
    #Transform the outputlayer features , looping through each feature
    infeature=inlayer.GetNextFeature()
    while infeature:
        geom=infeature.GetGeometryRef()
        #Reproject the geometry
        geom.Transform(transformation)
        #create a new feature that is transformed
        newfeature=ogr.Feature(outlayerdefn)
        #Set the geometry for the givenfeature
        newfeature.SetGeometry(geom)
        #Set the fields
        for i in range(outlayerdefn.GetFieldCount()):
            newfeature.SetField(outlayerdefn.GetFieldDefn(i).GetNameRef(),infeature.GetField(i))
        #Add the feature to the shapefile
        outlayer.CreateFeature(newfeature)
        #Dereference the outfeature for the next loop
        newfeature=None
        inlayer.GetNextFeature()

The size of the input .shp file is 300KB and for the .dbf file it is 1KB


texts:
I'm trying to join touching lines within a linestring dataset, so that I end up with a single feature for each set of touching features (eg the set outlined in green below). I've tried pretty much every algorithm I can find, including dissolve, merge, union (in QGIS, including the GDAL and SAGA options), ST_LineMerge in PostGIS, but the lines refuse to join, even though the end nodes of the segments are coincident. Can anyone suggest a process which would do this? I'm running this as part of a processing model in QGIS which has several hundred thousand records, so it needs to be a process I can include in this.
The next part of the process uses 'split with lines' to split the feature into components at line joins, so in the example below I need to end up with three features.

I can manually merge the features into one using the 'merge selected features' tool, but I can't see how to run this as a process, and even so, running 'split with lines' on the result splits the lines as they were initially, resulting in four features not three.
My guess is that there is something about the data which is causing this, but I can't see a way of fixing it, or running the merge.


texts: Spatial Data from Postgis to Java/REST with RestEasy + Jackson

texts:
I'm developing an application where I take a route from google maps, save it in postgresql/postgis. I'd like to take this route (LineString geometry) and send it back to Google Maps.

I already can save the map in the database, now I want to show the routes on a map.

On the client side I'm using JQuery (AJAX) + Google Maps Api.
On the server side I'm using REST services with RestEasy.

I tried to work with Jackson to transform the Entity to JSON, but, the Jackson fails to transform geometric data in json. Does anyone know how to transform this spatial data into JSON?

My Entity:

public class Rota {
...
private LineString caminho;
public Rota() {
}
public Rota(Integer id, Usuario usuario, String descricao, LineString     lineString) {
    this.id = id;
    this.usuario = usuario;
    this.descricao = descricao;
    this.caminho = lineString;
}
// getters and setters
@Type(type = "org.hibernate.spatial.GeometryType")
public LineString getCaminho() {
    return caminho;
}
public void setCaminho(LineString caminho) {
    this.caminho = caminho;
}
// hashcode and equals method
}

My REST service:

@Inject
private RotaBC rotaBC;

@GET
@Path("/rotas/{idrota}")
@Produces("application/json")
public Rota show(@PathParam("id") String id) {
    return rotaBC.show(Integer.parseInt(id));
}

The error:

  15:05:16,447 ERROR
  [org.apache.catalina.core.ContainerBase.[jboss.web].[default-host].[/sharecar].[default]]
  (http--0.0.0.0-8080-1) Servlet.service() for servlet default threw
  exception: java.lang.ClassNotFoundException:
  org.codehaus.jackson.map.JsonMappingException$Reference from [Module
  "org.codehaus.jackson.jackson-mapper-asl:main" from local module
  loader @1d332b (roots: /opt/demoiselle/server/jboss-7.1/modules)] ...
  at
  org.jboss.modules.ConcurrentClassLoader.loadClass(ConcurrentClassLoader.java:120)
  at
  org.codehaus.jackson.map.ser.std.BeanSerializerBase.serializeFields(BeanSerializerBase.java:166)
  ... at
  org.codehaus.jackson.map.ser.BeanPropertyWriter.serializeAsField(BeanPropertyWriter.java:446)
  at
  org.codehaus.jackson.map.ser.std.BeanSerializerBase.serializeFields(BeanSerializerBase.java:150)
  at
  org.codehaus.jackson.map.ser.BeanSerializer.serialize(BeanSerializer.java:112)


texts: Displaying addresses (point) relative to quantity at address using ArcGIS Desktop?

texts:
I just geocoded 40,000 plus addresses.

However there are several addresses that have 100's of people living at a given address (ex. apartment building).

How do I show larger values as larger symbols using ArcGIS Desktop 10.2?


texts:
I was wondering why sf::st_intersection() takes too much time, in comparison with sp::over(); I use sf for pretty much all of my geospatial workflows (I was used to sp as well), and more and more packages are building support for sf objects; is there a way for speeding up st_intersection? According to st_intersection man page: "A spatial index is built on argument x", so I can't build a spatial index to speed it up.

In the following example eke is a point (sp or sf) object and ver_poly is a polygon object, even coerced to "Spatial", sp::over() goes faster (just using 10000 points, whole object is 150000 features):

> system.time(en_ver <- over(eke[1:10000,], as(ver_poly, "Spatial")))
   user  system elapsed
  0.573   0.004   0.576
> system.time(eke[1:10000,] %>% st_as_sf(coords = c("longitud", "latitud")) %>% st_set_crs(4326) %>%
+   st_intersection(ver_poly))
although coordinates are longitude/latitude, st_intersection assumes that they are planar
   user  system elapsed
 83.552   0.004  83.585


texts:
I've created a python plugin for QGIS. Before I upload it to the QGIS plugins repository, I was doing some testing to make sure it works. In Github, I forked the master repository and created a working branch. My development version in QGIS has been this branch I created. The plugin works, exists in my QGIS plugins folder, and is called "eco-valuator-CreateAndExportPrintLayout", as this is the name of the GitHub branch I was working in :

I merged my development branch with the master branch (link above). Now if I download this repository from the master branch, the folder is called "eco-valuator-master". If I put this folder in with the rest of my plugin folders, QGIS doesn't like it. The plugin is listed as "Installed" if I go to "Manage/Install Plugins" but doesn't show in the processing toolbar (where this plugin should be) The photo below shows the one in question, highlighted:

Lastly, if I manually change the name of the folder from "eco-valuator-master" to "eco-valuator-CreateAndExportPrintLayout", the plugin works. So obviously the issue is with the folder name.

However, I can't find any file within that folder that contains the folder's name, so why is this affecting the way the plugin works in QGIS?

When looking at the config file, it lists the name of the folder will be "ecosystem_service_valuator". Does this impact what is happening? Maybe after I upload the plugin to the QGIS plugins repository, this will be a non-issue?


texts: Print with different settings of layer view?

texts:
Using QGIS is it possible to print on the same paper two windows with different settings of layer?

I mean:

window1: layer1 ON - layer2 OFF
window2: layer1 OFF - layer2 ON


texts: Difference of Mapzen Vector Tiles (mvt) to pbf extracted from .mbtiles

texts:
I was wondering what is the difference between a mvt tile and a pbf extracted from a .mbtiles container using mbutil?

Mvt is in pbf as well. The pbf is supposed to be a vector tile as well.

They are not exactly the same, since I tried to read the pbf like I read a mvt file, but I was not successful. Could anyone explain the difference to me?


texts: Trimming Specific Text Values in QGIS Field Calculator

texts:
Sorry if this is a dumb question, but I am wondering if there is a way to copy the text from one field to a new field, but eliminate specific words.  For example:

**Old Field**                         **New Field**

Killbear Provincial Park               Killbear
Quetico Provincial Park                Quetico
La Cloche Provincial Park              La Cloche
Wildcat Hill Provincial Park           Wildcat Hill

I thought of trying a trim function, but not all the values have the same start location or trim length.

Thanks...


texts: Extracting feature classes from file geodatabase to convert to CAD?

texts:
I am trying to convert multiple feature classes located within a file geodatabase into a CAD friendly format.

Without extracting the files the "Export to CAD" tool does not recognize the feature class.


texts: How to georeference an image dynamically

texts:
I am looking for dynamic approach how to georeference an image of a radar. This is a data that is going to be fetch every hour. Of course I will only need radar data (without base layer and borders) which means I will have to extract it, but before that (as I am going to move marks and lables of cities) I would like to georeference it. Image looks like this:

Main cities and smaller cities are marked with diamonds (bigger for main cities/smaller for smaller cities).

In general, what I am able to find out (for example over google maps) are coordinates of this cities as points of interests. I assume this task could be somehow done with GDAL or/and python but I am pretty much lost since I only use PHP and for this kind of tasks until now I only used Desktop GIS since I havent had need to use it for dynamic tasks.

Any help will be preety much appriciated!


texts: Bulk changing symbology using ArcGIS Pro

texts:
I'm currently processing a land classification map which is in feature class format. This map is symbolised by a defined color for each land class (yellow for sand, green for vegetation and blue for water).
Every time I conduct a process, I need to re-assign the symbology to the output.
Is there anyway I can re-apply the same colors from one feature class to another without having to manually selecting each colour one-by-one?


texts: Creating a subset of points shapefile with minimum distance between points in ArcGIS 10

texts:
I have a sample of ~700 animal locations that I will be using for habitat selection analysis. I am trying to figure out a way to subsample my locations in a way that will maximize my sample size for analysis, but minimize spatial autocorrelation. I envision taking a subsample of the points such that no point is within X distance (the distance at which autocorrelation is no longer a problem) of another. I don't necessarily need to select the "maximum" number of points, it is more of a conceptual problem of balancing sample size with autocorrelation. Is there a way to do this in ArcGIS 10? The images below show the type of data and spatial clustering I am working with. The autocorrelation distance is roughly equal to that of the gray buffer (570 m), so subsampling based on that will severely limit my sample size.  Is there a better approach that I am not thinking of?

Thanks for any input!


texts: How do you remove isolated roads from a bigger network?

texts:
I'm trying to perform an OD Cost Matrix calculation on a set of points using a road network from OpenStreetMap data. The network contains lots of roads that don't connect to the rest of the network. In the case where a point snaps to such a place on the network, it can't reach any other sample points.

How do I remove these sections of road? Or if that's not possible, how could I systematically identify those areas and connect them to the rest of the network?


texts:
Queen and rook neighborhoods are two common ways to calculate statistics for a focal cell. They are also known an Moore and (von) Neumann neighborhoods. I think a "bishop version", i.e. from lower left to upper right corner also exists. But why are these two methods still implemented in software packages?

Source: Lloyd, C. (2010). Spatial data analysis: an introduction for GIS users. Oxford university press.

I always thought that rook was used in the "old" days when computation power was low and that queen´s case should be used in GIS analysis wherever possible. I wonder whether there is any advantage of rook vs queen?


texts: Using .NET Runtime SDK or .NET ArcObjects?

texts:
I want to build a desktop GIS application with following features in it:

It can be able to read all the attributes of any digitized map
It can be able to draw polygon or line
With the help of that line, that app can calculate the distance between two points
it can be able to calculate area in case of polygon

I have two options and I do not know which to choose or what is best.

Should I use NET Runtime SDK or .NET ArcObjects to develop my GIS application?

My level of C# is moderate. I am not a beginner nor expert.


texts: I want to have the sum of at least 2000 raster combined to one raster layer using ModelBuilder (Plus tool). I think I'm missing something here. This model does not work.


texts: Creating ArcGIS add-ins for both 10.0 and 10.1

texts:
Has anyone tried to deploy the ArcMap add-in for both ArcGIS 10.0 and ArcGIS 10.1? There is a target version parameter in the add-in config XML file that needs to be specified. After upgrading the code to VB MS VisualStudio 2010 and ArcObjects 10.1 it can no longer be used in ArcGIS 10.0.
What are your strategies for upgrading your add-ins from ArcGIS 10.0 to 10.1 while maintaining the compatibility with 10.0?


texts: ArcGIS rest API error The given key was not present in the dictionary

texts:
I'm using the rest API to add a feature to a public feature layer and it's giving me the error "Error: 1000, Description: The given key was not present in the dictionary". What does this mean? How do I fix?

My feature service: http://services2.arcgis.com/lfD49qylDsb6IIcE/arcgis/rest/services/PropertiesWordpressLoading/FeatureServer/0

I add my feature using this form: /services/PropertiesWordpressLoading/FeatureServer/0/applyEdits

My json:
[
  {
    "geometry" : {"x" : -86.288146369965, "y" : 32.842269287433},
    "attributes" : {
      "CATEGORY" : 'For Sale',
      "LIST_ID" : '00-000',
      "CMN_NAME" : 'Toly Test Tract'
    }
  }
]

I get the same result if I use the form located in http://services2.arcgis.com/lfD49qylDsb6IIcE/arcgis/rest/services/PropertiesWordpressLoading/FeatureServer/0/addFeatures


texts:
I'd like to negate some values from my map using the NOT operator.

However, i need to deselect a whole field which contains over 600 entries so I don't want to select each. Is there a way to do this?

Is it possible to use: NOT SPA_NAME, without specifying the values

Added from answer:

I have a merged layer with many different soil and land use properties. And I'm trying to select areas that meet certain soil and land use conditions.

So far my expression is: SOIL_TEXTURE = 'clay' OR SOIL_TEXTURE = 'loam' AND LAND_USE = 'arable' AND LAND_USE = 'grassland'

So now I'd like to make sure special protected areas (SPA_NAME field) are excluded from the selection, but there are so many values that is wouldn't make sense to individually select each.

On closer look I am seeing that the command SPA_NAME IS NULL didn't work. Areas that are under special protection are still being selected


texts: Normalized Difference Built Up Indice with a threshold value in QGIS modeler

texts:
I want to use the modeler in QGIS, model the NDBI with a threshold value in order to produce a binary image. positive values = 1 and negative values = 0
my question is about which tool I can use in the modeler to create the binary image?


texts:
I have a map with several moving markers on it and have been trying to remove each individual one when their respective animation ends. I am using Leaflet and Leaflet's plugin Leaflet.MovingMarker. I tried using methods isEnded() and end without much sucess.

With the latter method (end), the following code removes only but one marker (full code here):

marker.on('end', function() {
  map.removeLayer(this);
  });

With a click event instead of an end event, each marker is made easily removable, but then this is not what I am looking for:

marker.on('click', function() {
  map.removeLayer(this);
  });

With the method isEnded(), have not had any luck at all either (full code here):

var ended = marker.isEnded();

  if (ended = true) {
  console.log("this marker ended!");
  map.removeLayer(this);
  }

I wonder what I am missing here?


texts: Allign shapefile polygons and underlying image in MATLAB spatial plots?

texts:
I am plotting a map in MATLAB. For this I have values at corresponding Latitude-Longitude.

I am creating meshgrids for Lat, Lon, and values for plotting the map.

Somehow plotted map does not exactly follow the polygon boundaries whereas polygon boundaries are accurate (I have verified it).

How to resolve this error?

Script and image are shown below, I am using geographic coordinate system.

axesm('MapProjection','mercator','AngleUnits', 'degrees');

geoimg = geoshow(lat,lon,data_plot, 'DisplayType', 'texture');
% lat lon, and data_plot are meshgrids
geoimg.AlphaDataMapping = 'none';
geoimg.FaceAlpha = 'texture';
alpha(geoimg,double(~isnan(data_plot)))
colormap('winter');
hold on

shapefile='area.shp'
shape = shaperead(shapefile,  'UseGeoCoords', true);
geoshow(shape, 'FaceAlpha', 0)


texts: Using Photo points to "geotag" old images in QGIS

texts:
I have a bunch of old aerial images(in jpg format), however they do not contain any coordinates.
I have a shape file which has the coordinates of the images and one of the attributes is the image file name. ( i assume that it will be the center point of these images)
Is there a way to "geotag" or to pin the photos to these points using QGIS?
There are about 6000 aerial photos.


texts: Transforming to WKTLine using JSTS library is very slow

texts:
I have a performance issue with JSTS library. When my OpenLayers.Geometry.LineString has more than 1000 components (OpenLayers.Geometry.Point), transforming theses points with the JSTS library is very slow (often it kills the browser)

    function TransformLine(feature) {
        var line = feature;
        var OLParser = new OpenLayers.Format.WKT();
        var parsedLine = OLParser.write(line);
        var WKTReader = new jsts.io.WKTReader();
        WKTLine = WKTReader.read(parsedLine);
        var WKTtoOLParser = new jsts.io.OpenLayersParser();
        WKTLine = WKTLine.buffer(20000, 2, 3)  // <== this line is very slow

        // ... drawing on map, skipped
}

Is there any way to boost perfomance of this operation?


texts: Extending line to feature programmatically

texts:
I'm trying to programmatically extend hundreds of lines (blue in image) to a single polyline feature (red in image). Sometimes the lines intersect each other, but I'd like the lines to extend, not to other existing lines, but to the polyline feature mentioned before. Even if i select 1 line and the polyline feature, the 1 line is stopped at the other unselected lines.
~ arcpy.ExtendLine_edit(lines, "2000 Meters", "FEATURE")

I can accomplish the process manually in edit mode by selecting the feature, clicking the extend line tool in the advanced editor toolbox, then clicking the line that needs extending. Is there a way to do this part of the selection in python?


texts: onSelect event on dynamically created layers

texts:
I'm creating different layers dynamically full of features, everything on-the-fly, and I need to do some stuff when each specific feature is selected. I can't paste the actual source code here as it has dependencies and it's quite long also, but you can imagine something like the following:

if (e.Bubbles) {
     jQuery.each(e.Bubbles, function (i, obj) {
         jQuery.each(obj, function (n, el) {

             /* features and layers ... */
             var point = new OpenLayers.Geometry.Point(el[1],el[0]);

             /* some code here ... */
         });

         layer = new OpenLayers.Layer.Vector(i, {styleMap: e.Style});
         map.addLayer(layer);
         e.Layers.push(layer);

         /* etc. */

When I'm trying to bind the onSelect event to the layer, it works but it replaces the previous layer behavior, so the last layer in the loop works fine when you clicking on the features, but only the last layer. Here is my code (inside the loop):

selectControl = new OpenLayers.Control.SelectFeature(
    Layer,
    {
        clickout: false,
        toggle: false,
        multiple: false,
        hover: false,
     }
 );

selectControl.activate();

layer.events.on({
    "featureselected": function(event) {
    console.log(event.feature);
    }
});

So, the selectControl will be replaced inside the loop each time and it will lose its own functionality for the previous ones, and just works for the latest layer.

I tried pass all the layers as an array as well (which is supposed to work) but without any success.

selectControl = new OpenLayers.Control.SelectFeature(
    e.Layers,
    {
        clickout: false,
        toggle: false,
        multiple: false,
        hover: false,
        onSelect: function(event) {
            console.log(event.feature);
        },
    }
);

selectControl.activate();

I hope somebody has been solve the problem!


texts: Both dashed line and stroke line in parallel in the SLD

texts:
I just started in GIS. How can I draw in the GeoServer's SLD editor, both dashed and stroke line in parallel?
like below:

-----------------------------------
___________________________________


texts: Building footprints of Atlantic city, NJ

texts: does anyone know where to get information of building footprints information on Atlantic city?


texts: Finding minimum edge to edge distance of polygons using ArcGIS Desktop?

texts: I have a a map of about 3000 polygons in ArcGIS 10. I'm looking to find the distance between each one of them. I know how to do it using the lat and long coordinates of the centroid, but I'm looking for the shortest straight line distance from the closest edge of one polygon to the closest edge of the other polygon. Any ideas?


texts: Query to find all lines connected to other

texts:
I have problem with one spatial (Oracle) query. Simplified overview of the problem is:

I have two tables:

CREATE TABLE thin_lines (
  id NUMBER PRIMARY KEY,
  geom SDO_GEOMETRY),
  color VARCHAR2(32);

CREATE TABLE thick_lines (
  id NUMBER PRIMARY KEY,
  geom SDO_GEOMETRY),
  color VARCHAR2(32);

Both tables contain 2D lines. The example picture:

The problem:

I want to create SQL query or PL/SQL procedure to find ids of all blue thin lines (whole blue branch comprising of blue lines) connected to the thick green lines.


texts: How to use buffer() in QGIS Value Relation field

texts:
When adding a point, I'd like to auto-populate one of its fields (the street address) setting it to a field in the nearest point in another layer (a layer containing all existing addresses).

It looks like the Value Relation edit widget could do the trick, if I could figure out how to use the buffer() function.

Right now, I can get the list of all addresses by selecting the Layer, Key column and Value column. I thought I could filter that list so that it has only nearby points (ideally just one), using something like:

within($geometry_A, buffer($geometry_B,50)) = 1

where $geometry_A is the new point and $geometry_B is the existing point.

$x, $y and $geometry seem to refer to the existing points in the selected Layer.

Assuming I am on the right track, how do I refer to the geometry of the new point (i.e. $geometry_A)?


texts:
I have a layer with some geometries and a data layer that defines an M:N relation. Is there any way I could select (on the map) objects that are in relation with specific id? Basically I would like to select by expression like

object_id in (select object_id from relations where relative_id = 12345)

I can do that from the db manager select * from objects where object_id in..., but that means creating new layer each time, which is extra work and needs to be deleted or there would soon be hundreds of them. And the qgis joins don't seem to support M:N relations. Is there a way to select (in qgis, not sql sense) by such query?

The layers happen to be in spatialite database. They look something like

create table objects (object_id integer primary key, geometry, attributes...);
create table relations (object_id integer references objects (object_id),
                        relative_id integer references something);


texts: Is it possible to catch mouse events from a hidden feature layer?

texts: I want to avoid drawing a feature layer. I want it to be hidden and load as fast as possible. I have tried with opacity zero but i think it still draws.  Is there any way to minimize the time of being put on the map? I care only for mouse events at the position of the geometries.  I need this so as to to keep the feature layer in a map containing the same poins that has been designed with a dynamic or tiles service for better performance. So i will have my points drawn as images altogether and then add the feature layer to as t catch the mouse over event and show a popup .


texts: How do I modify features during runtime, in a GeoJSON feature group on Leaflet?

texts:
I have added some GeoJSON to my Leaflet powered map as below. (Altered version)

var array_of_layers = L.layerGroup([]);
var styleholder = [];

var belgium_style = {color: "#d3e46f"};
var bosnia_style = {color: "#f3c1d3"};
var (...)

function styling(name, geoname, color){
    this.numid = count;
    styleholder[count] = this;

    this.name = name;
    this.geoname = geoname; //the geojson filename, just in case.
    this.label = name;
    this.color = color;
    this.opacity = 1;
    this.stroke = false;
    this.fillOpacity = 1;
    this.fillColor = color;

    count++;
}

var belgium_style = new styling("Belgium", "belgium", "#52EB8D");
var bosnia_style = new styling("Bosnia", "bosnia", "#EB9752");
(...)

var everything = L.geoJson(myfile, {
    onEachFeature: function(feature){
        array_of_layers.addLayer(feature);
    },
    style: function(feature){
            switch(feature.properties.name){
            case "belgium": return belgium_style; break;
            case "bosnia": return bosnia_style; break;
            case "denmark": return denmark_style; break;
            case "great_britain": return britain_style; break;
            case "greece": return greece_style; break;
            case "italy": return italy_style; break;
            case "serbia": return serbia_style; break;
            case "spain": return spain_style; break;
            }
    }

});

So the next thing I need to do is, for instance, make a country become blue and all the others gray, when I mouse over it. And, make that possible as well by a function, that I can call whenever I want, regardless of mouse position.

Therefore, during runtime, how do I change the color of any country as I wish? In other words, what is the correct syntax for re-defining the style function for a particular feature?

I have tried to iterate over the styleholder array with a for(var x in styleholder) loop until I get to find the style with the name I want to make blue. By looking at the debugging console, I can tell that I get to find the country, but there's no use changing the color value of a country once it's been created on the map. So the color doesn't change.

The code below sums up what I'm talking about:

var change_color_to_blue = function(country){

for(var a in styleholder){
   if (a.geoname === country){
       a.color = "#0000FF";
       a.fillColor = "#0000FF";
   }
}

change_color_to_blue("belgium");

You see, although it runs without an error, the color remains the same.

The same thing for a mousing-over feature, could be possibly something like this:

everything.on('mouseover', function(e){

var moused_over_country = e.layer;
moused_over_country.color = "#0000FF"; //this is not how it's done

MagicallyMakeAllOtherLayersGray();

refreshgeoJson(); //that does not exist.
refreshMap();     //this neither, most probably.

});

Even this works, and the target country gets blue, how do I make the rest gray?
Is there a way to iterate over every feature in geoJson and reset every style function respectively?


texts: How to set position of label around a point in qgis according to a field

texts:
I'm working on a shape file in qgis (2.8.4 Wien for Mac). In this file I have several points with about ten fields. What I'm interest in is to use two of these fields to set position of the label and the rotation of the symbol of the point: one field is the id of the point, the other one is the rotation.

I represent the point symbol as a triangle with a little line, and I rotate it according to the rotation field (I have no problem with this).
The main problem is to place the label of the ID around the point in a specific position, according to the rotation, but without rotate it.

The rotation field contains, obviously, values from 0 to 359.

In the picture below, there is a sample of what I want to obtain:

I think is possible to use a offset function with sinus and cosinus, but the corresponding buttons are disabled on qgis interface.

There is a way to do this? Unfortunately this is a requirement of the work I'm doing.


texts: QGIS GPS Connection with PyQgis causes minidump

texts:
After some experimenting and reading source codes I tried to connect to my GPS device as follows (from a QGIS plugin):

self.pushButton_conGPS.clicked.connect(self.connectGPS)

[...]
def connectGPS(self):
    port = u'\\\\.\\COM4'
    #port = u'localhost:2947:'
    self.detector = QgsGPSDetector(port)
    print 'detected, connection functions'
    self.detector.detected[QgsGPSConnection].connect(self.connectionSucceed)
    print '---'
    self.detector.detectionFailed.connect(self.connectionFailed)
    print 'advance'
    self.detector.advance()

def connectionSucceed(self, connection):
    print 'connection Succeed'
    QgsGPSConnectionRegistry.instance().registerConnection(connection)
    print QgsGPSConnectionRegistry.instance().connectionList()[0].currentGPSInformation().utcDateTime.toString()

def connectionFailed(self):
    print 'connection failed.'

Pressing my above referred button results in

Note the output on the python console, minidump occurs at QgsGPSConnectionRegistry.instance().registerConnection(connection)

Related post: Connect / Disconnect GPS device via PyQGIS

What is it that causes the minidump? I think similiar c++ code of the GPS Information panel works pretty well...


texts:
I have a minLat, minLong, maxLat and maxLong from a box drawn on a map. Right now, my query to get all the points in the box looks like this:

SELECT *
FROM geomTable
WHERE (longitude BETWEEN minLon AND maxLon)
AND (latitude BETWEEN minLat AND maxLat)

I want to use the geometry column instead of lat/long columns to get the results.
 I tried this query:

SELECT *
FROM mytable
WHERE mytable.geom && ST_MakeEnvelope(minLon, minLat, maxLon, maxLat, 4326);

from this post: Select bounding box using postgis but it is returning no results.

Does anyone have an example of how to select all the points within a box created by min and max lat/longs using geometry in postgis?


texts: Avoiding Feature limit exceeded on layer in ArcGIS JavaScript API?

texts:
My code as below:

require([
    "esri/Map",
    "esri/views/SceneView",
    "esri/layers/FeatureLayer",
    "dojo/domReady!"
], function (Map, SceneView, FeatureLayer) {

var featureLayer = new FeatureLayer("https://services.arcgis.com/V6ZHFr6zdgNZuVG0/ArcGIS/rest/services/votes_popup_final3/FeatureServer/0");
var map = new Map({
    basemap: "streets",
    ground: "world-elevation",
    layers: [featureLayer]
});
var view = new SceneView({
    container: "viewDiv",
    map: map,
    scale: 5000000,
    center: [-110, 32]
});

});

The result as below:

Only show a little part of all features. And show the below error in Firebug:

[esri.layers.graphics.controllers.SnapshotController] Feature limit exceeded on layer " Votes popup final3 - All Votes ". Not all features are shown.

But when use MapView, it can show all features. The code as below:

require([
  "esri/Map",
  "esri/views/MapView",
  "esri/layers/FeatureLayer",
  "esri/widgets/Legend",
  "esri/widgets/LayerList",
  "dojo/domReady!"
], function (Map, MapView, FeatureLayer, Legend, LayerList) {

var featureLayer = new FeatureLayer("https://services.arcgis.com/V6ZHFr6zdgNZuVG0/ArcGIS/rest/services/votes_popup_final3/FeatureServer/0");
var map = new Map({
    basemap: "streets",
    layers: [featureLayer]
});

var view = new MapView({
    container: "viewDiv",
    map: map,
    zoom: 4,
    center: [-110, 32]
});

});

Of course need include webgl as below:

var dojoConfig = {
  has: {
    "esri-featurelayer-webgl": 1
  }
};

This time it can show all of features. The page shows as below:


texts: Combining features with multiple similar field values to one and adding unique values as new fields using ArcGIS Pro?


texts:
How can I create a custom projection with the following parameters as per ArcGIS Desktop?

PROJCS["Lambert",GEOGCS["GCS_NZGD_2000",DATUM["D_NZGD_2000",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Azimuthal_Equal_Area"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",175.0],PARAMETER["Latitude_Of_Origin",-35.0],UNIT["Meter",1.0]]

I need to create one in QGIS.


texts: GeoServer WCS geotiff style not loading in ArcMAP

texts:
I have published a GeoTIFF file using GeoServer's (2.4.5) WCS service.

I also defined a colormap style on GeoServer.  I know the style works, because, when I load the map through WMS the colormap is there and the legend image is also generated.

But when I try to load it into ArcMAP or QGis using WCS, I just get a stretched gray scale.

Opening the same GeoTIFF file directly from the catalog in arcmap, lets me calculate the statistics and histogram, and then I can define and edit a colormap for that layer, but this option to calculate the statistics for the WCS raster is not available.

I also tried to create a classified or unique value colormap in the symbology tab at the properties dialog, but those options are not available for the WCS, only for the file opened directly.

How can I make ArcMap take into account the style information I created using GeoServer or, alternatively, how can I make ArcMap let me create and use a classified colormap for this WCS raster.

Thanks to anyone who tries to help.


texts: How can i keep browser from caching map services?

texts: I am using esri's js api with arcserver 10.1 map services. My issue is that sometimes another application deletes a feature from the geodatabase that is the basis for my map service. But my browser does not react quickly to the change. It seems to be caching the map service response. This is evident in ie and firefox. Even a browser refresh will still show the feature that should have been deleted already. Is there a way to force the cache to flush and get new map service data?


texts: Checking layers with scale based visibility for visibility in QGIS

texts:
I am writing a Python method which I want the mouse click event to connect to. The method selects objects from layers that satisfy certain conditions. I want to select from any visible layer, that has features in the current map extent.

To evaluate which layer is visible I check the property:

QgsProject.instance().layerTreeRoot().findLayer(layer.id()).isVisible()

There are some layers with scale dependent visibility. I don't want to select objects from layers, that are not visible due to the current map scale.

isVisible() from the line of code above returns True, even if the layer is not visible due to the current scale, but is checked in the legend tree.

My question: is there another way to check if a layer with scale based visibility is out of the visible scale interval, then comparing the map scale with minimumScale() and maximumScale() of each layer?


texts: Gravity/Huff model tools in QGIS without rasterlang plugin?

texts:
As the person of this question:

Gravity/Huff model tools

I am looking for a tool to simulate a gravity model using a point-based layer.

As in the case there, my points have a z-value and the higher this value, the larger is their "sphere of influence". Their influence is inversely proportional to the distance to the center.

The answer given there by jake would be perfect, but unfortunately it requires the rasterlang plugin which, as far as I understand, is not supported anymore (or at least it does not work on my QGIS 2.8.2)

So.. does anyone have any other idea?


texts:
I'm adding some functionality to an existing QGIS plugin (EcoValuator from Key-Log Economics). After playing with the existing version of the plugin, it works on a PC, but not on a Mac. The issue is with a step in the code which calculates values in an attribute table.

Before diving too deeply into the code, is there a reason why the plugin would work on PC but not on Mac?  What are some things I should look for?

I was unable to find mention of this on StackOverflow, QGIS documentation, or elsewhere.


texts: PyQGIS select layer by original table name

texts:
As far as I see, on QgsMapLayerRegistry.instance() there is a method to select layers:

mapLayersByName(Qstring)

and the method

mapLayer(QString)

respectively, wich accesses a layer by its id.

As the name says, mapLayersByName() selects layers by their name in the legend, wich causes scripts or plugins to fail if a user renames a layer/table (i.e. renaming a table '160520_parcels' to 'Parcels').

So is there a way to access a layer by the name of it's origining table in PyQGIS?

update

something like ...

def selectLayerByOriginalName(tablename):
    # return layers whose orininal names start with 'tablename'
    llist = []
    for l in QgsMapLayerRegistry.instance().mapLayers().values():
        if tablename in l.id():
           if l.id().index(tablename) == 0:
               llist.append(lyr)
    return llist

... works, but is there a more convenient approach?


texts: Configure shortest path calculation with specific origin/destination points

texts:
I would like to calculate the shortest path from a point to another following a network. I'm currently using QNEAT3 for that purpose, but I can either calculate one path from two specific points or a full OD matrix (so based on all unique origin/destination combinations). What I'm trying to do is set specific OD couples (for example by first calculating a distance matrix with the 5 closest points) then calculate the distance between those combinations (and not all combinations, it is much too heavy for my computer).

Does anyone know how to handle this?


texts:
I am trying to find digitizing errors in a line-layer. Obviously the snapping was not activated when the person who digitized created the lines.

As far as I know QGIS cannot find undershoots, so I try to find close vertices.
OpenJump has a tool called "Find Close Vertices" which helps finding some of the digitizing-errors. How could I archieve something like this in QGIS?

I tried to use v.clean with the break option to find the errors but this did not help much.

Test-data can be found here: https://gist.github.com/rdbath/104876de54fdfe0fc3ceed521220673e?short_path=f75c9b4


texts: Export all visible layers in map view window to DXF (as one)

texts: I'm looking to find a way of exporting all that is visible in the map view window to a dxf file. That is to say, multiple features from multiple layers, all at once. Can this be done?


texts: Timeout while uploading large files in GeoNode / GeoServer

texts:
I am getting timeout while uploading large files in GeoNode ( System Configuration : GeoNode 3.1.0 , GeoServer 2.16.2 installed on Ubuntu 18.04. ).
I tried compressing large files to decrease their size (E.g. gdal_translate -of GTiff -co "COMPRESS=DEFLATE" -co "PREDICTOR=2" uncompressed.tiff compressed.tiff ), this helped for some cases only. But I have other files which are still large.
So I wish to increase the timeout. But, I could not find some exact reference which can guide me to increase the timeout for this scenario.
While exploring geonode/geoserver settings, I found different timeout settings and tried to increase timeout period in different settings as follows:

Timeout increased in /etc/nginx/nginx.conf

keepalive_timeout increased to 300 from 65

Timeout increased in GeoNode UWSGI config
(/etc/uwsgi/apps-enabled/geonode.ini)

env = OGC_REQUEST_TIMEOUT increaded to 600 from 60.

harakiri increased to 600 from 300

Still timeout occurs while uploading large files. Is there some settings I am missing?


texts: What is a relevant geospatial modelling approach when no spatial autocorrelation is found?

texts:
Although I have found on StackOverflow all answers I ever had on GIS matters using R, ArcGIS and QGIS in the past, I am currently faced with one question I can't find an answer to online.

For my final dissertation project, I am working on water quality data based on sampled I have collected across an aquifer. I am exploring faecal contamination patterns, so I am mostly interested in faecal bacteria count and my method (using R) is as follows:

1. ESDA

2. Stepwise logistic regression for dimensionality reduction. That's because my dataset is very wide given that I only have 97 samples, and I am trying to predict a contaminated/non-contaminated status, not a level of contamination (hence logistic over linear regression)
--> the best possible result is still a rather poor model with low predictive power and no spatial autocorrelation found in residuals (that excludes using a Geographically weighted Regression)

3. Geospatial modelling: still trying to run some form of interpolation
--> I tried Moran's I, LISA, and even Mantel test (distance matrix against faecal contamination matrix) and found no spatial autocorrelation for the variable of interest.
--> I used IDW, because unlike Kriging it doesn't make assumptions about spatial autocorrelation

4. Hierarchical clustering (unsupervised machine learning) to still extract some form of classification from this complex dataset despite the difficulty to predict anything with supervised methods tested above.

My main question right now is about step 3. Given that I found NO spatial autocorrelation, am I right to still go ahead and run an IDW interpolation? Or should I just stop my analysis right there and say: look, my result is "no significant result from a geospatial point of view"? I have read so many articles, tutorials and handbooks, that this useful question Choosing IDW vs Kriging interpolation summarizes pretty well, but it doesn't really lift my doubts.

@Spacedman has provided a very useful answer.

My sampling pattern was indeed quite poor due to more or less limited access to groundwater sources across the study area.

Initially I was testing a new real-time method for faecal matter detection, but it turns out it just doesn't work in this area and is shockingly uncorrelated with the actual contamination. That's what I am finding in step 1. So my study has become an exploration of contamination patterns across the aquifer. In step 2 I try to work out whether other hydrochemical parameters I registered would work as decent predictors and if I can build a regression model. Turns out: not really. So in step 3 I was hoping that I could still build some sort of classification of the study area and display contaminated vs non-contaminated areas.
In step 4 I rely on Machine Learning, and I use Gower distance (which combines all parameters into a single distance metric) to classify different "families" of contamination, if that make sense?

Could @Spacedman clarify what he means by "you have demonstrated that a non-spatially correlated model is adequate"?

Does it mean I should just stick with my logistic regression or the clustering?


texts: Creating polyline segments by points in maximum distance to each other along another polyline using ArcGIS Desktop

texts:
I have a polyline and point layer. The polyline layer contains different roads. The point layer contains points along these roads with random distances between the points. The points are also listed randomly and not in order along these roads.

I want to connect with polylines the points along the roads which are close to each other at a maximum (e.g. 35.0 m) or smaller distance.

The roads turn even to 180° at some locations (as shown in the attached figure) and the distance between the roads can be smaller than 35.0 m.

My ultimate purpose is to get the road segments where the points are distant in maximum 35m.

Any suggestion how to work these out?

I use ArcGIS 10.4.1 for Desktop


texts: Can qgis_core work with external python distribution (such as Anaconda)?

texts:
I've built a python application that has its own GUI. I am trying to use qgis.core/PyQGIS as map display only. Specifically, I have anaconda 2.3 and its python distribution (2.7.11).

After research on postings on this thread, it seems that this scheme can only work if the application I am building (or anyone else' application) would need to "live inside QGIS world", i.e. move all dependencies to QGIS' python distribution. Could somebody (preferably from QGIS team) to give a definitive answer to this.


texts: Is there anything faster than PostGIS ST_Intersects for seeing where features intersect (or something wrong with my approach)?

texts:
Basically have 2 tables.  1 is buildings (1809555 rows called oregon).  Another is flood zone (164 rows called oregonslr1).  Want to see where they intersect to get number of buildings flooded.  I am doing:

SELECT COUNT(*)
FROM oregon, oregonslr1
WHERE ST_Intersects(oregonslr1.geom, oregon.wkb_geometry);

I want a count of the geometry in oregon that intersects with the oregonslr1 geom.  Problem is it is taking forever.  Ran for like 90 minutes got nowhere.  The reason I switched to PostGIS was because this data was too big to work with in QGIS.  I'm wondering if I am going about this wrong or need to make so optimizations because I may have to scale this for many more states and it would be incredibly time consuming at the rate it is going? Also using pgAdmin, not sure if makes slower....

Added results from EXPLAIN, couldn't run explain analyze too long

Subdivided stuff with this command:

create table oregonslr2_sub as select st_subdivide(geom, 25) as geom from oregonslr2;


texts: Running ArcGIS RunTime Map With C# Code on WPF?

texts:
I have a WPF application called WpfApplication1. In ManiWindow.xaml I have a

Snippet

<esri:MapView x:Name="MyMapView">

</esri:MapView>

and in ManiWindow.xmal.cs I only have

namespace WpfApplication1
{
    /// <summary>
    /// Interaction logic for MainWindow.xaml
    /// </summary>
    public partial class MainWindow : Window
    {
        public MainWindow()
        {

            InitializeComponent();

        }

    }
}

According to ESRI we can add the Map to the MapView as as long as we have an empty MapView has been added to the page at design time and has been given the name MyMapView (<esri:MapView x:Name="MyMapView"/>.

Now my question is how/ where should I add this code to get map running?

// create a new Map
var myMap = new Esri.ArcGISRuntime.Controls.Map();

// create a new layer (world street map tiled layer)
var uri = new Uri("http://services.arcgisonline.com/arcgis/rest/services/World_Street_Map/MapServer");
var baseLayer = new Esri.ArcGISRuntime.Layers.ArcGISTiledMapServiceLayer(uri);
// (give the layer an ID so it can be found later)
baseLayer.ID = "BaseMap";

// add the layer to the Map
myMap.Layers.Add(baseLayer);

// set the initial view point
var mapPoint = new Esri.ArcGISRuntime.Geometry.MapPoint(-117.445, 33.9,
    Esri.ArcGISRuntime.Geometry.SpatialReferences.Wgs84);
var initViewPoint = new Esri.ArcGISRuntime.Controls.ViewpointCenter(mapPoint, 250000);

myMap.InitialViewpoint = initViewPoint;

// assign the Map to the MapView's Map property
this.MyMapView.Map = myMap;


texts:
I have one string with points:

var points =
  [-417494.99760966154,4921341.910474674,-416272.0051570987,4922111.058071794]

And I am trying to print a Box

var Square = new ol.geom.Polygon([points]);
var SquareFeature = new ol.Feature(Square);

source = new ol.source.Vector();
source.addFeature(SquareFeature);

But it not working.

How should I parse the points?


texts: Converting Hebrew Excel file into Shapefile in QGIS

texts:
I have an Excel that was written in Hebrew that I need to convert into Shapefile file. I am working with Windoes 7 and with  QGIS 2.18.17.
When I was converting the Excel into CSV I used "CSV (Comma delimited" format.
I also took care to be sure that the numbers and the words in the Excel columns are in the correct encoding.

The result was an error with the layout of the Hebrew words.

I head tried other kinds of encoding for the Excel conversion and for the Delimited text file as well Including, but it did not help.


texts: Making data joined from Excel independent in ArcMap?

texts:
I have joined some data from Excel (.csv) "text" to my attributes table of one feature class. However I can't edit it. I found out it is still connected to Excel file from where I added it, and I need to edit changes in original Excel file to make them appear in ArcMap 10.2.

I'd like to know, how can I break the link from the original file and make the data appear independently in a 'field' inside the feature class.


texts: Can POSTGIS 2.0 be install with POSTGRES 9.2 on Ubuntu 12.04 OS at the moment?

texts:
This question is more about environment set up.

At the current moment, March 4th 2013, can POSTGIS2.0 be install with POSTGRES 9.2?

I check their website out and to my understanding it is not possible...

 http://trac.osgeo.org/postgis/wiki/UsersWikiInstall#PostGIS2.0

I hope that's not the case. Any one can tell and point out the instruction how to install POSTGIS 2.0 on POSTGRES 9.2 on Ubuntu?

Thanks


texts: Reducing the accuracy of coordinate data

texts:
Let's say that I have a large set of latitude, longitude pairs which are accurate to 5 meters. I want to reduce their accuracy to 100 meters so that they render quicker and consume less storage space. Is there a name for this sort of lossy compression of coordinate data?

For instance, the three yellow points in the image below are part of a larger polyline. The red line represents the polyline at 5 meter accuracy. The thick blue line represents what the polyline would look like at 100 meter accuracy and the thin blue line represents the margin of error on that thick blue line (+/- 100 meters). Because the middle point falls within that margin of error, it can be eliminated entirely, thus reducing the storage size of the polyline.


texts:
I am looking  for comprehensive address data (incl. house numbers) for Germany. Is there some kind of government agency that provides free downloads?

I know that the BEV (Bundesamt für Eich- und Vermessungswesen) provides this service for Austria.

Are other countries also that generous with their data?

As mentioned. I am particularly interested in a dataset for Germany, but also information about other countries.


texts: ogr2ogr from MIF to KML including MID generating split kmls

texts:
I have a set of MIF  and MID files that I would like to translate to KML.

I used ogr2ogr to transform the MIF into one big KML which worked ok like so:

ogr2ogr -f "KML" geo-kml geo-mif

I have two questions:

Can I incorporate the accompanied MID file into the KML?
Can I make seperate KMLS for each placemark in the MIF?


texts:
I'm making a Map Visualization which requires an overlay that shows the light/dark patterns across the planet (like this).  I accomplished this by building a simple image layer that loads an image with a cosine wave shaded black on the top half (I'm not worried about seasonal adjustments yet).  I can set the extent on the  according to the current time to get a map that shows the light patterns in real time.  However, I can't figure out how to move the Image Layer without moving the base Layer to show the light/dark overlay in the past or future.

My first inclination was to use some sort of move() method, but I couldn't find it in the API.  So I tried to adjust the Layer's extent on each change like this:

sunLayer({extent : new OpenLayers.Bounds(newLong, -90, newLong+360, 90)});
sunLayer.redraw();

But nothing happens.  I was thinking that maybe I need to build a Vector Layer with just one Feature (the image), and then tie the dimensions of the graphic to the dimensions of the map and move the feature around, but this approach feels wrong to me.

Has anyone else done something like this (where you entirely decouple an Image Layer from the base layer)?


texts: Draw arc in Geo map WGS84 coordinate system

texts:
im trying to draw an arc in OpenLayers, in non Geo map it works well and here is my code:

function calculateArcPointsNonGeo(center, radius,startAngle,endAngle,segments)
 {
            var pointList =  [],
                Angle,
                x,
                y,
                point,
                i;

            var dAngle = segments + 1;
            for (i = 0; i < dAngle; i++) {
                Angle = startAngle - (startAngle - endAngle) * i / (dAngle - 1);
                 x = center.X + radius * Math.cos(Angle * Math.PI / 180);
                 y = center.Y + radius * Math.sin(Angle * Math.PI / 180);
                point = new MapPoint(x,y);
                pointList.push(point);
            }

            return pointList;
        }

in addition there is a great Example that demonstrates it , but when im trying to do this on Geo map (WGS84 coordinate system) its start to acting weird, i assume that i need to convert degrees to radinas or something like this.
please help


texts: Creating coordinates from beg and end points on a line

texts:
Is it possible to create WGS-1984 coordinates from the start and end points along a line segment?

The line features were created using linear referencing but now I need to create lat and long coordinates for the start and end of each line segment.

Here is some of the data from the attribute table:

OBJECTID *  CORRIDOR_CODE_RB    FRFPOST TRFPOST Shape *
   9        C001976N               0      0.76  Polyline M
   10       C002068N             0.432    0.984 Polyline M
   11       C002111N               0      0.432 Polyline M


texts: Check if custom CRS already exists, find EPSG

texts:
I created a "custom" CRS with specifications I got from someone who needs the data in that specific system. It worked well and I now have a *.prj file, so basically a WKT-representation of that CRS. How can I find out if it already exists as some "official" CRS which might even have an EPSG code?

I tried this with Python / OSR:

from osgeo import osr
sr = osr.SpatialReference()
wkt = '''MY_WKT'''
sr.ImportFromWkt(wkt) # returns 0, so that worked
sr.AutoIdentifyEPSG() # returns 7, of which I don't know the meaning, but it does not seem to work
sr.GetAttrValue('AUTHORITY', 1) # returns '4326', but MY_CRS is not the standard WGS84, so that's definitely wrong

Does that mean that my CRS is "truly custom" since there is no EPSG for that?  Is there another way of finding out whether it already exists or not?


texts: Merging GeoTiffs with GDAL causes grid-like artifacts?

texts:
I'm processing ASTER DEM tiles, merging them into larger blocks (to minimize visible tile edges) and importing the contours into a postgis database. Workflow:

gdal_merge.py -of GTiff -o merged.tif tile1.tif tile2.tif
gdal_contour -i 10 -a ele -snodata -9999 merged.tif contours.shp
ogr2ogr -f PostgreSQL -nln srtm4 -append -t_srs 'EPSG:3785' contours.shp

For some reason the resulting "merged" contours have grid-like artifacts:

When I skip the merge step and process the same two GeoTiffs separately, there are no artifacts. I've tried merging with gdalbuildvrt and got the same grid effect.

Why is this happening and how can I avoid this?


texts: mod_tile not using all available threads

texts:
I'm setting up an OpenStreetMap server with mod_tile, renderd and mapnik. My problem its low performance. The machine is set up with 8 cores and 16GB RAM.

I have 2 use cases:

render_list -af -n 8 <coordinates> - renders like a wonder, all 8 cores working fully, getting the tiles is a snap
renderd -f -c <path to conf> and then requesting tiles via slippymap.html equivalent in browser, then served by mod_tiles - it renders with average speed on levels 15-17 but it's getting worse on lower levels and on 10th level it's really very slow and there are white holes where 404 was given instead of tile. What's strange: out of 8 cores only 1-3 are working and mostly not at 100%. Memory usage isn't a problem also, only <0.5GB out of 16GB are used.

As render_list is working quite good, I assume that mapnik, renderd and postgres' database are working OK, but mod_tiles lags. Here is /etc/apache2/mods_enabled/tile.conf:

ModTileTileDir /var/lib/mod_tile
LoadTileConfigFile /usr/local/etc/renderd.conf
ModTileRenderdSocketName /var/run/renderd/renderd.sock
ModTileEnableStats Off
ModTileBulkMode Off
ModTileRequestTimeout 0
ModTileMissingRequestTimeout 30
ModTileVeryOldThreshold 31536000000000

ModTileCacheDurationMax 604800
ModTileCacheDurationDirty 900
ModTileCacheDurationMinimum 10800
ModTileCacheDurationMediumZoom 13 86400
ModTileCacheDurationLowZoom 9 518400
ModTileCacheLastModifiedFactor 0.20
ModTileEnableTileThrottling Off

LogLevel debug

renderd outputs standard PNGs in 8x8 metatiles. If there is something else I should post, please let me know.

Is there anything I can do to get similar performance of mod_tile's serving as I have with render_list?


texts: ArcPy merge using a list from GetParameterAsText

texts:
I am pretty new to ArcPy so bear with me.  I am having a hard time understanding how to feed a list of FCs from a "GetParameterAsText" into something like the merge tool.  I have my GetParamterAsText window set to Feature Class and multivalue so I can drag a group of layers into it in the tool dialogue.

inputLayer = arcpy.GetParameterAsText(1)

When I run the merge command, it looks like it is messing up the expected syntax by putting ";" as a separator instead of ",".

Doing this:

arcpy.Merge_management(inputLayer, outputName, "")

Throws me this error:

"C:\Development\Merge_Dissolve\GA0combined.shp;C:\Development\Merge_Dissolve\GA25combined.shp" does not exist

I have tried to create a function to replace ";" with ",", but it doesn't seem to work on my inputLayer variable, but I'm guessing there's a more standardized way to working with a list of FCs with tools that accept multiple inputs such as Merge.


texts: I have two TIN surfaces, one being a pit shell and a second being a horizontal layer which intersects the shell on all sides. I need to find out the volume between the intersection horizontal layer and the base of the pit, but I'm having trouble clipping the layer to the pit edges in order to do so,


texts: QGIS 3.4.3 - Input dialog for a column in attribute table

texts:
I’m working with QGIS 3.4.3 and Python and I’m looking for a better mechanism in attribute table (vector layer) which tells me which fid and which index of a column I’m currently working with.

For example, double click on a column, please tell me fid and index. Similar to signal setAttributeValueChanged but in my case the response is too late. I need the information before entering a value.

My idea is to create a dialog with possible input values and descriptions.

Currently my mechanism works if the user enters value “?” in a column. I get all input values and descriptions from an external dll. This works fine.

I found out a way which could work. I have created an action "My help" using python code and have made visible it in attribute table. Fid and Index is shown correct.

Now I have two questions:

1) Can I create my action "My help" in python automatically. Currently I did it in QGIS interactively. I need a sample.

2) The action "My help" needs access to my python plugin. I need to call some functions. Is it possible? How? (Ok, problem solved by myself):

fid = [% $id %]
idx =  [% @field_index %]
aLayer = qgis.utils.iface.activeLayer()
myhelp_dlg = myhelp.classFactory(qgis.utils.iface)
value = myhelp_dlg.ExampleDialog(aLayer,int(fid),int(idx),None)
aLayer.changeAttributeValue(int(fid), int(idx), str(value) )
myhelp_dlg = None


texts: I need to convert a lot of garmin GDB files to GPX in order to convert it to another vector format. Is there any tools available can do a batch convert from garmin GDB files to GPX? Using MapSource allow only one file at a time.


texts: Purchasing NDVI imagery for South Western England?

texts:
I am working on a project developing an application for livestock farmers where they can monitor their feed using satellite data.

The main aim is to calculate the biomass of a field of grass.  The spectral resolution can be 10m but we are also considering better resolutions.  The temporal resolution should be high enough to account for cloud cover in the UK.  Where do you suggest starting?

The project is still in proposal stage so I need to consider all potential budgets.  We would ideally like to use only satellite data for commercial reasons.

I am new to satellite imaging and GIS so a bit of extra explanation might be necessary!


texts: Creating field indicating if two layers intersect

texts:
I have the similar problem to this question but in QGIS: Create a field indicating if the polyline intersects the shapefile boundaries
I want to select the areas intersecting with the layer station_300m I used this expression:
intersects($geometry, geometry('station_300m_5142f160_96d1_491b_836f_9bc56c884a03'))

and I got this error message

I also tried this expression:
intersects($geometry,geometry(get_feature('station_300m_5142f160_96d1_491b_836f_9bc56c884a03')))

but it also didn't work. What goes wrong in my expression?
I'm using QGIS version 3.16.6-Hannover.


texts: How to assess significant land cover changes using Chi-square test?

texts:
I have a series of land use maps at different times as a result of a classification process on Landsat images. I am mapping 3 different thematic classes : bare soils, vegetation and croplands.

I would like to answer the following question : "are land use change values significant ?" I believe chi-square test might help here.

How does the test work for my use case, and how should I run it ? Should i sample some pixels (if so, how much ?), or or work with the entire rasters ?


texts: Need python code example - setting a combo listbox to Null

texts:
I want to be able to set the value of a field to Null using python, but I get an error.
I am using arcpy.CalculateField_management


texts:
I'm aware that this site is intended for cartographers, geographers and GIS professionals. I must admit that I'm absolute beginner in the area of maps, but not in the area of programming (40+). I hope somebody will be kind enough to point me in the right direction.

I have map consisting of a set of 256x256 jpg orthophoto tiles. For each tile, I have WGS84 lat/long coordinates and info that 1 pixel is 0.5 m. First I would like to create local Leaflet map on my PC (concept from http://gis.stackexchange.com/questions/82936/how-i-can-load-tilelayer-in-leaflet-framework-using-local-tiles) and then also on Android (concept from https://android.stackexchange.com/questions/8312/how-can-i-can-open-local-files-in-the-default-android-browser).

How far did I manage to get? I found an online map which uses Leaflet: https://mr.si/topo/. One of the layers on this map (Topo) is using tiles the way I want to:
L.tileLayer("https://s3-eu-west-1.amazonaws.com/topo-slovenia/z{z}/{y}/{x}.png", {
            minZoom: 10,
            maxNativeZoom: 15,
            detectRetina: !0,
            attribution: '© GURS',
            unloadInvisibleTiles: !1
        })

Question for me is, of course, how do I name my folders and tiles to correspond to {y} and {x}?

I obtained one of the tiles from above map at address https://s3-eu-west-1.amazonaws.com/topo-slovenia/z15/11652/17705.png. From another map I got WGS84 (approximate) coordinates of upper left corner of this tile: 46° 2' 34,28" N, 14° 30' 47,06" E. With the help of proj4.js I converted this from EPSG:4326 to EPSG:3857 (or at least I hope I did) and got 5125467.235607444, 1633146.5805398274. I don't have the slightest idea how this relates to 11652, 17705.


texts: Looking for a shapefile of the US Southwest Border Sectors

texts:
I'm looking for a shapefile similar to the following image:

I'd just like the outlines of the sectors. I'm only concentrating on the southwest border which is San Diego all the way to Rio Grande Valley, but I'm ok with data on all of the sectors.

Does anyone know where I can find something like this?


texts: Trouble getting area of census block with OGR

texts:
I'm writing a small ogr/gdal script to get the area of census blocks:

in_file = data/zips/tabblock2010_01_pophu.zip
ds = ogr.Open(input_file)
lyr = ds.GetLayerByIndex(0)
lyr.ResetReading()
for feat in lyr:
  geom = feat.GetGeometryRef()
  print geom.GetArea()

This gives me numbers like 0.0003370560154999829 for block 010310106002015.  Is this an area in meaningful units?  If so, which ones?  Alternatively, is it something that doesn't make sense like area in lat-long space, where I need to do some kind of projection before I can get the area in real distance units?

EDIT: If I call DumpReadable() I get:

POLYGON ((-86.121456 31.303452,-86.121207 31.304302,-86.120986 31.305055, ...

so it does look like this is in lat/long units, and so probably gives me area in lat-long space.  How do I get real area?


texts: Click and Double click, deselect interaction OpenLayers 3

texts:
I have tow ol.interaction.Select, one to the single click and the other one to double click.

var selectDoubleClick = new ol.interaction.Select({
    multi: true,
    condition: ol.events.condition.doubleClick,
    style: function(feature) {
        //stuff here;
    }
    return feature;
}

selectDoubleClick.on('select',function(event) {
    //stuff here
}

var selectSingleClick = new ol.interaction.Select({
    multi: true,
    condition: ol.events.condition.singleClick,
    style: function(feature) {
        //stuff here
    }
    return feature;
}

selectSingleClick.on('select',function(event) {
    //stuff here
}

But when i want to deselect the double click, i have to click twice in the map. Is any way to deselect in double click with just one click?


texts:
I've been using GIS for a dissertation project and am having a little trouble with using a hillshade on a specific DEM.

My project is looking at Cyprus, and as such I've looked around and found a really helpful Mosaic to Raster DEM file that looks to be of high quality.

When I run a hillshade on it however, it looks awful and doesn't really resemble the original DEM. I know it's not an issue with this DEM, and definitely an issue with me as I have the same problem with Bathymetry data too.

I've attached a couple of images of what I mean to hopefully better explain, but the gist of it is the hillshade looks awful when coloured, and nothing like it should. Similarly, this is the case for Bathymetry data, even when I alter the normal 345 angle to account for reversal.

Wondering if anyone could help, would be massive boon to this project getting a few really high quality maps in there.

Above is an album to my issue. Maybe it's a colour problem? But no matter what I try, it looks ugly.

Big props to the creator of the mosaic to raster - Paraskeva Charalambos, University of Edinburgh


texts: Writing numpy arrays to irregularly shaped, multiband raster

texts:
I have the following numpy arrays:

predictions; type is float32; shape: (94, 1461)

lats; float32; shape: (94,); centroid of cell in vector layer, decimal degs.

lons; float32; shape: (94,); centroid of cell in vector layer, decimal degs.

94 corresponds to the number of cells in the vector grid below
1461 corresponds to the number of bands

I've looked at similar questions but most of them discuss exporting full (i.e. no missing, rectangular) rasters or single-band rasters. I was wondering how could I export a multi-band raster with multiple noData cells similar to the gridded vector layer here:

Additional info:
spatial ref: epsg:4326, GTiff, cell size 0.0270

My attempt:

xmin,ymin,xmax,ymax = [lons.min(),lats.min(),lons.max(),lats.max()]
xres = 0.0270
yres = 0.0270
geotransform=(xmin,xres,0,ymax,0, -yres)
output_raster = gdal.GetDriverByName('GTiff').Create(output_file, 9, 17, 1461, gdal.GDT_Float32)  # Open the file
output_raster.SetGeoTransform(geotransform)
for i in range(1461):
        output_raster.GetRasterBand(i+1).WriteArray(predictions)
        output_raster.GetRasterBand(i+1).SetNoDataValue(noData)
output_raster.FlushCache()

Throws a: line 229, output_raster.GetRasterBand(i+1).WriteArray(predictions)
ValueError: array larger than output file, or offset off edge


texts: Best practices when employing file geodatabase that is going to be used and edited by multiple users?

texts:
I have a file geodatabase that contains 3 feature classes in it and each feature class contains a lot of data that will be viewed and edited daily by 6 other users.  I am looking for any advice from people who have been in the same situation and the best practices to avoid edit lock outs and to keep everything running smoothly as possible. Is these any open source tools I can use that will enable multi user editing on a file geodatabase? Thank you!

My work only has ArcView level licencing with a couple extensions (Business Analyst, Network Analyst and Spatial Analyst)


texts: Change default "My template" directory for Arcmap layout view?

texts: So the default directory for "My templates" for switching the layouts in layout view is C:\Users\xxxx\appdata\raoming\esri\desktop..... so premade layouts can be accessed quickly. I need these same templates across a couple computers with access to the same network. Is it possible to point this default template folder to another drive, so that they all share the same templates?


texts: Generating flow map with varying width and alpha level which accumulates using ArcMap

texts:
I'm using ArcMap.
I have a very large line shapefile with ~300'000 links which all have the attribute flow (a numeric variable from 0-200).
My goal is to display a map like in this example (produced in QGIS):

This map is from the following this link. There are more
Basically I want to have a varying (and adjustable) width of the links, according to the value of flow. Ideally, the color of the links should contain an alpha value, to see overlapping links accumulate.
I have not been able to find anything on codeshare that could help me with this issue.


texts: gIntersection error: linearRing not created

texts:
I am trying to do something relatively simple but with large shapefiles. I have a shapefile of protected areas in the USA:

and watershed boundaries:

I want to calculate the percentage and proportion of protected areas in each watershed. I want to retain the designation of each protected area as well (i.e., federal, state, etc.)

To do this I started to try to intersect the two shapefiles:

  pas <- readOGR(dsn = ".", layer = "PADUS1_4Combined") # protected areas

  huc <- readOGR(dsn = ".", layer = "HUC2_USA_PRJ") # hydrologic units (2 digit, USGS)

Checked the CRS

  crs(huc)

  crs(pas)

  crs(huc)
  CRS arguments:
   +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0

  crs(pas)
  CRS arguments:
   +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0

Intersection:

  test <- gIntersection(pas, huc)

I get this error:

  Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, drop_lower_td, unaryUnion_if_byid_false,  :
    rgeos_crdMat2LinearRing: linearRing not created

I searched for this error and cannot find anything. I am dealing with a very large dataset here so I am wondering if this is the issue.

Is there a workaround using simple features (sf)?


texts: R: How can I create a polygon from several polygons and a line?

texts:
I have some spatial objects that I want to combine, namely three polygons and a line. The ultimate goal is to create a polygon that I can overlay on an area of a map of administrative divisions of the UK to highlight a certain zone.

If it were a case of combining four polygons, I would have a number of options, but as one of the objects is a line I'm at a bit of a loss as to how to tackle it. I would have thought combining polygons and lines would be a common operation but it seems not.

Here's the map produced by my code (see below).

The red line starts and ends on the coastline or, strictly speaking, on the bank of a river leading to an estuary at the south end, as can be seen in the following images.

So, because there is an intersection I assume it is possible to group the elements together, I just don't know now. To summarise...

Code follows:

library(httr)
library(ggplot2)
library(scales)

response <- GET(url="https://dl.dropbox.com/u/4115850/spobjects.rda")
load(rawConnection(response$content))

class(dist.shape) # map with administrative divisions
class(terr.shape) # line defining the landward side of the zone
head(dist.df) # data frames created using fortify() function
head(terr.df) # ditto

ggplot(data = dist.df, aes(x = long, y = lat, group = group)) +
     geom_polygon(colour = "grey50", fill = "grey80", size = 0.5) +
     geom_path(data = terr.df, aes(x = long, y = lat, group = group),
               colour = "red", size = 0.5) +
     theme()


texts: ArcGIS Server security and Windows domains - issues for non-IE browsers

texts:
I'm looking for some advice / experience in ArcGIS Server security with a Windows domain. Normally in this situation I would do some prototyping but I don't have access to a domain.

I'm planning to create ArcGIS map services that are restricted by roles stored in a domain.  These map services will be used by a map page in a custom .NET application. Both the application and ArcGIS Server will run on the same domain and potentially in the same IIS site (e.g. Default Web Site).

The .NET application will immediately require Windows authentication, so I know that IE users will pass through seamlessly but other browsers will prompt for a login. It's not clear to me if non-IE users will also be prompted for a login when first accessing a restricted map service.

Q1: In the above configuration (same domain, same IIS site) would non-IE users get prompted a second time to access a restricted map service (e.g. http://domain_name/ArcGIS/.../service_name), or would the credentials they provided for http://domain_name/application/ be used for ArcGIS?
Q2: Would this change if a different IIS site (but still the same domain) was used for ArcGIS Server and the .NET application?

Another (perhaps simpler) option I'm considering is having my map page's server-side code generate a security token on load by making a request to the ArcGIS token service with expiry time calculated as # hours until 17:30. This token would be passed to the client in JavaScript. In this approach only a single account needs access to the map service. The security token contains

  the user's name, expiration time and other information (source)

Q3 does this 'other information' include anything that restricts the token to the IP or user-agent of the requestor - anything that will stop my server creating it on behalf of my client?

EDIT from further research it seems that a token can be explicitly bound to a given IP (provided when creating it) or to a referrer. In this case Q3 is answered 'no'

Any thoughts on these three questions would be super welcome. Thanks


texts: Combining data from overlapping polygons and assign

texts:
I have two multipolygons layers:

Layer A: Census blocks with data from 2010 (last census).

Layer B: A study about poor neighborhoods with data from 2016 (POOR NH).

Both layers have a "population" attribute.

Aim: to update the CENSUS data with the POOR NH information (more updated and accurate).

I want to merge and replace the data from layer B to A, but assigning the "population" value proportionally (and only in cases where census blocks contains "poor places").

As you will see, being two different data matrices (pointing to different objectives), the sizes of the polygons do not match (the POOR NH are generally larger). Therefore, what I want is to assign the values from Layer B to Layer A, but respecting the proportions in which they overlap, simultaneously. That is, if a 30% (example: 1000) of Layer B is within 50% (2000) of Layer A, that '300'  replaces the '1000' from the other layer. The new "population" value of layer A, instead of being '2000', must be '1300'. Of course, assuming a homogeneous distribution of the population throughout the census blocks.

I've tried with intersection, split with lines (converting the polygon Census layer to polylines) and spatial join, but all generates multiples duplicated rows, because each combination creates a new polygon with all the features from the previous ones.

I assume that what I want to achieve requires several steps, but I'm stuck.

I've already checked:

How to combine data from overlapping polygons?