Using GeoJSON from an ESRI REST Service in Leaflet

I love Leaflet for making web maps but sometimes linking Leaflet with external spatial data resources, particularly those from ESRI, is difficult at best and almost impossible at worst. Tools like ESRI-Leaflet (https://esri.github.io/esri-leaflet/) have been a life saver in most cases, particularly for using ESRI basemaps. However, recently, I needed to use a Web Feature Service from an ArcGIS REST service to display data in leaflet.  It took a while and I thought I'd share the results.


Understanding spatial REST services endpoints
First, all RESTful services terminate in an 'endpoint'.  That is the final URL to access a particular layer of spatial data.  For example, here is the ESRI REST endpoint for the US National Weather Service's Watches and Warnings data.

https://idpgis.ncep.noaa.gov/arcgis/rest/services/NWS_Forecasts_Guidance_Warnings/watch_warn_adv/MapServer/1

Depending on their configuration, ArcGIS Server-based REST layers can expose a 'query' interface. If available, it's linked on the bottom of the endpoint bottom.  Look for this text:

Supported Operations:   Query   Generate Renderer   Return Update

Clicking the Query page will bring you to an intimidating Query form.  However, the form itself can be bypassed using an HTTP POST.

https://idpgis.ncep.noaa.gov/arcgis/rest/services/NWS_Forecasts_Guidance_Warnings/watch_warn_adv/MapServer/1/query?where=phenom%3D%27FW%27&outFields=*&f=geojson

Basically, you call the query endpoint and pass a 'where' clause (in this case it's: where=phenom='FW').  The URL will get mangled to replace the = (%3D) and ' (%27) but basically calling with a where clause and passing it an outFields=* will give you a filtered set of elements and all their data fields. Finally, including the "f=geojson" will return the results in geojson format.

Getting GeoJSON from the REST service into Leaflet
In leaflet, just define a variable to to that url something like:

myURL = "https://idpgis.ncep.noaa.gov/arcgis/rest/services/NWS_Forecasts_Guidance_Warnings/watch_warn_adv/MapServer/1/query?where=phenom%3D%27FW%27&outFields=*&f=geojson";

The layer itself is loaded using JQueries AJAX functionality as follow:

var geojsonLayer;
$.ajax({
dataType: "json",
url: myURL,
success: function(data) {
geojsonLayer = new L.GeoJSON(data,{style: style}).addTo(map);}
}).error(function() {});


function style() {
return {
weight: 1,opacity: 0.65,dashArray: '1',fillOpacity: 0.65,color: '#000000',fillColor: '#d53e4f'
};
}

This will call the style() function for each GeoJSON element loaded and then add the result to a Leaflet map called 'map'.

Then you can add and use the geojsonLayer as needed within Leaflet.  There are lots of other things you can do with that Layer, like zoom on click, highlighting, etc.., and those are well described in the Leaflet Chloropleth tutorial here: https://leafletjs.com/examples/choropleth/

Iterate through POST in PHP

When linking PHP scripts with HTML Forms, sometimes it is useful to iterate through all the elements passed to the script from the POST event.  Here is a snippet of code that does that:

foreach($_POST as $key=>$value)
{
echo "$key: $value\n";
}
php?>

This little code snippet can really come in handy when debugging a PHP script that receives data by POST.

Obtaining OGR driver information

Like GDAL, OGR requires the use of a 'driver' to properly format an imported or exported file. To determine what drivers are available as part of your your OGR installation, type the following at the command line:

ogrinfo --formats

Here are the results for my system running MS4W:

Supported Formats:
-> "ESRI Shapefile" (read/write)
-> "MapInfo File" (read/write)
-> "UK .NTF" (readonly)
-> "SDTS" (readonly)
-> "TIGER" (read/write)
-> "S57" (read/write)
-> "DGN" (read/write)
-> "VRT" (readonly)
-> "REC" (readonly)
-> "Memory" (read/write)
-> "BNA" (read/write)
-> "CSV" (read/write)
-> "GML" (read/write)
-> "GPX" (read/write)
-> "KML" (read/write)
-> "GeoJSON" (read/write)
-> "GMT" (read/write)
-> "SQLite" (read/write)
-> "ODBC" (read/write)
-> "PGeo" (readonly)
-> "PostgreSQL" (read/write)
-> "MySQL" (read/write)
-> "XPlane" (readonly)
-> "AVCBin" (readonly)
-> "AVCE00" (readonly)
-> "Geoconcept" (read/write)

Getting started with Python-based open source geoprocessing

Well, I've had this Blog configured for some time and never found the inspiration that I needed to get started until today. In my profession, I seem to be constantly evaluating new geospatial technologies to find the right tool for the right job and as it would happen, many of these tools lack documentation so I am left to my own vices. Often I can find a good tutorial or some example code but sometimes I can't. My new-found inspiration for this blog is to share my misery with others and hopefully help others find the information they need to make progress.

This week's exploration is all about the SWIG-based Python bindings for GDAL or the Geospatial Data Abstraction Library (http://gdal.org/). GDAL is a miracle tool for handling raster-based spatial data. You can translate data between various formats such as GeoTIFF, GRID, etc..., and you can also reproject the rasters to a new coordinate system. Combining GDAL and Python is a natural marriage because it allows you to script these processes to simplify repetitive tasks. The main problem is that there is little documentation for beginners such as myself. The following is a simple script that I have been working through slowly to learn the ins and outs of the Python GDAL interface. I know it must be painfully simple for some but this is how I work things out. As I learn more, I'll keep posting things here so that others can learn along with me.

# Import the gdal and osr modules
from osgeo import gdal,osr

# GDAL works on the concept of the "Driver" which automagically recognizes the type of file you are trying to open.
# The simplest method is to register all the drivers so that it has the maximum amount of options for files that you are opening.

gdal.AllRegister()

# Open a new GDAL dataset
dataset = gdal.Open("c://tmp//gdal//usgsicurrent.bil")
if dataset is None:
print "Could not open dataset"

# Print some information about the file
print 'Driver: ', dataset.GetDriver().ShortName,'/',dataset.GetDriver().LongName
print 'Size is ',dataset.RasterXSize,'x',dataset.RasterYSize,'x',dataset.RasterCount
print 'Projection is ',dataset.GetProjection()

# Get some information about the image like the corner coordinates and cell sizes
geotransform = dataset.GetGeoTransform()
if not geotransform is None:
print 'Origin = (',geotransform[0], ',',geotransform[3],')'
print 'Pixel Size = (',geotransform[1], ',',geotransform[5],')'

xOrigin = geotransform[0]
yOrigin = geotransform[3]
pixelWidth = geotransform[1]
pixelHeight = geotransform[5]

# Create new Spatial References using OSR
outproj = osr.SpatialReference()
ingeog = osr.SpatialReference()

# Set the output projection to the US National Atlas Lambert Equal Area Azimuthal Projection
outproj.ImportFromEPSG(2163)

# Set the input projection to WGS84
ingeog.SetWellKnownGeogCS("WGS84")

# Use the ExportToWkt() to set the input raster's coordinate system using the spatial reference created earlier.
dataset.SetProjection(outproj.ExportToWkt())

# Create a new Coordinate Transformation object
ct = osr.CoordinateTransformation(ingeog, outproj)

# Transform a point from Geographic to Projected Coordinates
(MyX,MyY,MyZ) = ct.TransformPoint(-100.0,45.00)

# Compute the offset
xOffset = int((MyX - xOrigin) / pixelWidth)
yOffset = int((MyY - yOrigin) / pixelHeight)

# Read the raster value for the tranformed coordinate
band = dataset.GetRasterBand(1) # 1-based index
data = band.ReadAsArray(xOffset, yOffset, 1, 1)
value = data[0,0]

# Print the value
print value

I know, it doesn't seem like much but this little program has many uses and the individual components are are useful are the whole. As I learn more and have more examples, I'll continue to share them.

This site helped with the coordinate system tranformation: http://help.maptiler.org/coordinates/batch-transformation

Here is a GREAT tutorial on Python-based geoprocessing using GDAL/OGR: http://www.gis.usu.edu/~chrisg/python/

Here is the official GDAL API tutorial with some information about the Python implementation: http://www.gdal.org/gdal_tutorial.html