From ac3a760ea7def0977163c2ecb01afbc4cc4bd9af Mon Sep 17 00:00:00 2001 From: Sean Gillies Date: Sat, 8 Mar 2014 14:49:31 -0700 Subject: [PATCH] Finish reprojection doc. --- docs/reproject.rst | 49 +++++++++++++++++++++++++++++++++++++++---- examples/reproject.py | 12 +++++++---- rasterio/warp.py | 2 +- 3 files changed, 54 insertions(+), 9 deletions(-) diff --git a/docs/reproject.rst b/docs/reproject.rst index 5cc77e88..23aab695 100644 --- a/docs/reproject.rst +++ b/docs/reproject.rst @@ -5,13 +5,54 @@ Rasterio can map the pixels of a destination raster with an associated coordinat reference system and transform to the pixels of a source image with a different coordinate reference system and transform. This process is known as reprojection. -Rasterio's ``rasterio.warp.reproject()`` is a very geospatial specific analog to +Rasterio's ``rasterio.warp.reproject()`` is a very geospatial-specific analog to SciPy's ``scipy.ndimage.interpolation.geometric_transform()`` [1]_. -Result ------- +The code below reprojects between two arrays, using no pre-existing GIS datasets. +``rasterio.warp.reproject()`` has two positional arguments: source and destination. +The remaining keyword arguments parameterize the reprojection transform. -https://a.tiles.mapbox.com/v3/sgillies.hfek2oko/page.html?secure=1#6/0.000/0.033 +.. code-block:: python + + import numpy + import rasterio + from rasterio.warp import reproject, RESAMPLING + + with rasterio.drivers(): + + # Consider a 512 x 512 raster centered on 0 degrees E and 0 degrees N + # with each pixel covering 15". + src_shape = (512, 512) + src_transform = [-256.0/240, 1.0/240, 0.0, 256.0/240, 0.0, -1.0/240] + src_crs = {'init': 'EPSG:4326'} + source = numpy.ones(src_shape, numpy.uint8)*255 + + # Prepare to reproject this rasters to a 1024 x 1024 dataset in + # Web Mercator (EPSG:3857) with origin at 0.0, 0.0. + dst_shape = (1024, 1024) + dst_transform = [-237481.5, 425.0, 0.0, 237536.4, 0.0, -425.0] + dst_crs = {'init': 'EPSG:3857'} + destination = numpy.zeros(dst_shape, numpy.uint8) + + reproject( + source, + destination, + src_transform=src_transform, + src_crs=src_crs, + dst_transform=dst_transform, + dst_crs=dst_crs, + resampling=RESAMPLING.nearest) + + # Assert that the destination is only partly filled. + assert destination.any() + assert not destination.all() + +See `examples/reproject.py `__ for code that writes the destination array to a GeoTIFF file. I've +uploaded the resulting file to a Mapbox map to demonstrate that the reprojection is +correct: https://a.tiles.mapbox.com/v3/sgillies.hfek2oko/page.html?secure=1#6/0.000/0.033 + +References +---------- .. [1] http://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.interpolation.geometric_transform.html#scipy.ndimage.interpolation.geometric_transform diff --git a/examples/reproject.py b/examples/reproject.py index a6e61d6f..a938851d 100644 --- a/examples/reproject.py +++ b/examples/reproject.py @@ -5,7 +5,7 @@ import tempfile import numpy import rasterio -from rasterio.warp import reproject +from rasterio.warp import reproject, RESAMPLING tempdir = '/tmp' tiffname = os.path.join(tempdir, 'example.tif') @@ -27,9 +27,13 @@ with rasterio.drivers(): destination = numpy.zeros(dst_shape, numpy.uint8) reproject( - source, destination, - src_transform, src_crs, - dst_transform, dst_crs) + source, + destination, + src_transform=src_transform, + src_crs=src_crs, + dst_transform=dst_transform, + dst_crs=dst_crs, + resampling=RESAMPLING.nearest) # Assert that the destination is only partly filled. assert destination.any() diff --git a/rasterio/warp.py b/rasterio/warp.py index 4cfc3a69..b0ee65d9 100644 --- a/rasterio/warp.py +++ b/rasterio/warp.py @@ -1,4 +1,4 @@ """Raster warping and reprojection""" import rasterio -from rasterio._warp import reproject +from rasterio._warp import reproject, RESAMPLING