2014 (old posts, page 2)

The GeoJSON media type

It's official: application/vnd.geo+json is the media type for GeoJSON registered with the IANA. I had been holding out for a profile of application/json, but it seems that the internet has resigned itself to a Cambrian explosion of JSON media types.

Start serving your GeoJSON up with the new content type

HTTP/1.1 200 OK
Content-Type: application/vnd.geo+json

and start asking for it by name.

Accept: application/vnd.geo+json; q=1.0, application/json; q=0.5

The GeoJSON-LD project

JSON-LD is JSON for linking data and it's a web standard: http://json-ld.org/. It offers a smooth, even zero-edit in some cases, path for upgrading existing GeoJSON data and makes web-friendly extension of GeoJSON possible.

Previously, if you added a "timestamp": "2014-03-27" to a GeoJSON feature object's "properties", a user has to guess what you mean. Playing guessing games with my kids is fun. Playing pointless guessing games with data sucks. One of the most important things JSON-LD offers is a framework in which you can unambiguously identify "timestamp" as "http://purl.org/dc/terms/created". Now users of your data have a solid start toward understanding what you mean by "timestamp".

There's more to JSON-LD than defining your feature properties. You can also bring in terms from other vocabularies to express that your features are, say, schema.org AdministrativeAreas or Landforms.

If you'd like to help shape the way we use GeoJSON in the future, please follow https://github.com/geojson/geojson-ld and jump into discussion on the open vocabulary and context issues.

Warping images with rasterio

Warping imagery from one spatial basis to another is a crucial function. I've been determined that rasterio be able to do this. Scipy provides a geometric_transform() function that I like, but I'd like rasterio to operate at a higher level and take advantage of GDAL's fast transformers. GDAL's gdalwarp program only operates on files, and I'm aiming lower than this (so rasterio can let us build better scripts). The ReprojectImage() function in osgeo.gdal, as pointed out in http://jgomezdans.github.io/gdal_notes/reprojection.html, only works on GDAL dataset objects. And the code to make them gets messy.

My requirements for a rasterio warper are:

  • Speed approaching gdalwarp.
  • Same results (same algorithms) as gdalwarp.
  • Sources and destinations can be ndarrays as well as bands in files.
  • Simplicity and clarity.

I think I've managed to satisfy these requirements. The warp module will be the key feature of the next rasterio release.

Here's an example within the interactive dataset inspector.

$ rasterio.insp rasterio/tests/data/RGB.byte.tif
Rasterio 0.7 Interactive Inspector (Python 2.7.5)
Type "src.name", "src.read_band(1)", or "help(src)" for more information.

The reproject() function takes source and destination ndarrays and an affine transformation vector and coordinate reference system for each. In this example, I'm warping the UTM source to a Web Mercator destination of slightly smaller extent. I computed the extents of the destination using pyplot, rasterio has yet no method to come up with an optimal destination extent (to do).

>>> source = src.read_band(1)
>>> dst_transform = [
...     -8789636.708,
...     300.0,
...     0.0,
...     2943560.235,
...     0.0,
...     -300.0]
>>> dst_crs = {'init': 'EPSG:3857'}
>>> destin = np.empty(src.shape, dtype=np.uint8)
>>> from rasterio import _warp
>>> _warp.reproject(
...     source, destin,
...     src.transform, src.crs,
...     dst_transform, dst_crs)
>>> import matplotlib.pyplot as plt
>>> plt.subplot(121)
<matplotlib.axes.AxesSubplot object at 0x1075362d0>
>>> plt.imshow(source)
<matplotlib.image.AxesImage object at 0x1078b1450>
>>> plt.gray()
>>> plt.subplot(122)
<matplotlib.axes.AxesSubplot object at 0x1075ffe50>
>>> plt.imshow(destin)
<matplotlib.image.AxesImage object at 0x1078b19d0>
>>> plt.gray()
>>> plt.savefig('/tmp/warp.png')

It's pretty much just three statements. One to get a source ndarray from a dataset. One to create a destination ndarray. One to call reproject().

I'm out of time to properly transform the axes in the plots to show coordinates in their proper CRS, but they're warped correctly.


Rasterio cookbook recipe number one

I got a nice message today from a developer who found rasterio to be the right tool for scaling the values of a 15 GB LiDAR DEM. Here's the code.

import numpy as np
import rasterio

Bryan Luman
Use it however you like at your own risk

I have a huge DEM converted from LiDAR LAS points. I'd like to make
it slightly less massive for further processing.

All my z values are between 682 and 869 feet stored as float32. If
I convert them to centimeters and they will fit well within a 16 bit
unsigned integer.

This saves a significant amount of storage space and greatly reduces
any type of future analysis.


CM_IN_FOOT = 30.48

with rasterio.drivers():
    with rasterio.open('massive_dem_as_float') as src:
        kwargs = src.meta
            bigtiff='YES' # Output will be larger than 4GB

        windows = src.block_windows(1)

        with rasterio.open(
                **kwargs) as dst:
            for idx, window in windows:
                src_data = src.read_band(1, window=window)

                # Source nodata value is a very small negative number
                # Converting in to zero for the output raster
                np.putmask(src_data, src_data < 0, 0)

                dst_data = (src_data * CM_IN_FOOT).astype(rasterio.uint16)
                dst.write_band(1, dst_data, window=window)

This uses the windowed read and write and windows iterator features of rasterio that I worked on with Asger Skovbo Petersen.

Do you have any more short rasterio code examples? I'd love to start collecting these into a cookbook of sorts.