2

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
  • 2
    Is there a reason why you want to do this with rasterio, instead of GDALwarp? Curiosity and learning are good things but considering the size of the raster and your moderate (at best) hardware a compiled tool might be a good idea and leave the learning experience for a more easily accomplished task. – Michael Stimson Apr 5 '18 at 4:32
  • @Michael I'd be happy to do it with GDALwarp, I just don't know how. Would that require less memory? – Adrian Apr 5 '18 at 4:38
  • I do have access to a university cluster, so I could ask for e.g. 16 gigabytes of RAM for a few hours if my personal laptop is too puny. – Adrian Apr 5 '18 at 4:40
  • 1
    Oh, that would work much better, does it have a RAID? GDALwarp gdal.org/gdalwarp.html can be called from QGIS or on the command line (provided you're not afraid of CMD) use parameters -s_srs EPSG:4326 and -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" (your PROJ.4 string with quotes). Choose a good resampling method, the default is nearest which isn't so good for continuous rasters but works fine for classified/monochrome images. – Michael Stimson Apr 5 '18 at 4:49
  • 1
    Reprojecting a 2GB raster with gdalwarp is trivial, even if you don't mess around with multiple processes. However, if your raster is compressed and you want compressed output, gdalwarp to VRT first then gdal_translate to your final output with your -co compression= options - gis.stackexchange.com/a/89549/2856 – user2856 Apr 5 '18 at 5:44
3

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.

  • For compressed output, I recommend to gdalwarp to VRT first then gdal_translate to the final output with the -co compression= options - gis.stackexchange.com/a/89549/2856 – user2856 Apr 5 '18 at 8:19
  • @Luke, I fear that it does not work in this case because VRT refers to original tiff file that has unsuitable structure with 123904 pixels wide stripes. See my edit. Have you tried to make the conversion with the provided data? – user30184 Apr 5 '18 at 9:00
  • Thank you, I am trying this now. If it works I'll let you know and mark your answer as accepted. – Adrian Apr 6 '18 at 2:30
0

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.