Sunday, 31 July 2016

Retrieving start and end point coordinates with ArcPy?



How do I retrieve start and end point coordinates with ArcPy for a polyline feature class?


I expect to pass a segment identifier to a subroutine and have it pass back start and end coordinates. The Field Calculator method doesn't work for me, because I need the value for other calculations that can't be performed within it. (I also prefer not to change the data to store these coordinates as attributes.) I am attempting to calculate address ranging for a "center out" addressing scheme. The address value depends on the distance to the "county center".




javascript - Leaflet load JSON data


Does my code have problem? Because if I run in Visual Studio, it doesn't show error, but the file didn't load on my map. I put the JSON file and HTML on my desktop.


This is my code :


    var osmUrl = 'https://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png',
osmAttrib = '© OpenStreetMap contributors',
osm = L.tileLayer(osmUrl, { maxZoom: 20, attribution: osmAttrib });
var map = L.map('map').setView([24.151687694799833, 120.64116954803465], 15).addLayer(osm);

L.control.scale().addTo(map);


$.getJSON("taiwan.json", function (data) {
var datalayer = L.geoJson(data, {
onEachFeature: function (feature, featureLayer) {
featureLayer.bindPopup(feature.properties.date);
}
}).addTo(map);
map.fitBounds(datalayer.getBounds());
});




This is my json (part of one)


