1Reprojection 2============ 3 4Rasterio can map the pixels of a destination raster with an associated 5coordinate reference system and transform to the pixels of a source image with 6a different coordinate reference system and transform. This process is known as 7reprojection. 8 9Rasterio's :func:`rasterio.warp.reproject()` is a geospatial-specific analog 10to SciPy's ``scipy.ndimage.interpolation.geometric_transform()`` [1]_. 11 12The code below reprojects between two arrays, using no pre-existing GIS 13datasets. :func:`rasterio.warp.reproject()` has two positional arguments: source 14and destination. The remaining keyword arguments parameterize the reprojection 15transform. 16 17.. code-block:: python 18 19 import numpy as np 20 import rasterio 21 from rasterio import Affine as A 22 from rasterio.warp import reproject, Resampling 23 24 with rasterio.Env(): 25 26 # As source: a 512 x 512 raster centered on 0 degrees E and 0 27 # degrees N, each pixel covering 15". 28 rows, cols = src_shape = (512, 512) 29 d = 1.0/240 # decimal degrees per pixel 30 # The following is equivalent to 31 # A(d, 0, -cols*d/2, 0, -d, rows*d/2). 32 src_transform = A.translation(-cols*d/2, rows*d/2) * A.scale(d, -d) 33 src_crs = {'init': 'EPSG:4326'} 34 source = np.ones(src_shape, np.uint8)*255 35 36 # Destination: a 1024 x 1024 dataset in Web Mercator (EPSG:3857) 37 # with origin at 0.0, 0.0. 38 dst_shape = (1024, 1024) 39 dst_transform = A.translation(-237481.5, 237536.4) * A.scale(425.0, -425.0) 40 dst_crs = {'init': 'EPSG:3857'} 41 destination = np.zeros(dst_shape, np.uint8) 42 43 reproject( 44 source, 45 destination, 46 src_transform=src_transform, 47 src_crs=src_crs, 48 dst_transform=dst_transform, 49 dst_crs=dst_crs, 50 resampling=Resampling.nearest) 51 52 # Assert that the destination is only partly filled. 53 assert destination.any() 54 assert not destination.all() 55 56 57See `examples/reproject.py <https://github.com/mapbox/rasterio/blob/master/examples/reproject.py>`__ 58for code that writes the destination array to a GeoTIFF file. I've uploaded the 59resulting file to a Mapbox map to show that the reprojection is 60correct: https://a.tiles.mapbox.com/v3/sgillies.hfek2oko/page.html?secure=1#6/0.000/0.033. 61 62Estimating optimal output shape 63------------------------------- 64 65Rasterio provides a :func:`rasterio.warp.calculate_default_transform()` function to 66determine the optimal resolution and transform for the destination raster. 67Given a source dataset in a known coordinate reference system, this 68function will return a ``transform, width, height`` tuple which is calculated 69by libgdal. 70 71Reprojecting a GeoTIFF dataset 72------------------------------ 73 74Reprojecting a GeoTIFF dataset from one coordinate reference system is a common 75use case. Rasterio provides a few utilities to make this even easier: 76 77:func:`~rasterio.warp.transform_bounds()` 78transforms the bounding coordinates of the source raster to the target 79coordinate reference system, densifiying points along the edges to account 80for non-linear transformations of the edges. 81 82 83:func:`~rasterio.warp.calculate_default_transform()` 84transforms bounds to target coordinate system, calculates resolution if not 85provided, and returns destination transform and dimensions. 86 87 88.. code-block:: python 89 90 import numpy as np 91 import rasterio 92 from rasterio.warp import calculate_default_transform, reproject, Resampling 93 94 dst_crs = 'EPSG:4326' 95 96 with rasterio.open('rasterio/tests/data/RGB.byte.tif') as src: 97 transform, width, height = calculate_default_transform( 98 src.crs, dst_crs, src.width, src.height, *src.bounds) 99 kwargs = src.meta.copy() 100 kwargs.update({ 101 'crs': dst_crs, 102 'transform': transform, 103 'width': width, 104 'height': height 105 }) 106 107 with rasterio.open('/tmp/RGB.byte.wgs84.tif', 'w', **kwargs) as dst: 108 for i in range(1, src.count + 1): 109 reproject( 110 source=rasterio.band(src, i), 111 destination=rasterio.band(dst, i), 112 src_transform=src.transform, 113 src_crs=src.crs, 114 dst_transform=transform, 115 dst_crs=dst_crs, 116 resampling=Resampling.nearest) 117 118 119See ``rasterio/rio/warp.py`` for more complex examples of reprojection based on 120new bounds, dimensions, and resolution (as well as a command-line interface 121described :ref:`here <warp>`). 122 123It is also possible to use :func:`~rasterio.warp.reproject()` to create an output dataset zoomed 124out by a factor of 2. Methods of the :class:`rasterio.Affine` class help us generate 125the output dataset's transform matrix and, thereby, its spatial extent. 126 127.. code-block:: python 128 129 import numpy as np 130 import rasterio 131 from rasterio import Affine as A 132 from rasterio.warp import reproject, Resampling 133 134 with rasterio.open('rasterio/tests/data/RGB.byte.tif') as src: 135 src_transform = src.transform 136 137 # Zoom out by a factor of 2 from the center of the source 138 # dataset. The destination transform is the product of the 139 # source transform, a translation down and to the right, and 140 # a scaling. 141 dst_transform = src_transform*A.translation( 142 -src.width/2.0, -src.height/2.0)*A.scale(2.0) 143 144 data = src.read() 145 146 kwargs = src.meta 147 kwargs['transform'] = dst_transform 148 149 with rasterio.open('/tmp/zoomed-out.tif', 'w', **kwargs) as dst: 150 151 for i, band in enumerate(data, 1): 152 dest = np.zeros_like(band) 153 154 reproject( 155 band, 156 dest, 157 src_transform=src_transform, 158 src_crs=src.crs, 159 dst_transform=dst_transform, 160 dst_crs=src.crs, 161 resampling=Resampling.nearest) 162 163 dst.write(dest, indexes=i) 164 165.. image:: https://farm8.staticflickr.com/7399/16390100651_54f01b8601_b_d.jpg) 166 167Reprojecting with other georeferencing metadata 168------------------------------------------------ 169 170Most geospatial datasets have a geotransform which can be used to reproject a dataset 171from one coordinate reference system to another. Datasets may also be 172georeferenced by alternative metadata, namely Ground Control Points (gcps) or 173Rational Polynomial Coefficients (rpcs). For details on gcps and rpcs, see 174:doc:`georeferencing`. A common scenario is using gcps or rpcs to geocode 175(orthorectify) datasets, resampling and reorienting them to a coordinate 176reference system with a newly computed geotransform. 177 178.. code-block:: python 179 180 import rasterio 181 from rasterio.warp import reproject 182 from rasterio.enums import Resampling 183 184 with rasterio.open('RGB.byte.tif') as source: 185 print(source.rpcs) 186 src_crs = "EPSG:4326" # This is the crs of the rpcs 187 188 # Optional keyword arguments to be passed to GDAL transformer 189 # https://gdal.org/api/gdal_alg.html?highlight=gdalcreategenimgprojtransformer2#_CPPv432GDALCreateGenImgProjTransformer212GDALDatasetH12GDALDatasetHPPc 190 kwargs = { 191 'RPC_DEM': '/path/to/dem.tif' 192 } 193 194 # Destination: a 1024 x 1024 dataset in Web Mercator (EPSG:3857) 195 destination = np.zeros((1024, 1024), dtype=np.uint8) 196 dst_crs = "EPSG:3857" 197 198 _, dst_transform = reproject( 199 source, 200 destination, 201 rpcs=source.rpcs, 202 src_crs=src_crs, 203 dst_crs=dst_crs, 204 resampling=Resampling.nearest, 205 **kwargs 206 ) 207 208 assert destination.any() 209 210.. note:: 211 When reprojecting a dataset with gcps or rpcs, the src_crs parameter should 212 be supplied with the coordinate reference system that the gcps or rpcs are 213 referenced against. By definition rpcs are always referenced against WGS84 214 ellipsoid with geographic coordinates (EPSG:4326). 215 216 217 218References 219---------- 220 221.. [1] https://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.geometric_transform.html#scipy.ndimage.geometric_transform 222