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