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