1Virtual Warping
2===============
3
4Rasterio has a ``WarpedVRT`` class that abstracts many of the details of raster
5warping by using an in-memory `Warped VRT
6<http://www.gdal.org/gdal_vrttut.html#gdal_vrttut_warped>`__. A ``WarpedVRT`` can
7be the easiest solution for tiling large datasets.
8
9For example, to virtually warp the ``RGB.byte.tif`` test dataset from its
10proper EPSG:32618 coordinate reference system to EPSG:3857 (Web Mercator) and
11extract pixels corresponding to its central zoom 9 tile, do the following.
12
13.. code-block:: python
14
15  from affine import Affine
16  import mercantile
17
18  import rasterio
19  from rasterio.enums import Resampling
20  from rasterio.vrt import WarpedVRT
21
22  with rasterio.open('tests/data/RGB.byte.tif') as src:
23      with WarpedVRT(src, crs='EPSG:3857',
24                     resampling=Resampling.bilinear) as vrt:
25
26          # Determine the destination tile and its mercator bounds using
27          # functions from the mercantile module.
28          dst_tile = mercantile.tile(*vrt.lnglat(), 9)
29          left, bottom, right, top = mercantile.xy_bounds(*dst_tile)
30
31          # Determine the window to use in reading from the dataset.
32          dst_window = vrt.window(left, bottom, right, top)
33
34          # Read into a 3 x 512 x 512 array. Our output tile will be
35          # 512 wide x 512 tall.
36          data = vrt.read(window=dst_window, out_shape=(3, 512, 512))
37
38          # Use the source's profile as a template for our output file.
39          profile = vrt.profile
40          profile['width'] = 512
41          profile['height'] = 512
42          profile['driver'] = 'GTiff'
43
44          # We need determine the appropriate affine transformation matrix
45          # for the dataset read window and then scale it by the dimensions
46          # of the output array.
47          dst_transform = vrt.window_transform(dst_window)
48          scaling = Affine.scale(dst_window.height / 512,
49                                 dst_window.width / 512)
50          dst_transform *= scaling
51          profile['transform'] = dst_transform
52
53          # Write the image tile to disk.
54          with rasterio.open('/tmp/test-tile.tif', 'w', **profile) as dst:
55              dst.write(data)
56
57
58Normalizing Data to a Consistent Grid
59-------------------------------------
60
61A ``WarpedVRT`` can be used to normalize a stack of images with differing
62projections, bounds, cell sizes, or dimensions against a regular grid
63in a defined bounding box.
64
65The `tests/data/RGB.byte.tif` file is in UTM zone 18, so another file in a
66different CRS is required for demonstration.  This command will create a new
67image with drastically different dimensions and cell size, and reproject to
68WGS84.  As of this writing ``$ rio warp`` implements only a subset of
69`$ gdalwarp <http://www.gdal.org/gdalwarp.html>`__'s features, so
70``$ gdalwarp`` must be used to achieve the desired transform:
71
72.. code-block:: console
73
74    $ gdalwarp \
75        -t_srs EPSG:4326 \
76        -te_srs EPSG:32618 \
77        -te 101985 2673031 339315 2801254 \
78        -ts 200 250 \
79        tests/data/RGB.byte.tif \
80        tests/data/WGS84-RGB.byte.tif
81
82So, the attributes of these two images drastically differ:
83
84.. code-block:: console
85
86    $ rio info --shape tests/data/RGB.byte.tif
87    718 791
88    $ rio info --shape tests/data/WGS84-RGB.byte.tif
89    250 200
90    $ rio info --crs tests/data/RGB.byte.tif
91    EPSG:32618
92    $ rio info --crs tests/data/WGS84-RGB.byte.tif
93    EPSG:4326
94    $ rio bounds --bbox --geographic --precision 7 tests/data/RGB.byte.tif
95    [-78.95865, 23.5649912, -76.5749237, 25.5508738]
96    $ rio bounds --bbox --geographic --precision 7 tests/data/WGS84-RGB.byte.tif
97    [-78.9147773, 24.119606, -76.5963819, 25.3192311]
98
99and this snippet demonstrates how to normalize data to consistent dimensions,
100CRS, and cell size within a pre-defined bounding box:
101
102.. code-block:: python
103
104    from __future__ import division
105
106    import os
107
108    import affine
109
110    import rasterio
111    from rasterio.crs import CRS
112    from rasterio.enums import Resampling
113    from rasterio import shutil as rio_shutil
114    from rasterio.vrt import WarpedVRT
115
116
117    input_files = (
118        # This file is in EPSG:32618
119        'tests/data/RGB.byte.tif',
120        # This file is in EPSG:4326
121        'tests/data/WGS84-RGB.byte.tif'
122    )
123
124    # Destination CRS is Web Mercator
125    dst_crs = CRS.from_epsg(3857)
126
127    # These coordiantes are in Web Mercator
128    dst_bounds = -8744355, 2768114, -8559167, 2908677
129
130    # Output image dimensions
131    dst_height = dst_width = 100
132
133    # Output image transform
134    left, bottom, right, top = dst_bounds
135    xres = (right - left) / dst_width
136    yres = (top - bottom) / dst_height
137    dst_transform = affine.Affine(xres, 0.0, left,
138                                  0.0, -yres, top)
139
140    vrt_options = {
141        'resampling': Resampling.cubic,
142        'crs': dst_crs,
143        'transform': dst_transform,
144        'height': dst_height,
145        'width': dst_width,
146    }
147
148    for path in input_files:
149
150        with rasterio.open(path) as src:
151
152            with WarpedVRT(src, **vrt_options) as vrt:
153
154                # At this point 'vrt' is a full dataset with dimensions,
155                # CRS, and spatial extent matching 'vrt_options'.
156
157                # Read all data into memory.
158                data = vrt.read()
159
160                # Process the dataset in chunks.  Likely not very efficient.
161                for _, window in vrt.block_windows():
162                    data = vrt.read(window=window)
163
164                # Dump the aligned data into a new file.  A VRT representing
165                # this transformation can also be produced by switching
166                # to the VRT driver.
167                directory, name = os.path.split(path)
168                outfile = os.path.join(directory, 'aligned-{}'.format(name))
169                rio_shutil.copy(vrt, outfile, driver='GTiff')
170