Reprojecting a large raster using Python and rasterio

Question

I'm trying to reproject a raster that describes land cover over time in the Brazilian Amazon. The raster can be downloaded at http://storage.googleapis.com/mapbiomas-public/COLECAO/2_3/CONSOLIDACAO/AMAZONIA.tif. It is roughly 2.1 gigabytes on disk. The metadata looks like this (there is one band per year):

{'count': 17, 'crs': CRS({'init': u'epsg:4326'}), 'dtype': 'uint8', 'affine': Affine(0.0002694945852358565, 0.0, -75.00007357655362,
       0.0, -0.0002694945852358565, 6.000027445691108), 'driver': u'GTiff', 'transform': (-75.00007357655362, 0.0002694945852358565, 0.0, 6.000027445691108, 0.0, -0.0002694945852358565), 'height': 90112, 'width': 123904, 'nodata': None}

You can visualize the data at http://mapbiomas.org/map.

The raster is originally in http://spatialreference.org/ref/epsg/wgs-84/; I'd like to reproject it to an equal area projection in meters, e.g. http://spatialreference.org/ref/sr-org/brazil-albers-equal-area-conic-wgs84/proj4/, at a resolution similar to the original (I think 30 meters would do).

My personal laptop has only 6 gigabytes of RAM. Can I reproject this raster using rasterio and Python?

I tried the following code, but it requires me to load the original raster values into RAM, and I get a memory error:

import numpy as np
import pyproj
import rasterio
from rasterio.warp import reproject, Resampling
from rasterio import crs, transform

def main(infile_path='./mapbiomass/AMAZONIA.tif',
         outfile_path='./mapbiomass/AMAZONIA_reprojected.tif',
         origin_lat=-11.188404,
         origin_lon= -58.338657,
         resolution=30.0):
    destination_crs_string = ' '.join(['+proj=aea +lat_1=10 +lat_2=-40 +lat_0=-25 +lon_0=-50',
                                       '+x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_def'])
    pyproj.transform(pyproj.Proj(init='EPSG:4326'), pyproj.Proj(destination_crs_string), origin_lon, origin_lat)
    origin_x, origin_y = pyproj.transform(pyproj.Proj(init='EPSG:4326'),
                                          pyproj.Proj(destination_crs_string),
                                          origin_lon,
                                          origin_lat)
    # See https://github.com/mapbox/rasterio/blob/master/examples/reproject.py
    dst_shape = (64, 64)  # Test
    dst_transform = transform.from_origin(origin_x, origin_y, resolution, resolution)
    dst_crs = crs.CRS.from_string(destination_crs_string)
    reprojected_mapbiomass_values = np.zeros(dst_shape, np.uint8)

    original_mapbiomass_raster = rasterio.open(infile_path)
    original_mapbiomass_values = original_mapbiomass_raster.read()  # Memory error
    reproject(original_mapbiomass_values,
              reprojected_mapbiomass_values,
              src_transform=original_mapbiomass_raster.transform,
              src_crs=original_mapbiomass_raster.crs,
              dst_transform=dst_transform,
              dst_crs=dst_crs,
              resampling=Resampling.nearest)
    pdb.set_trace()

    with rasterio.open(outfile_path,
                       'w',
                       driver='GTiff',
                       width=dst_shape[1],
                       height=dst_shape[0],
                       count=reprojected_mapbiomass_values.shape[0],
                       dtype=np.uint8,
                       nodata=None,
                       transform=dst_transform,
                       crs=dst_crs) as dst:
        dst.write(reprojected_mapbiomass_values)

if __name__ == '__main__':
    main()

---

Edit: I have

$ gdalwarp --version
GDAL 1.11.3, released 2015/09/16

and based on the suggestions in the comments I am trying

gdalwarp -tr 30.0 30.0 -s_srs EPSG:4326 -t_srs "+proj=aea +lat_1=10 +lat_2=-40 +lat_0=-25 +lon_0=-50 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs" AMAZONIA.tif AMAZONIA_reprojected.tif

Answer 1

The original tiff is written in one row high stripes as can be seen from gdalwarp report Block=123904x1. That is suboptimal structure if amount of memory is a bottleneck but with your 6 GB it should not be. If you want to make the source data more gdalwarp friendly you can make an interim version with gdal_translate:

gdal_translate -of GTiff -co tiled=yes -co compress=LZW AMAZONIA.tif AMAZONIA_tiled.tif

Gdalwarp command to use would then be like

gdalwarp -of GTiff -co TILED=YES -co COMPRESS=LZW -tr 30.0 30.0 -s_srs EPSG:4326 -t_srs "+proj=aea +lat_1=10 +lat_2=-40 +lat_0=-25 +lon_0=-50 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs" AMAZONIA_tiled.tif AMAZONIA_reprojected.tif

With these parameters the process is using one CPU core and on my computer it seems to take about 890 MB of memory. Conversion with gdalwarp takes 40 minutes with this laptop but I did not measure the time of the first gdal_translate run. Make some tests and add your timings, it would be interesting to see your results.

EDIT With untiled input gdalwarp takes a bit more memory but the difference is not huge: 930 MB vs. 890 MB. However, with 123904 wide stripes gdalwarp is extremely slow, I would say unusable. So I would say that preparing data with gdal_translate into tiled tiff is a compulsory step.

Answer 2

if you want to re-project, re-sample and stack a bunch of rasters together you can use the following gist that i wrote using rasterio:

https://gist.github.com/prakharcode/b83caaaa2fc6d2d62b7fe558656df0d1

the logic to do any of this for a single raster would remain the same.