{"type": "FeatureCollection","features": [
{
"type": "Feature",
"properties": {
"date": "1\/5\/1995 6:14:00",
"id": "第001號",
"lat": "24.96",

"lng": "121.7",
"depth": "91.5",
"magnitude": "4.5",
"desc": "台北市地震站東偏南方20.2公里",
"_id": 49710583
},
"geometry": {
"type": "Point",
"coordinates": [ 121.7, 24.96 ]
}

},
{
"type": "Feature",
"properties": {
"date": "1\/10\/1995 15:55:00",
"id": "第002號",
"lat": "23.68",
"lng": "121.43",
"depth": "3.8",
"magnitude": "5.1",

"desc": "花蓮西林地震站南方14.6公里",
"_id": 49710584
},
"geometry": {
"type": "Point",
"coordinates": [ 121.43, 23.68 ]
}
},

Answer



You need to define styles or pointToLayer function in L.geojson to create markers



var geojsonMarkerOptions = {
radius: 8,
fillColor: "#ff7800",
color: "#000",
weight: 1,
opacity: 1,
fillOpacity: 0.8
};
$.getJSON("taiwan.json", function (data) {
var datalayer = L.geoJson(data, {

onEachFeature: function (feature, featureLayer) {
featureLayer.bindPopup(feature.properties.date);
},
pointToLayer: function (feature, latlng) {
return L.circleMarker(latlng, geojsonMarkerOptions);
}
}).addTo(map);
map.fitBounds(datalayer.getBounds());
});

accuracy - If the GPS navigation message takes 12 1/2 minutes to cycle, how can receivers update your position every second?


According to Wikipedia's article on GPS,



Each GPS satellite continuously broadcasts a navigation message... Each complete message takes 750 seconds (12 1/2 minutes) to complete.




My understanding of how GPS works is that the receiver receives the broadcast message and calculates location based on that. So how can a GPS system start delivering reasonably accurate locations within a few seconds?


Wikipedia mentions that the message is divided into frames and sub-frames:



The message structure has a basic format of a 1500-bit-long frame made up of five subframes, each subframe being 300 bits (6 seconds) long... [A] complete data message requires the transmission of 25 full frames.... [T]his gives 750 seconds to transmit an entire almanac message (GPS). Each 30-second frame begins precisely on the minute or half-minute as indicated by the atomic clock on each satellite.



But the iPhone's GPS, for example, delivers about a location point a second. Even if a frame or a subframe is enough, that's 6-30 seconds. How can the hardware update its position with reasonable (reported) accuracy every single second? Is the iPhone lying to me?




qgis - Removing 'tails' from line features?



I have a set of polygon centerlines created using voronoi polygons. The method seems to work well for what I want to accomplish, but these 'tail' remnants are throwing off successive analysis (attached an image). Is there a way to remove them from the lines automatically? Note: I would prefer if this could be done with QGIS tools, but can create a python script if necessary.


I would be open to switching methods for obtaining polygon centerlines. I know there is a 'Polygon Centerline' tool in ArcMap, but I don't have access to it.





wmts layer in Openlayers from Geoserver


I have some layers in Geoserver(latest stable) that came from PostGIS data (they belong to a store with PostGIS type and gets data from a PostGIS database)


In those layers there are only WMS settings, not WMTS



I went to Caching Defaults in Geoserver, the WMTS Service is already enabled.


So I went to my Openlayers and try this, so I can get a PostGIS-based layer, from Geoserver to the map.


var projection = ol.proj.get('EPSG:3857');
var textent = ol.proj.transformExtent([2297128.5, 4618333, 2459120.25, 4763120], 'EPSG:900913', 'EPSG:3857');
var projectionExtent = projection.getExtent();
var size = ol.extent.getWidth(projectionExtent) / 256;
var resolutions = new Array(14);
var matrixIds = new Array(14);
for (var z = 0; z < 14; ++z) {
// generate resolutions and matrixIds arrays for this WMTS

resolutions[z] = size / Math.pow(2, z);
matrixIds[z] = z;
}

var ait = new ol.layer.Tile({
opacity: 0.7,
extent: textent,
source: new ol.source.WMTS({
url: 'http://localhost:8080/geoserver/mymap/wmts?',
layer: 'mymap:planet_osm_polygon, mymap:planet_osm_line, mymap:planet_osm_roads, mymap:planet_osm_point',

matrixSet: 'EPSG:3857',
format: 'image/png',
projection: projection,
tileGrid: new ol.tilegrid.WMTS({
origin: ol.extent.getTopLeft(projectionExtent),
resolutions: resolutions,
matrixIds: matrixIds
}),
style: 'default'
})

})

planet_osm_polygon, planet_osm_line, planet_osm_roads and planet_osm_point are composing the final layer. In simple WMS settings would be params: {'LAYERS': 'mymap:planet_osm_polygon, mymap:planet_osm_line, mymap:planet_osm_roads, mymap:planet_osm_point'.... but, in the WMTS case I dont know how I should set them, so that why I set them in layer


This does not work, I get


GET http://localhost:8080/geoserver/mymap/wmts?layer=mymap%3Aplanet_osm_polygon%2C%20mymap%3Aplanet_osm_line%2C%20mymap%3Aplanet_osm_roads%2C%20mymap%3Aplanet_osm_point&style=default&tilematrixset=EPSG%3A3857&Service=WMTS&Request=GetTile&Version=1.0.0&Format=image%2Fpng&TileMatrix=9&TileCol=286&TileRow=196 404 (Not Found)


Am I using the wrong settings or Geoserver does not support WMTS for PostGIS-based layers?




Saturday, 30 July 2016

arcgis desktop - Calculate frequency of a attribute field per polygon



I have spatially joined a polygon file with census tracts of areas in Houston with a theft occurrence point shapefile. The attribute table contains a field that represents different classes of theft (ie. Class A: $0-$1000, Class B: $1001 -$2000 etc). I want to calculate the frequency of classes in each of my census tracts. So in essence how many different theft classes occur in each of the census tracts. I have figured out how to get the count, but this is not what I want. How would I go about obtaining and representing the frequency?




Getting user's QGIS version using PyQGIS?



I am trying to write a Qgis plugin, which should work on Qgis2 and Qgis3. How can I find out which Qgis version the user is using, so the right imports are made?



Answer



In QGis < 3, qgis.core has a QGis object:


>>> from qgis.core import QGis

with various formats of the version string:


>>> QGis.QGIS_VERSION
u'2.18.5'
>>> QGis.QGIS_VERSION_INT
21805


and the release name, fwiw:


>>> QGis.QGIS_RELEASE_NAME
u'Las Palmas'

In versions 3 and above, this module is renamed Qgis for consistency.


There's a bit of a bootstrap problem here, since you can't get the version to figure out if you need to import Qgis or QGis without knowing. So you probably have to wrap it in a try and catch the exception.


try:
from qgis.core import Qgis
except ImportError:

from qgis.core import QGis as Qgis

Restricting Tin interpolation in QGIS



I am trying to do a bathymetric map for a lake using QGIS. However, I see no option to impose a mask so interpolation does not extending beyond the boundaries of my lake. I have traced around the lake (polygon) to which I would like to restrict the interpolation. Is there a "set mask" function in QGIS that would restrict all subsequent calaculations? I know how to convert the polgon to raster if needed. I have cliped to the lake boundaries but it still runs over. Do I need to go to GRASS to do this? Using version 1.8 on Windows 7.


Thanks


Alyre




Is ServerVector broken in OpenLayers 3.5?


Is ServerVector broken in OpenLayers 3.5?


Just running the following gives an error:


var wfsTest = new ol.source.ServerVector({
format: new ol.format.GeoJSON()
});

The error is: "Object doesn't support this action"


Am I doing something wrong?? It used to work fine in previous OpenLayers 3.x versions.



Answer




According to the change log of 3.5.0:



The ol.source.ServerVector class has been removed. If you used it, for example as follows:


var source = new ol.source.ServerVector({
format: new ol.format.GeoJSON(),
loader: function(extent, resolution, projection) {
var url = …;
$.ajax(url).then(function(response) {
source.addFeatures(source.readFeatures(response));
});

},
strategy: ol.loadingstrategy.bbox,
projection: 'EPSG:3857'
});

you will need to change your code to:


var source = new ol.source.Vector({
loader: function(extent, resolution, projection) {
var url = …;
$.ajax(url).then(function(response) {

var format = new ol.format.GeoJSON();
var features = format.readFeatures(response,
{featureProjection: projection});
source.addFeatures(features);
});
},
strategy: ol.loadingstrategy.bbox
});

Exposing GRASS GIS add-on in QGIS Processing framework?


Having QGIS and GRASS GIS 7 installed, how can I expose a GRASS GIS add-on, which has been installed using g.extension, in QGIS' Processing framework?




arcgis online - Creating simple maps to view on Mobile devices?



My organization wants to create a few maps that will be interactive and viewable via phones (android, iphones,etc) and ipads.


We already have a solution using Geocortex (Silverlight) for browser which cannot display on most devices. I've been looking into HTML5, Javascript and ArcGIS Online.


What I need is something that displays nicely for mobile use and that is easy to set up. Looking at Javascript is slightly intimidating to me and ArcGIS Online maps don't display well in the phone's browser.


I do not come from a Comp Sci or development background so if there's already an existing interface that makes it easily to display for mobile I'd prefer that rather than messing with code. The map would contain point features with hopefully a pop-up containing a picture and a short description of the property.


I'm pretty computer savvy but I'd like this to be as painless as possible.


Are there any other interfaces that exist that make this easier?


The mobile map will be only for viewing purposing, no editing or capture.



I'm open to looking at other platforms outside of ArcGIS as well.


As an added bonus (but not mandatory), I've been looking into Story Maps and plan on using one in the future, specifically the Short List. However, I don't see this working out for a mobile device as its very picture heavy and would crowd the display. If anyone has any other alternative ideas I'm open to them!


I've been looking at Leaflet and Mapbox as a possible solution for mobile use - there is plenty of documentation - but still a lot of code!


For now we've gained access to HTML5 viewer for geocortex and will be using it for our mobile sites. If there are any other alternatives I haven't already listed above - that are about on-par with ArcGIS Online for ease of use for non-programmers, I would be interested to hear about them.



Answer



You're right to look into Mapbox, despite not having a Computer Science background. Mapbox has a utility called TileMill which is essentially a web map development studio that has a familiar GIS-like interface. It is based in Javascript, but is designed to make development easy for non-developers. It gives you all of the code as well, though you may or may not actively use it.


https://www.mapbox.com/tilemill/


Friday, 29 July 2016

raster - Cannot properly import MODIS MOD07_L2 .hdf files into R using rgdal


I've been trying to accurately import and process some .hdf files from the MODIS Atmospheric Profile product (MOD07_L2) in R for several days now. Still, there's something going wrong during data import. The below code can be reproduced using one exemplary .hdf file (MOD07_L2.A2013001.0835.006.2013001192145.hdf) that can be downloaded from Dropbox.


library(rgdal)

# Extraction of metadata via `GDALinfo`
filename <- "MOD07_L2.A2013001.0835.006.2013001192145.hdf"
gdalinfo <- GDALinfo(filename, returnScaleOffset = FALSE)

metadata <- attr(gdalinfo, "subdsmdata")

# Extraction of SDS string for parameter 'Skin_Temperature' (formerly 'Surface_Temperature')
sds <- metadata[grep("Skin_Temperature", metadata)[1]]
sds <- sapply(strsplit(sds, "="), "[[", 2)

# Raster import via `readGDAL`
sds.rg <- readGDAL(sds)

So far, so good, but here comes the confusing part:



> summary(sds.rg$band1)
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
-14870 -14850 -14850 -14840 -14840 -14820 53529

Considering the fact that Skin_Temperature has an official valid range from 150 to 350 K (see MOD07_L2:Format & Content), the mean would inherit a value of


> (-14840 - (-15000)) * 0.01
[1] 1.6

after considering the corresponding add_offset (-15000) and scale_factor (0.01). Note that we're still talking about Kelvin, not °C. Extracting SDS No. 8, i.e. Skin_Temperature, using


library(gdalUtils)

gdal_translate(filename, dst_dataset = "tmp.tif", sd_index = 8)

and opening the resulting file called "tmp.tif" in QGIS resulted in seemingly reliable values centered around 15000, i.e. roughly 27 °C. However, importing "tmp.tif" back into R using raster again resulted in values comparable to the ones shown above:


> summary(raster("tmp.tif"))
tmp
Min. -14867.31
1st Qu. -14848.13
Median -14845.89
3rd Qu. -14840.53
Max. -14819.93

NA's 0.00

I've been searching the internet and stumbled across similar problems related to rgdal. However, when I tried to cast toUnSigned on band 1 of my previously generated 'SpatialGridDataFrame', I received the following error message:


> toUnSigned(sds.rg$band1, 16)
Error in toUnSigned(sds.rg$band1, 16) : band not integer

Apparently, the data imported into R is not even of type integer (what it is supposed to be), but numeric:


> sds.rg$band1[1:5]
[1] NA NA -14839.40 -14840.25 -14839.26


Is there an apparent mistake in my code, or is there any point I miss when importing the .hdf and .tif files using rgdal? I would be extremely grateful for any kind of help.



Answer



The answer is surprisingly simple.


sgr_lst <- readGDAL(sds, as.is = TRUE)

solves the issue. The only thing left to do then is transform the resulting SpatialGridDataFrame to a proper Raster* object using raster(). For this purpose, it is necessary to retrieve the west, east, south, and north bounding coordinates (see ?extent: xmin, xmax, ymin, ymax) from the metadata e.g. via


meta <- attr(gdalinfo, "mdata")

## string patterns
crd_str <- paste0(c("WEST", "EAST", "SOUTH", "NORTH"), "BOUNDINGCOORDINATE")

## search for patterns in metadata
crd_id <- sapply(crd_str, function(i) grep(i, meta))
## extract information
crd <- meta[crd_id]
## create 'extent' object
crd <- as.numeric(sapply(strsplit(crd, "="), "[[", 2))
ext <- extent(crd)

The thus determined extent (together with the EPSG code) can then be passed on to the finally created RasterLayer which, after applying the scale factor and offset, definitely looks fine now.


## rasterize, apply offset and scale factor

rst_lst <- (raster(sgr_lst) - -1.5e+04) * 1.0e-02
## set extent and coordinate reference system
extent(rst_lst) <- crd
projection(rst_lst) <- "+init=epsg:4326"

lst


arcpy - What is the right SQL query for this case?




I want to use input string in part of the SQL query, and the code has been posted below.


import arcpy
import os

in_features = arcpy.GetParameterAsText(0)
target_county = arcpy.GetParameterAsText(1)
out_path = os.path.dirname(in_features)
out_name= "%s.shp" % target_county
where_clause=""""CNTY_NM" = %s""" % target_county

arcpy.FeatureClassToFeatureClass_conversion(in_features, out_path, out_name, where_clause)

When I run the script, it gives the error. enter image description here enter image description here


I don't know what is wrong with the code. Any suggestions?


I am using ArcGIS 10.3.1 Desktop Advanced License.



Answer



In a SQL where clause, single quotes ("FIELD" = 'provided_value') are needed around the match string in the where clause. You are missing the single quotes around the provided value. It should be like this:


where_clause="\"CNTY_NM\" = '%s'" % target_county

arcgis desktop - Why is the result of a raster calculation integer when it should be decimal?


I'm trying to calculate the NDVI using raster calculator, but the result is integer instead of decimal.


NDVI = (NIR - RED)/(NIR + RED)

The result is normally in the range -1 to 1, but when i tried to calculate de result is -1, 0 and 1, and it does not have have any number in between like, -0.9 or 0.1. And i need those values so i can classify the result in term of soil occupation





One-to-many join in QGIS



I am looking for a possibility to make a one-to-many table-join in QGIS.


I have a shapefile with postal codes and would like to join them with survey results (more than one). In other words, there is a specific postal code (shapefile) on one side and there are many survey participants on the other side.


postal **code1** ----- survey participant **1**
postal **code1** ----- survey participant **2**

postal **code1** ----- survey participant **3**
postal **code2** ----- survey participant **1**
postal **code2** ----- survey participant **2**
and so on ...

If I do a normal join (right click on Layer --> Properties --> Join) , I just get the first entry of the table with the postal codes.


I found out a way to do it ArcGIS with the tool "Make Query Table", but I´m looking for a possibility to do it in QGIS.




Create common vertexes between shared polygon boundaries in PostGis?


Given two polygon tables L1 (blue) and L2 (red) that share a common boundary but different amount of vertexes, how can one make the shared boundary identical (they should share the same vertexes) ?


EDIT:


By identical I mean that 'If we say the Blue Layer is the reference layer, then the edge of the Red that is common with the blue should afterwards have 3 vertexes instead of two that is now'


enter image description here



Answer



What you have in PostGIS for that is ST_Snap http://postgis.net/docs/ST_Snap.html


This is a four-vertex polygon


POLYGON (( 420 420, 420 520, 700 520, 700 420, 420 420 ))


This is an adjacent five-vertex polygon that has one vertex in the middle of the common boundary at location (582 420)


POLYGON (( 420 420, 582 420, 700 420, 700 320, 420 320, 420 420 ))

With the following SQL you can make the four-vertex polygon to snap to the five-vertex polygon:


select ST_AsText(ST_Snap(ST_GeomFromText('POLYGON (( 420 420, 420 520, 700 520, 700 420, 420 420 ))'),ST_GeomFromText('POLYGON (( 420 420, 582 420, 700 420, 700 320, 420 320, 420 420 ))'),1))

The result which now contains a new vertex at (582 420)


POLYGON((420 420,420 520,700 520,700 420,582 420,420 420))


Looks perfect, doesn't it? Unfortunately it is not so simple with real world data. For the first you must first snap A to B and then B to A for getting the same vertices to both polygons. And if the polygon boundaries make a curve then there will be small overlaps and gaps between the polygons which are created from different vertices and the result of using ST_Snap function may not be satisfying.


For cleaning bigger datasets I would use OpenJUMP which has quite a clever tool just for this purpose. It is processing the whole layer and it takes care of adding vertices to both parties at one run. With very complicated polygons the result may still not be perfect but from my experience the result is better than by using ST_Snap.


enter image description here


enter image description here


EDIT


An advanced example about snapping A to B and then B to snapped-A


Source geometries


POLYGON (( 280 580, 660 580, 660 460, 460 480, 280 460, 280 580 ))
POLYGON (( 280 460, 520 480, 660 460, 660 380, 280 380, 280 460 ))


enter image description here


SQL to use


SELECT ST_AsText(ST_Snap(ST_GeomFromText('POLYGON (( 280 580, 660 580, 660 460, 460 480, 280 460, 280 580 ))'),ST_GeomFromText('POLYGON (( 280 460, 520 480, 660 460, 660 380, 280 380, 280 460 ))'),20));

SELECT ST_AsText(ST_Snap(ST_GeomFromText('POLYGON (( 280 460, 520 480, 660 460, 660 380, 280 380, 280 460 ))'),
ST_Snap(ST_GeomFromText('POLYGON (( 280 580, 660 580, 660 460, 460 480, 280 460, 280 580 ))'),ST_GeomFromText('POLYGON (( 280 460, 520 480, 660 460, 660 380, 280 380, 280 460 ))'),20),20));

The result:


enter image description here


Thursday, 28 July 2016

pyqgis - S57/ENC QgsVectorLayer with stand-alone python script


I'm running the OSGeo4W64 Qgis 3.2.0.


In the gui I can add a vector-layer using the same data as in the example below without any problems. Adding the layer using the python console works too.


When trying to do the same from a stand-alone script I get "ERROR 4: S57 Driver doesn't support update."


Adding other kml files works, but not S57/ENC


Looking in the code where the error is written, it looks like the difference may be the file mode used. I don't want or need to edit, so read-only/update=False is fine.


import os

from qgis.core import *

regi = QgsProviderRegistry.instance('C:\\OSGeo4W64\\apps\\qgis\\plugins')

QgsApplication.setPrefixPath(r'C:\OSGeo4W64\apps\qgis')
app = QgsApplication([],False)
app.initQgis()

l1=QgsVectorLayer(r"D:\ENC_NOAA\All_ENCs\ENC_ROOT\US1AK90M\US1AK90M.000|layername=M_COVR","COVR1","ogr")
print("Valid" if l1.isValid() else "Not Valid")


I have set CPL_DEBUG=ON and the output shows



ERROR 4: S57 Driver doesn't support update.


GDAL: GDALOpen(D:\ENC_NOAA\All_ENCs\ENC_ROOT\US1AK90M\US1AK90M.000, this=0000018CCEB33550) succeeds as S57.


GDAL: GDALClose(D:\ENC_NOAA\All_ENCs\ENC_ROOT\US1AK90M\US1AK90M.000, this=0000018CCEB33550)


Not Valid



Can I control the open mode via the uri?


Why the difference between interactive and stand-alone?





field calculator - Does QGIS support parsing of date strings in tables?


It is straightforward to load GPX files (GPS downloaded tracks / points) into QGIS, then save and manipulate them as shapefiles. However, the recording date (in the comment (cmt) field) of waypoints is in the format YY-MMM-DD HH:MM:SS (e.g. 28-AUG-11 11:57:51) and it would be possible to sort by date only if the format could be parsed to something closer to the ISO 8601 standard (YYYY-MM-DD ...).


String manipulation in the field calculator can take you partway: '20' || substr(cmt,8,2) || '-' || substr(cmt, 4, 3) || '-' || substr(cmt, 1, 2) || ' ' || substr(cmt, 11, 8) However, you're then left with YYY-MMM-DD HH:MM:SS - i.e. AUG instead of 08.


Is there a straightforward way to parse the date string directly, providing integer months? (Ideally, something like the R strptime() function)



Answer



If you're willing to use PostGIS, you can convert the GPX date field to a Postgres timestamp field like so (this works as long as the months are always abbreviated to 3 letters, which I believe is the case):



to_timestamp('28-AUG-11 11:51:57', 'DD-MON-YY HH:MI:SS');

Or, in the field calculator, you can do something like this:


replace(replace(replace(substr(date, 4, 3), 'JAN', '01'), 'FEB', '02'), 'MAR', '03')

(of course, put in more nested replaces for all the months). It isn't pretty but it works.


labeling - Colliding labels for point features in QGIS


I have a layer (point features) containing various points of interest (POIs), some of them have the exact same location. I would like to map these POIs and add a label showing their names. However, for those POIs that share the same location, QGIS just display the overlapping labels and it is not readable.


Is there a way in QGIS to automatically place labels for point features so that they never overlap, even though they belong to the same layer?


I found this relevant thread, but no solution to the problem.





arcgis desktop - can I extract overlap of polygons within the same layer?


I am trying to extract the overlap of polygons that are within the same layer in 10. The layer came from a buffered street centerline and I'm trying to extract and delete the intersections from the layer. I have tried 'intersect' but it's virtually impossible to tell which records are an overlap and which ones are not. We are trying Shape tools but so far the computer gave an 'overflow' error.




gdal - How can I get the proj4 string or EPSG code from a shapefile .prj file?



I have to work with a shapefile defined in an unusual projection, which my GIS software does not know about. How can I get the proj4 definition or the EPSG code for the projection?




Answer



The shapefile should have a .prj file which defines the projection. You can use it together with one of the following 3 options to get either the proj4 string, WKT definition or EPSG code.


To get proj4 definition:


If you have gdal installed on your system, you can use the gdalsrsinfo command line application to get the proj4 definition as the OGC WKT definition:


$ gdalsrsinfo .prj

Exampe:


$ gdalsrsinfo 7490.prj 

PROJ.4 : '+proj=aea +lat_1=-18 +lat_2=-32 +lat_0=0 +lon_0=24 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs '


OGC WKT :
PROJCS["Albers_Equal_Area_Conic_South_Africa",
GEOGCS["GCS_WGS_1984",
DATUM["WGS_1984",
SPHEROID["WGS_84",6378137,298.257223563]],
PRIMEM["Greenwich",0],
UNIT["Degree",0.017453292519943295]],
PROJECTION["Albers_Conic_Equal_Area"],
PARAMETER["False_Easting",0],

PARAMETER["False_Northing",0],
PARAMETER["longitude_of_center",24],
PARAMETER["Standard_Parallel_1",-18],
PARAMETER["Standard_Parallel_2",-32],
PARAMETER["latitude_of_center",0],
UNIT["Meter",1]]

To get EPSG number:


You can use the prj2epsg.org service to upload the .prj file and get a list of matching EPSG codes. Warning: prj2epsg.org returns the list of codes which match best. It may be that none of them are correct.


Using python:



If you are comfortable with Python, you can use the very useful script at Shapefile PRJ to PostGIS SRID lookup table? to get the WKT, proj4 string and EPSG code.


See also:


You can find definitions for most standard projections at spatialreference.org.


Use Orfeo Toolbox in R


I was wondering if it is possible to use the algorithm of Orfeo Toolbox in RStudio. I have been looking for some information but didn't find anything relevant. From the Orfeo toolbox I see that it can be access through command line so that may be an option?


I am working with Ubuntu Xenial 16.04



Answer



OTB algorithms are in OTB_installation_path/bin. So, as you can use them in command line (see more in documentation), the same way could be used trough R by system function.


You can create your owns functions to work with OTB in R. An example for Mean Shift Image Segmentation:



meanshift.segm <- function(otb.path = "", raster.in = "", out.path = "", name ="", filter.meanshift.spatialr = "5",
filter.meanshift.ranger = "15", filter.meanshift.thres = "0.1",
filter.meanshift.maxiter = "100", filter.meanshift.minsize = "100",
mode.vector.outmode = "ovw", mode.vector.inmask = "", mode.vector.neighbor = "false",
mode.vector.stitch = "true", mode.vector.minsize = 1, mode.vector.simplify = 0.1,
mode.vector.layername = "layer", mode.vector.fieldname = "DN", mode.vector.tilesize = 1024,
mode.vector.startlabel = 1){
# Set configuration
conf <- paste("-in",raster.in,"-filter meanshift","-filter.meanshift.spatialr",filter.meanshift.spatialr,
"-filter.meanshift.ranger",filter.meanshift.ranger,"-filter.meanshift.thres",filter.meanshift.thres,

"-filter.meanshift.maxiter",filter.meanshift.maxiter,"-filter.meanshift.minsize",filter.meanshift.minsize,
"-mode vector","-mode.vector.out",paste(out.path,"/",name,".shp",sep=""),"-mode.vector.outmode",mode.vector.outmode,
ifelse(missingArg(mode.vector.inmask),"",paste("-mode.vector.inmask",mode.vector.inmask)),
"-mode.vector.neighbor", mode.vector.neighbor,
"-mode.vector.stitch",mode.vector.stitch,
"-mode.vector.minsize",mode.vector.minsize,
"-mode.vector.simplify",mode.vector.simplify,
"-mode.vector.layername",mode.vector.layername,
"-mode.vector.fieldname",mode.vector.fieldname,
"-mode.vector.tilesize",mode.vector.tilesize,

"-mode.vector.startlabel",mode.vector.startlabel)
# apply function in command line
system(paste(otb.path,"/otbcli_Segmentation"," ",conf,sep=""))
# save configuration for futher use
write.table(x = conf,file = paste(out.path,"/",name,"_conf.txt",sep=""),row.names = F, col.names = F)
}

# usage, you can set any option listed above
meanshift.segm(otb.path = "OTB-5.8.0-Darwin64/bin", raster.in = "path to/rater_in.tif", out.path="/out/path", name = "test")


# or, to read the output into R
out_path = '/out/path'
lyr_name = 'test'

# usage
meanshift.segm(otb.path = "OTB-5.8.0-Darwin64/bin", raster.in = "path to/rater_in.tif", out.path=out_path, name = lyr_name)

shp <- readOGR(dsn=out_path, layer = lyr_name, driver = 'ESRI Shapefile')

Using QGIS libraries to write standalone program in C#?


I have investigate some liberaries which can help me to write stand alone applications.


Is there any way to use qgis liberaries in an standalone application which is written by C#?



Answer




Not really no. You only options for QGIS based applications are Python and C++, unless you want to build a web app which calls QGIS server to show the map then you can use what ever you want as it's just making WMS calls to the server.


Creating a standalone application with Python using QGIS is pretty easy.


Getting started with a standalone app in Python:



If you are using QGIS 2.4 there is some handy application bootstrap logic:


from qgis.core.contextmanagers import qgisapp

def main(app):
pass


if __name__ == "__main__":
with qgisapp() as app:
main(app)

That will take care of the bootstrapping logic normally needed to load QGIS into the right state.


(Sorry about all the plugs to my stuff, others feel free to add others)


python - qgis.core import error


I would want to use the Qgis api for custom desktop application. When I imported the qgis.core module, I get the following error


>>> import qgis.core 
Traceback (most recent call last):
File "", line 1, in
ImportError: DLL load failed: The
specified module could not be found.

A similar error was reported in the forums,im not sure if the exact solution was reached.


Python :ImportError: DLL load failed: The specified module could not be found



Forums do suggest that I would have to use the OSGeo Version of python. How would i get to choose the version of python to use?


I think this could be a path or a pythonpath.Have given my path and python path information below.Let me know if I had missed any of the directories that should be added for qgis core to get its dependencies met.


PATH=C:\Python26;C:\Python26\Scripts;C:\Python26\ArcGIS10.0;C:\Python26\ArcGIS10.0\Scripts;C:\OSGeo4W\apps\qgis-dev\bin


PYTHONPATH=C:\Python26;C:\OSGeo4W\apps\qgis\python;C:\OSGeo4W\apps\python27;C:\Python26\ArcGIS10.0;C:\Python26\Lib\site-packages



Answer



Try using the hints provided in command-line-environment-setup to set up your shell properly before launching your application.


Regards


Tim


qgis - Handle add new feature event and/or access feature before commit?


Is it possible with PyQGIS to catch "add new feature" event before commiting changes on the layer? or at least to access new features before commit.



Answer



The QgsVectorLayer has quite a number of signals that will may suite your needs. Documentation for 1.8 and master on these links.



Note there are minor changes to the API between 1.8 and master, and it appears that the python bindings have a few issues that are still unresolved.


This sample code defines a couple of listeners, then attaches them to a memory layer. If you start editing the layer you will see each signal triggers a line in the message log.


This code is (mostly) compatible with 1.8 and master and is intended to be pasted into the python console.


from PyQt4.QtGui import *

def logLayerModified( onlyGeometry = None ):
QgsMessageLog.logMessage( "layer modified" )
QApplication.beep()

def logFeatureAdded(fid):

QgsMessageLog.logMessage( "feature added, id = " + str(fid) )
QApplication.beep()

def logEditingStarted():
QgsMessageLog.logMessage( "editing started" )
QApplication.beep()

def logCommittedFeaturesAdded( layerId, addedFeatures ):
message = layerId + " has features added: "
for feature in addedFeatures:

message += str( feature.id() ) + ", "
QgsMessageLog.logMessage( message )
QApplication.beep()

layer = QgsVectorLayer( "Point", "Layer Signal Demo", "memory" )
QgsMapLayerRegistry.instance().addMapLayers( [layer] )

layer.layerModified.connect( logLayerModified )
layer.featureAdded.connect( logFeatureAdded )
layer.editingStarted.connect( logEditingStarted )

layer.committedFeaturesAdded.connect( logCommittedFeaturesAdded )

EDIT: Even with an incomplete set of signals you can probably get the access you require by performing you own analysis.


I have a naive implementation here that will slow down for large feature counts. Depending on your needs you may be able to optimize by reducing the amount of data being handled during the cross referencing. You may get speed gains by omitting the geometry or full attribute table from the query until you have identified the features you are interested in.


Note that this is designed for 1.8, feature iteration has changed significantly for master.


from PyQt4.QtGui import *

def layerModified( onlyGeometry = None ):
global establishedFeatureIds
newFeatureIds = []

layer.select()
feature = QgsFeature()
while layer.nextFeature( feature ):
id = feature.id()
if not id in establishedFeatureIds:
newFeatureIds.append( id )
message = "New features: "
for id in newFeatureIds:
message += str(id) + ", "
QgsMessageLog.logMessage( message )


def editingStarted():
global establishedFeatureIds
establishedFeatureIds = []
layer.select()
feature = QgsFeature()
while layer.nextFeature( feature ):
establishedFeatureIds.append( feature.id() )

layer = qgis.utils.iface.activeLayer()

if layer is None:
layer = QgsVectorLayer( "Point", "Cross-reference Demo", "memory" )
QgsMapLayerRegistry.instance().addMapLayers( [layer] )

layer.layerModified.connect( layerModified )
layer.editingStarted.connect( editingStarted )

Wednesday, 27 July 2016

arcgis server - Secure Print Service Failed to Execute



I have a web application that calls an ArcGIS 10.2.2 secure print service (accesses secured map services). I am calling the print service using a proxy page. Generally this works, but after roughly 30 minutes the print request fails for no apparent reason, and the ArcGIS Server logs contain the following:-


    Error executing tool.: Layer "internal": Unable to connect to map server at 
http://myserver:6080/arcgis/rest/services/Internal/Cadastral/MapServer.
Layer "jobs": Unable to connect to map server at
http://myserver:6080/arcgis/rest/services/Internal/General/MapServer.
Layer "wells": Unable to connect to map server at
http://myserver:6080/arcgis/rest/services/Internal/Oil_and_Gas/MapServer.
Failed to execute (ExportWebMap). Failed to execute (Export Web Map).

The error reported back from the http request is:




Error: Unable to complete operation



The rest of the application works perfectly. If I try printing later (say 20mins), or if I restart the Print Service Tool, it happily prints again for another ~30 mins. It is as if a connection is timing out. Token expiring? Has anyone had this problem before, or can someone suggest where I should direct my attention?




pyqgis - Making 'make'-command work for compiling QGIS-plugins?


How can I compile my plugin when the 'make'-command refuses to work in the OSGeo4W Shell?


Can I make the 'make'-command work?


Is there an alternative out there?


I'm trying to build a plugin for QGIS. It's plugin Plugin Builder creates a projectfolder. According to tutorials I have to open the OSGeo4W Shell and compile the resources of my plugin by navigating to the project folder and execute the make-command. However, this doesn't work for me, saying:



'make' is not recognized as an internal or external command,operable program or batch file.




I have a Windows sytem so that's not strange. Yet the OSGeo4W Shell should give me the abbility to use this command....I guesss. So I re-installed QGIS by using the Desktop Install by OSGeo4w in hope that it would configure the libraries correctly. Sadly with no effect.


Am I missing something?




openlayers 2 - How to set Zoom Level based on Overlays in Ordnance Survey Map?


I am using the Ordnance map as a backdrop in my website. when overlays the WMS from Geoserver, I am unable to increase the zoom level. Please see my code and kindly suggest it.


 var options = {

maxResolution: 0.55804443359375,
projection: "EPSG:27700",
units: 'm'
};
map = new OpenSpace.Map('map', options);

// Define or get the vector layer to be used
var vlayer = osMap.getVectorLayer();
var mlayer = osMap.getMarkerLayer();
vlayer.displayInLayerSwitcher = false;

mlayer.displayInLayerSwitcher = false;


// setup tiled layer
format = 'image/png';
var tiled = new OpenLayers.Layer.WMS(
"Oxford_Example",
"http://192.168.0.126:8080/geoserver/cite/wms",
{
transparent: 'true',

LAYERS: 'cite:Oxford_ExampleCapture',
STYLES: '',
format: format,
tiled: true,
tilesOrigin: osMap.maxExtent.left + ',' + osMap.maxExtent.bottom
},
{
buffer: 0,
displayOutsideMaxExtent: true,
singleTile: true,

ratio: 1,
isBaseLayer: false,
yx: { 'EPSG:27700': false }
}
);


osMap.addLayers([tiled]);



//Show Co-Ordinate on mouse over
osMap.addControl(new OpenLayers.Control.MousePosition());

//OverViewMap
var control = new OpenSpace.Control.OverviewMap();
osMap.addControl(control);
control.maximizeControl();

//ScaleBar
osMap.addControl(new OpenLayers.Control.ScaleLine({ displayClass: 'olControlScaleLine' }));

// osMap.addControl(new OpenLayers.Control.Scale({ displayClass: 'olControlScale' }));

//Layer Switcher
// var layerSwitch = new OpenLayers.Control.LayerSwitcher();
var layerSwitch = new OpenLayers.Control.LayerSwitcher({
div: OpenLayers.Util.getElement('layerswitcher')
});
osMap.addControl(layerSwitch)

But i can zoom in this level only. I am unable to the particular feature. See the below image:



enter image description here




Will ArcGIS for Desktop run on Windows 10?



Will ArcGIS 10 run on the new windows 10 technical preview?



I am setting up an interface on my Mac and there is the free option of downloading the Windows 10 technical preview, some sort of beta version, I'm assuming.


Just wanted to see if 10.1 will run ok on Windows 10 (ignoring the Mac aspect), as this is all I am using my windows interface for?




Tuesday, 26 July 2016

python - How to update only a single vector layer in QGIS map canvas?


I want to update a single vector layer as I am working on 8 vector layers and they are getting fast data from serial port so for updating I use mapCanvas().refresh() which refresh all the layer and I have almost 20 static S57 layer so the refresh takes a long time.




qgis - Load a specific spatial data from a postGIS table


A single postGIS table contains multiple type of spatial data(example: Linestring,Polygon,etc.,). How to load a specific spatial data from PostGIS table into QGIS by using pyqgis?




Sample code snippet is given below


uri = QgsDataSourceURI()    
uri.setConnection("xxx.xxx.xxx.xxx", "port", "dbname", "username", "password")
uri.setDataSource("public","tablename","geom")


uri.setWkbType(QGis.WKBPolygon) #geometry type has to be mentioned

uri.setSrid('4326')
vlay=QgsVectorLayer(uri.uri(),"test_poly","postgres")
QgsMapLayerRegistry.instance().addMapLayer(vlay)


How can I reproduce Getis-Ord GI* hot-spot analysis tool in QGIS?


I'm trying to find a hot-spot analysis plug-in for QGIS that would perform similarly to the Getis-Ord (Gi*) hot-spot analysis tool used in Arc, but have not been successful. Does anyone have a suggestion for where to find one?





landsat - Maximum ndvi for 3 months average in Google Earth Engine


I have a ndvi series which I am able to plot using below code:


var ndvi = l8.map(function(image) {
return image.select().addBands(image.normalizedDifference(['B5', 'B4']));

});
var ndviChart = ui.Chart.image.series(ndvi, point, ee.Reducer.mean(), 500);
print(ndviChart);

My goal is to calculate the average ndvi for every 3 months window and then calculate the maximum ndvi over that 3 months window. How to access data in ui.Chart.series and is there any way by which I can do the above calculation on Google server ?


In javascript it can be done using





I want to replace lista with ndviChart series. Can anybody help me in this direction?




Cube root on ArcGIS raster calculator?


How I could extract a cube root of the multiplication of 3 rasters in raster calculator?


Eg: cuberoot (ras1 * ras2 * ras3)




Python error can't find PyQt.Qtcore, when launching QGIS


I installed QGIS in Linux and it was working fine. I try to launch the program again and I get the following Python error:


Couldn't load plugin MetaSearch due to an error when calling its classFactory() method


Traceback (most recent call last):

File "/usr/lib/python2.7/dist-packages/qgis/utils.py", line 219, in startPlugin
plugins[packageName] = package.classFactory(iface)
File "/home/natalia/.qgis2/python/plugins/MetaSearch/__init__.py", line 29, in classFactory
from MetaSearch.plugin import MetaSearchPlugin
File "/usr/lib/python2.7/dist-packages/qgis/utils.py", line 478, in _import
mod = _builtin_import(name, globals, locals, fromlist, level)
File "/home/natalia/.qgis2/python/plugins/MetaSearch/plugin.py", line 28, in
from qgis.PyQt.QtCore import QCoreApplication
File "/usr/lib/python2.7/dist-packages/qgis/utils.py", line 478, in _import
mod = _builtin_import(name, globals, locals, fromlist, level)

ImportError: No module named PyQt.QtCore

However, I do have installed PyQt.QtCore. I think the problem is with my path and QGIS doesn't know how to look for the Python package.


I am really new to this, I am not sure how to fix it. The way that my path looks is:


PATH=/home/usr/.qgis2/python:$PATH
PATH=/home/usr/share/qgis/python:$PATH


kml - Visualization of pictures for points in QGIS


I'm working with QGIS Desktop 3.2.3.


A Colleague of mine worked on a Google Map and inserted in his Google Map Pictures. He exported the KML and gave it to me so that I can import it to QGIS.


The pictures are part of the attribute table as a URL. Is it possible to download the pictures into QGIS?




Monday, 25 July 2016

gdal - FWTools ogr2ogr in Python


I am trying to convert MapINFO tabs to ESRI Shapefiles using FW Tools. I downloaded the tools and can run the FW Tools Shell and convert no problem:


ogr2ogr -f "ESRI Shapefile" E:\CABWorking\MapINFO2SHP\SHP\LincsBoundary.shp E:\CABWorking\MapINFO2SHP\MapINFOData\Lincolnshire.tab

In the shell you can type python and use the shell as a python interpreter. Awesome, but I need this to work as I want to convert a ton of MapINFO files at once.


My problem is: I can't figure out how to get the following code to work in a python environment. I've been searching all morning and just have no clue. Any help would be awesome!


Thanks!



Answer



In Python the os module provides a quick and dirty way to run a command line argument such as the one you posted as part of a larger script as outlined in the Python documentation


import os

command = "ogr2ogr -f \"ESRI Shapefile\" E:\\CABWorking\\MapINFO2SHP\\SHP\\LincsBoundary.shp E:\\CABWorking\\MapINFO2SHP\\MapINFOData\\Lincolnshire.tab"
os.system(command)

Note the use of backslashes to escape characters such as " and \ in the command string.


As I said, this is quick and dirty, the Pythonic way to do this is to use the Subprocess library which provides the same functionality but in a more stable framework, and with some more features:


import os
command = "ogr2ogr -f \"ESRI Shapefile\" E:\\CABWorking\\MapINFO2SHP\\SHP\\LincsBoundary.shp E:\\CABWorking\\MapINFO2SHP\\MapINFOData\\Lincolnshire.tab"
process = subprocess.Popen(cmd, stdout=subprocess.PIPE, creationflags=0x08000000)
#the third command suppresses the opening of a new window for each process
#and can be removed if this is not the desired outcome

process.wait() #waits for one process to complete before launching another one



I used this answer to solve my issue. Thanks a ton @sgrieve. The final code ended up being:


import os, glob
inputDir = raw_input("==> ") # Had to have this code compatible with 2.3.4
outputDir = raw_input("==> ") # So that meant no easy file dialog boxes
TabFiles = glob.glob(str(inputDir) + '/*.tab')
for TabFile in TabFiles:
TabFileName = TabFile[(TabFile.rfind("\\"))+1:(TabFile.rfind("."))]

command = 'ogr2ogr -f "ESRI Shapefile" ' + outputDir + "/" + TabFileName + '.shp ' + TabFile
os.system(command)

I'd then open the FWTools Shell, enter python and call the script with:


execfile("script")

vector - Measuring the distance between lines and points in QGIS



Using QGIS - I have a vector of fault lines. I have point data in the form of a csv which I have also loaded into QGIS. I need to know how far each point is from the nearest fault line then enter this distance in a new column in the csv file. Any ideas?




data - Why can't I see all records in my personal GDB in QGIS?


I am running QGIS 1.8. I have a layer from an ESRI Personal GDB. There are 431 features in the layer. When you load the layer into QGIS, and look at the Metadata properties, it shows the number of features as 431. layer metadata


When I open the attribute table, the table only shows 180 features. Layer attribute table


This is a layer that had been loaded in straight from the Personal GDB. No filters or queries had been added to it, so all features should be showing.


I first came across this issue because another coworker was trying to access the same layer using QGIS 1.7. He was having the same issue.


A couple items of note:



  1. This layer was existing in the PGDB, and then using ArcGIS, I copied and replaced a new version from a File GDB into the PGDB.


  2. I am able to load other layers from the same PGDB into QGIS where all the features are displayed.

  3. The PGDB is a modified MS Access database. If you open the database with MS Access, then open the table for this layer, it contains all 431 records with attributes.


Here is a screenshot of the layer in QGIS: QGIS Screenshot


Here is a screenshot of the same layer in ArcGIS: ArcMap Screenshot


If anyone has some insight into why QGIS is able to see the actual number of features in the layer, but is only displaying a certain number, I would appreciate it.



Answer



Open the Personal GDB in MS Access to make sure the records not shown in QGIS all have valid data, including coordinates. I'd first check the 181st record (in your case), being the first record not shown in QGIS.


Given your confirmation in the comment above, you might file a bug with the QGIS team to request a visible error message in case invalid data is found.


arcgis desktop - How to display a .jpg file in ArcMap for georeferencing?



I'm struggling to georeference a .jpg map. When I load it in to ArcMap, I cannot see the picture at all, so I cannot georeference it. (I tried to zoom to the map layer but the map just does not appear).


Does anyone have a solution for this?



Answer



Instead of trying to zoom to layer you need to use the georeferencing toolbar to display your .jpg


Once you have the .jpg in your table of contents, navigate to where your image will be georeferenced (assuming you have a basemap) you can then open the drop-down and select fit to display:


Fit to display


This will center the image on your screen allowing you to start the georef process with it roughly in the right location.



clip - Create Polygon from holes in polygon ArcGIS


I have a polygon feature class with holes in it, which I'll refer to as Polygon 1:



enter image description here


I wish to create a polygon set of circles/clipped circles that are the holes as denoted 1->6 below. I believe I could accomplish this with the Clip tool and use of an overlaying polygon with the hole deleted, specifically this polygon, which I'll refer to as Polygon 2:


enter image description here


My understanding of the clip tool is that then I need to clip Polygon 2 by Polygon 1 leaving only the circles/clipped circles 1->6. As such I ran the Clip tool with the following inputs:


Input Features = Polygon2


Clip Features = Polygon1


I did not get the expected output, but rather a replica of Polygon1, see here:


enter image description here


I also tried reversing the inputs, the result is exactly the same.


What is the error in my understanding or my approach?




Answer



You should have used the ERASE (or difference) tool to obtain the holes:


inputFeature = Polygon 2
eraseFeature = Polygon 1

http://desktop.arcgis.com/en/arcmap/10.3/tools/coverage-toolbox/erase.htm


arcgis 10.0 - First attempt using ArcObjects: ESRI's Add-ins example unclear


I've come to realize that scripting ArcGIS 10.0 with Python imposes too many limitations for my goals with ArcGIS, so I've decided to bite the bullet and delve into ArcObjects.


Apparently, using ArcObjects in ArcGIS 10.0 has gotten somewhat easier because there is a new framework that avoids having to deal with COM objects. This seems good, but also means less examples for those starting at 10.0


The only detailed ESRI example of beginning with ArcObjects is a walkthrough of how to build an add-in. The instructions get somewhat vague in step 2a) and 3 of the section Adding ZoomToLayer functionality to the custom button:


< 2 > To implement ZoomToLayer functionality in your button, use the ArcGIS Snippet Finder tool.



  • a. Right-click the Visual Studio code editor window at the desired insertion point.

  • b, c, d



< 3 > Invoke the ZoomToLayer method upon OnClick() method implementation.


This is what my ZoomToLayerButton.cs file looks like before step 2:


using System;
using System.Collections.Generic;
using System.Text;
using System.IO;

namespace CustomUIElements
{
public class ZoomToLayerButton : ESRI.ArcGIS.Desktop.AddIns.Button

{
public ZoomToLayerButton()
{

}


protected override void OnClick()
{
ArcMap.Application.CurrentTool = null;

}
protected override void OnUpdate()
{
Enabled = ArcMap.Application != null;
}
}

}

and here is the snippet I grabbed:



 #region"Zoom to Active Layer in TOC"
// ArcGIS Snippet Title:
// Zoom to Active Layer in TOC
//
// Long Description:
// Zooms to the selected layer in the Table of Contents (TOC) associated with the active view.
//
// Add the following references to the project:
// ESRI.ArcGIS.ArcMapUI
// ESRI.ArcGIS.Carto

// ESRI.ArcGIS.Geometry
//
// Intended ArcGIS Products for this snippet:
// ArcGIS Desktop (ArcEditor, ArcInfo, ArcView)
//
// Applicable ArcGIS Product Versions:
// 9.2
// 9.3
// 9.3.1
// 10.0

//
// Required ArcGIS Extensions:
// (NONE)
//
// Notes:
// This snippet is intended to be inserted at the base level of a Class.
// It is not intended to be nested within an existing Method.
//

///Zooms to the selected layer in the TOC associated with the active view.

///
///An IMxDocument interface
///
///
public void ZoomToActiveLayerInTOC(ESRI.ArcGIS.ArcMapUI.IMxDocument mxDocument)
{
if (mxDocument == null)
{
return;
}

ESRI.ArcGIS.Carto.IActiveView activeView = mxDocument.ActiveView;

// Get the TOC
ESRI.ArcGIS.ArcMapUI.IContentsView IContentsView = mxDocument.CurrentContentsView;

// Get the selected layer
System.Object selectedItem = IContentsView.SelectedItem;
if (!(selectedItem is ESRI.ArcGIS.Carto.ILayer))
{
return;

}
ESRI.ArcGIS.Carto.ILayer layer = selectedItem as ESRI.ArcGIS.Carto.ILayer;


// Zoom to the extent of the layer and refresh the map
activeView.Extent = layer.AreaOfInterest;
activeView.Refresh();
}
#endregion


Regarding step 2.a)"Right-click the Visual Studio code editor window at the desired insertion point" where I should insert the snippet? Reading the comments in the snippet, I assume it should go after the first method:


public ZoomToLayerButton()
{

}


//here??

protected override void OnClick()

{
ArcMap.Application.CurrentTool = null;
}

Also, the comments in the snippet also say I should:


// Add the following references to the project:
// ESRI.ArcGIS.ArcMapUI
// ESRI.ArcGIS.Carto
// ESRI.ArcGIS.Geometry


Where do I add these? At the top of the ZoomTo LayerButton.cs file or my Config.Designer.cs file? (I have tried both of these locations and I get the same error: 'Carto' does not exist in the namespace ESRI.ArcGIS)


Regarding step 3, "Invoke the ZoomToLayer method upon OnClick() method implementation." I assume that the code shown later on after those instructions is what invokes the method:


protected override void OnClick()
{
ZoomToActiveLayerInTOC(ArcMap.Application.Document as IMxDocument);
}

I assume this also means I should erase the line that was in there before?


ArcMap.Application.CurrentTool = null;


I have not yet gone on to the next step of the walkthrough - Installing the custom button and tool - since I have not got the code so far to build. I knew I was in for a steep learning curve with ArcObjects, and this intro example certainly confirms this! Hopefully, with these clarifications, I will have my first ArcObjects code running soon enough.



Answer



Some clarifications:




  • You are still using COM objects, the object model on which ArcObjects is built. The add-in framework merely simplifies the registration of your COM object extensions into component categories and streamlining the install/uninstall process.


    There are a few things that you can't do with add-ins that you can only do with traditional COM object extensions, but typically add-ins will suffice for most desktop customization type projects.




  • References are set per-project in Visual Studio. Currently added references are listed under the "References" folder in the solution explorer. Right click References and click Add, then, under the .NET tab (not the COM tab), add the necessary assemblies.





Other than that it looks like what you're doing should work, although you will want to read up on object-oriented programming fundamentals like classes, interfaces, implementation and inheritance, GUI and user interface design, separation of concerns, etc.


In addition you'll want to start learning as much as you can about the practical basics: the .NET Framework, the Visual Studio IDE, generics, and language syntax (see the C# specification).


Perhaps on a more advanced level, you will do well to learn the ins and outs of COM (particularly marshaling and reference counting), serialization, threading, P/Invoke and the Windows API (for hacking things that ArcObjects or the add-in framework didn't provide).


Of course if you have questions about most of that you would want to look at stackoverflow.com instead of gis.se!


How to make a multi color map with one shapefile layer in QGIS?


I am working with a shapefile layer of the world. All the countries are in the Attribute table each one in a different row. I tried to make each polygon in a different color so the countries will be in a different color, but it did not work.enter image description here


In the Africa map I save each polygon separately and upload into a new QGIS. This method takes a lot of time. What would a better and faster method be?


enter image description here




Answer



You may right-click on the layer in the Layers Panel and go to Properties >> Style dialog.


Then, choose Categorized, set the field that stores the name of the countries and click the Classify button:


enter image description here


Finally, click the Apply button for applying the changes and you will see this:


enter image description here


qgis 2.10.1. mac - geotiff singleband pseudocolor



i've installed mac version from kyngchaos including all required packages ...
... and just noticed the following:


while windows version behaves like expected when trying to render an ordinary tif via "singleband pseudocolor" (see the image) qgis 2.10-win


the same procedure on the current mac-release returns: qgis 2.10-mac


interesting detail: pressing "apply" in mac version doesnt change anything and if you press "OK" and re-open layer properties the "render type"-setting has swiched back to "singleband grey".


even more interesting: "singleband grey" is faulty too: none of the options "color gradient", "min", "max", "contrast enhancement" does affect the way the tif is displayed, neither do the "load min/max values" settings !



anyone any ideas?




Sunday, 24 July 2016

Generating raster from isolines using SAGA or ArcGIS?


I have a vector dataset of mean annual precipitation polylines, and I would like to generate a raster from them. I'm not sure the best way to go about it.


ArcGIS has the Polyline to Raster tool but from reading its help file I get the impression it's not really designed to do what I need it to do. Obviously some kind of interpolation is required.


I am hoping that SAGA-GIS might have a useful module. Either that or I guess I could convert the lines to points and then krige them.




r - Transform coordinates from NOAA


I use data from NOAA for some analysis in R and I want to transform coordinates to EPSG:54004 or something useful.


There is something what they write about their coordinates.


...    

POS: 29-34
GEOPHYSICAL-POINT-OBSERVATION latitude coordinate
The latitude coordinate of a GEOPHYSICAL-POINT-OBSERVATION where southern
hemisphere is negative.
MIN: -90000 MAX: +90000
UNITS: Angular Degrees
SCALING FACTOR: 1000
DOM: A general domain comprised of the numeric characters (0-9), a plus
sign (+), and a minus sign (-).
+99999 = Missing


POS: 35-41
GEOPHYSICAL-POINT-OBSERVATION longitude coordinate
The longitude coordinate of a GEOPHYSICAL-POINT-OBSERVATION where values west from
000000 to 179999 are signed negative.
MIN: -179999 MAX: +180000 UNITS: Angular Degrees
SCALING FACTOR: 1000
DOM: A general domain comprised of the numeric characters (0-9), a plus
sign (+), and a minus sign (-).
+999999 = Missing

...

And my problem is, that I can't transform this coodrinates correctly. I use this R command:


    sc <- cbind(st$LAT, st$LON)
ptransform(sc/180*pi,
'+proj=latlong +ellps=sphere',
'+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs')

But coordinates are somehow wrong as you can see. The points should be in borders because they are meteostations from this countries. Map with points.



Answer




I solved the problem, it was in input to ptransform function in R.
Instead sc <- cbind(st$LAT, st$LON) I have to use sc <- cbind(st$LON, st$LAT). Then as you say I have to use EPSG:4326 as input CRS and EPSG:3857 as destination CRS. And the transformation with this command.


    tr <- ptransform(sc/180*pi, 
'+proj=longlat',
'+proj=merc')

Correct transformation


qgis - Using SVG from Wikipedia to fill polygon?


I'm trying to use the flag of each region in Italy to fill each region's polygon. I have downloaded the shapefile of Administrative areas for Italy from diva-gis and am now attempting to fill each of the 20 polygons with the corresponding flag image taken from Wikipedia at http://en.wikipedia.org/wiki/Regions_of_Italy


When you click on the flag image on the wikipedia page, and try to download a single flag, such as Sardinia:


http://en.wikipedia.org/wiki/Regions_of_Italy#mediaviewer/File:Flag_of_the_Italian_region_Sardinia.svg


you see that the images are already in svg format. However, when I attempt to use them in QGIS, all I see is blank space.



I downloaded Inkscape as I read somewhere this would be the appropriate tool for dealing with SVG, however opening the file in Inkscape and saving as SVG does not result in anything different.


There are 20 regions in Italy and therefore 20 flags, can someone lead me in the right direction?


My plan of attack so far is to change the Style to Categorized and classify the regions of Italy by the NAME_1 Column (which is the region name) and then change each symbol layer type to SVG fill. After that I was thinking I would download each SVG flag individually from wiki and chose it for the SVG fill by clicking on the ellipses (...) button and navigating to the appropriate file.



Answer



You have to create a folder called "svg" in your ".qgis2" directory. You should be able to find this in:


C:\Users\(your user name)\.qgis2

Then create the new svg folder, insert your flags in there and you should be able to see them in your Style > Symbol interface:


Symbol selector


For me, QGIS did tend to run a bit slow so it may be a good idea to first follow the answer provided by @hexamon.



GeoServer GeoWebCache seeding dies prematurely



I am using GeoServer 2.6 Final War File in Tomcat. I am setting up GeoWebCache and I am attempting to seed a table from Oracle that has about 600,000 records. (I used the Oracle NG provider to create the Data Source).
I select the following options from GeoServer under


"Tile Layers" >> "Seed/Truncate" for that layer:



  • Number of tasks to use: 01

  • Type of operation: Seed - generate missing tiles

  • Grid Set: EPSG:2236

  • Format: image/png


  • Zoom start: 00

  • Zoom stop: 18



I then click on Submit and see the following:


Id|Layer   |Status |Type|Estimated # of tiles|Tiles completed|Time elapsed |Time remaining|Tasks     
49|sf:MV0_1|RUNNING|SEED|2,857,446 |-1 |Estimating...|Estimating... |(Task 1 of 1)

Then I click on "Refresh list" and in about 10 seconds it is no longer there...I was trying to look in the logs folder and found no logging in:




C:\Program Files\Apache Software Foundation\Tomcat 8.0\logs



which is where GeoServer is installed.


Why is the seeding not working?




Here is the gridset 2236 defined in geowebcache.xml:



EPSG:2236
2236



306826.6397
36279.6630
982215.8007
2363339.4629



8192
4096

2048
1024
512
256
128
64
32
16
8
4

2
1
0.5
0.25
0.125
0.0625
0.03125

0.3048
256

256




Answering a question from comments:



"Does the preview of the integrated GWC show maps for you?"



I click on "Tile Layers", Then select the "EPSG: 2236 / png" to preview the cached tile. I get some tiles that display. (Only the ones generated before the process dies as mentioned above)





Answering another question from comments:



"You need to make sure that your Gridset in GeoWebCache matches what you have setup on the Layer in GeoServer."



Yes, the layer matches the gridset in GeoWebCache:


Data>>Layers>>Edit Layer>>Tile Caching>>Available Gridsets

One item shows:


Available gridsets 
Gridset: EPSG:2236

Published zoom levels: min/max
Cached zoom levels: min/max
Grid subset bounds: Dynamic

Also under "Tile Caching">>"Caching Defaults">>Default Cached Gridsets, I have only:


Gridset: EPSG:2236 
CRS : EPSG:2236
Tile Dimensions: 256 x 256
Zoom levels: 19
Disk Usage: 465.63 MB


Also under "Tile Caching">>"Disk Quota", disk quota is enabled and set to 1.953 GB.




Answering a question from comments:



"Do you have write permission for the logging folder?"



Yes, other log data is being updated in that folder such as :



C:\Program Files\Apache Software Foundation\Tomcat 8.0\logs




localhost*.log    
catalina*.log
etc



Answering a question from comments:



"Browse the folder and see if the images are stored? I remember that Geoserver couldn't read my folder size once, so I didn't think the images were cached but they were there when I browsed to the directory. Perhaps this is similar?"




Yes, there are some images maybe a handful which are created just before the seed tasks die. They are in the



C:\Program Files\Apache Software Foundation\Tomcat8.0\webapps\geoserver\data\gwc\



folder under subfolders.




I also want to mention that there is another oracle spatial table (setup under another layer) with same layout but with only 1300 records and that one does complete normally.




I found the error in the geoserver.log that is causing the gwc to stop:




org.geoserver.platform.ServiceException: This requested used more time than allowed and has been forcefully stopped. Max rendering time is 60.0s




Answer



I found the answer to the issue, although there is no available solution. It is a known bug where GeoWebCache stops forcefully. It looks like there is no current solution for bug # GEOS-6278:


https://jira.codehaus.org/browse/GEOS-6278


If a solution is available in the future for GeoWebCache, I will update this post.


Seeking point cloud (LiDAR) data?




Is there any freely available pointcloud data?


I'm especially interested in dense clouds of small areas(5m to 100m radius) but anything else would be fine too.


I've found that these pages have very nice and easy to access point cloud data.





Answer



I decided to merge other answers with mine and organize them into a tabular format. I think it is easier to read and manage for future visitors:


enter image description here


The table can be accessed from the following link in csv format:


View in tabular form: Free LiDar DataSources


Download (csv): Free LiDar DataSources


Please submit pull request if you intended to add to this list.


open source gis - Working with LiDAR data using other than Esri software?



I will be working with LiDAR data soon and I am wondering what the alternatives to using Esri software to do this are, including open source solutions?




Saturday, 23 July 2016

carto - Getting info from multiple rows in single popup of CartoDB?


If you have a table like this:


Name | Item | State ---------|--------|----------- Alice | Apple | Alabama Ben | Banana | Alabama Caitlin | Cherry | California


Assuming you have created polygons based on the states, how could you create an 'info window' which combines the details of all rows in the same state?



So for example, the popup for Alabama might be something like:


State: Alabama
People: 2
Items: Apple, Banana

As far as I can tell, the default behaviour for Cartodb would be to place one polygon on top of the other, and for the info window to display the info from the top most polygon only.




area - How does ST_Area in PostGIS work?


I am performing a simple calculation on a polygon known to be approximately 6226 km^2 in area. It is stored in a Geography (WGS84 SRID) column.


The query is:



select st_astext(col), st_area(col) area from table

and returns:


"POLYGON((-180 58.282525588539,-178.916399160189 57.4759784390599,-178.191728834624 58.5761461944577,-180 58.282525588539))" | 5807028547.33813

The area returned (5807028547.33813) appears to be mm^2 and not km^2? The documentation http://postgis.net/docs/ST_Area.html states "by default area is determined on a spheroid with units in square meters"


Is this a documentation error, or is the above correct and I'm fundamentally misunderstanding the functionality?



Answer



The calculation gives the right output. As you cited the documentation states "by default area is determined on a spheroid with units in square meters".


The result of your query is 5807028547.33813 m^2. To get the area in km^2 you have to devide the result by 1,000,000.



5807028547.33813 m^2 / 1,000,000 = 5807.02854733813 km^2


5807.02854733813 km^2 corresponds approximately to your expected 6226 km^2.


Friday, 22 July 2016

How to Join Attributes by Location in QGIS 1.8.0


I have a points layer with site attributes and a polygon layer with region attributes. I'd like to join the names of the regions to the point site locations using QGIS 1.8.0., on Windows XP. However, I keep generating an empty shapefile, with no records, but the polygon_region's field names are appended to the end of the points within the empty shapefile. I'm using the default settings for the "Join Attributes by Location" tool. The tool's progress bar indicator seems to advance forward to about 20% complete, then it instantly completes. I'm thinking that it's only processing the header of the shapefile then prematurely "completing"




Answer



My problem resolved using the Join Attributes by Points tool:


I had experienced a problem with the output (spatially joined) shapefile being written with header-column-names only, no data. It was skipping attribute creation even though the points geometries fell inside the polygon geometries. I had used and verified the same coordinate system for both layers and turned off On-The-Fly reprojection, to confirm they should join, and to see if the On-The-Fly option could be the culprit, but still no attributes were created for the join.


My solution was to export the data to a new CRS. I used the most common I knew of: [WGS84 (EPSG/SRID/ 4326)] Bingo--It worked!


If anyone experiences zero attribute creation in their spatial join, and they've made sure their CRS for both layers is exactly the same, and have already turned off On-The-Fly re-projection to confirm they should join, but it still fails, then export your data to WGS84 and try that CRS. That was the difference for me. Perhaps the tool has a problem with certain or custom CRS?


arcgis server - Specify geoprocessing service connection at publish


I have a geoprocessing service published from a tool written in Python and contained inside a Python toolbox.


I need to be able to publish this service in multiple environments (production, staging, testing, development) i.e. using different database backends.


I do not want to recreate the service definition package every time I need to move to a different environment.


In particular, I do not have direct access to production and do not want the people deploying the app to need to fire up ArcMap or ArcCatalog and risk some parameter not being set correctly.


Instead, I want to hand them a package and say "Publish this using the connection you created".


How can I accomplish that?



Edit:


Thanks to another question, I've run across this piece of documentation here:



If you are currently disconnected from the server or do not have access to a server connection, you can configure your service definition with no available connection to ArcGIS Server. Selecting this option creates a service definition file that will have to be configured to work with a server connection at the time of publishing.



So it seems to me that ESRI says this is possible, unless I've misinterpreted something.




Splitting shapefile into many shapefiles with open source?


I want to create a shapefile for each of the 72,000 census tracts in the United States that includes the census blocks in that tract using open source software.


I will probably start with a state shapefile that includes all of the blocks for that state and merge it with a GDB file. This will give me a list of blocks that are in the census tract.


But then I'm not sure how to split the state shapefile into smaller ones.


How can I do this?


Is there a way to automate this using open source software?


If so, what tools should I use?



I'm learning QGIS.



Answer



Basically, you can do the extract using ogr2ogr as long as you give the Census tract ID, so it's really an issue of getting 72,000 ogr2ogr calls.


ogr2ogr -where "tract = ''" /dest_folder /source_folder block_shapefile -nln block_shapefile_

Notes:



  • You don't have to specify the source format, ogr2ogr will figure it out. You specified shapefile, so that's what I'm assuming.

  • If you are not changing the format, you also don't have to specify the destination format. Otherwise, to get a shapefile, add `-f "ESRI Shapefile"

  • I'm using the -where switch to subset the data. Remember that attribute query is much faster than spatial query. Your question was a little vague, so I don't know if you intend to join the blocks to the tracts by attribute or spatially, but I highly recommend the former.


  • Note that I am assuming your blocks shapefile has a tract ID column in it. Most Census data sets will have a hierarchy of ID columns, i.e. the blocks shapefile will probably have a state, county, and tract, as well as block ID column, but may not have a column concatenating them all together. They must be concatenated, because counties are only unique with states, tracts are only unique with counties, and blocks are only unique within tracts. So you will have to either (a) create a new field with state & county & tract as your ID field, or add all three criteria to the where clause.


So how do you build 72,000 ogr2ogr calls programmatically? You can use any tool you want, but here's an example with R:


library(foreign)
dfBlocks = read.dbf("/source_folder/shapefile_name.dbf", as.is=TRUE)

strTract = unique(dfBlocks$tract_id)
for (i in 1:length(strTract)) {
strOGR = paste(
"ogr2ogr -where \"tract_id = '", strTract,

"'\" /dest_folder /source_folder layer_name -nln base_name_",
strTract, sep=""
)
system(strOGR)
}

You could also collect the system calls in one step, then iterate it to run the ogr2ogr call in a separate loop, or at a later time, or write it to a bash script that you run from the command line.


The major disadvantage is that it will scan the data source once for each ogr2ogr call, so actually importing the source data, iterating, and writing a shapefile for each row, would probably be more efficient. But I would recommend trying it in something other than R, which is somewhat slow at reading large spatial datasets (I ran an import just of the counties of the US, ~3000, and after 15 minutes I cancelled the import.)


arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...