I have several sets of rasters which represent a large portion of the world. Each set is defined from the same grid cell, but rasters in this set come from different sources, therefore these vary in size between them. I need to normalize each raster, average each set, and build maptiles with all grid cells. So far I reckon I could accomplish this following these steps:

  1. For each set of rasters, use coordinates from pixels from one of them as a "base grid" to sample all rasters in this set.
  2. Normalize the sampled data.
  3. Build a geodataframe with the averaged data of each set of rasters.
  4. Save said geodataframe as a raster image.
  5. Repeat with each set of rasters.
  6. Build a VRT with all rasters with averaged data.
  7. Use gdal to build maptiles from the VRT.

My biggest concern is how slow would it be to get coordinates for each pixel and sample on all rasters for these coordinates. Pixel by pixel on each raster.

Is there a faster way to merge different rasters into one while doing normalization and averaging the result?

Is this approach correct? How would you tackle it?


1 Answer 1


You do need to get the data onto a common grid. I think I'd try to do it all with gdal command line tools.

For step 1, I use my own gdlwarp2match.py program:

#!/usr/bin/env python
# https://github.com/drf5n/drf5n-public/blob/master/gdalwarp2match.py

from osgeo import gdal, gdalconst
import argparse

# some mappings per https://gdal.org/programs/gdalwarp.html and https://gdal.org/python/osgeo.gdalconst-module.html
resampling = { 'near': gdalconst.GRA_NearestNeighbour,
                   'bilinear': gdalconst.GRA_Bilinear,
                   'cubic': gdalconst.GRA_Cubic,}

parser = argparse.ArgumentParser(description='Use GDAL to reproject a raster to match the extents and res of a template')
parser.add_argument("source", help="Source file")
parser.add_argument("template", help = "template with extents and resolution to match")
parser.add_argument("destination", help = "destination file (geoTIFF)")
parser.add_argument("--resample", choices=resampling.keys(),
                    help="""Resampling/interpolation method """, default="near")

args = parser.parse_args()

# Source
src_filename = args.source
src = gdal.Open(src_filename, gdalconst.GA_ReadOnly)
src_proj = src.GetProjection()
src_geotrans = src.GetGeoTransform()

# We want a section of source that matches this:
match_filename = args.template
match_ds = gdal.Open(match_filename, gdalconst.GA_ReadOnly)
match_proj = match_ds.GetProjection()
match_geotrans = match_ds.GetGeoTransform()
wide = match_ds.RasterXSize
high = match_ds.RasterYSize

# Output / destination
dst_filename = args.destination
dst = gdal.GetDriverByName('GTiff').Create(dst_filename, wide, high, 1, gdalconst.GDT_Float32)
dst.SetGeoTransform( match_geotrans )
dst.SetProjection( match_proj)

# Do the work
gdal.ReprojectImage(src, dst, src_proj, match_proj, resampling[args.resample])

del dst # Flush

Then for step 3 I'd probably use gdal_calc.py to do the averaging, and for step 6&7 I'd use gdal_merge.py to assemble them.

But this procedure doesn't do the normalizing. If normalizing is needed, I'd look into https://jgomezdans.github.io/gdal_notes/reprojection.html and modify my script to do the normalization in-memory.

I'd also worry that gdal_merge.py wouldn't handle possible overlapping as desired, but I'm also not sure that VRT would handle overlaps any differently. The normalization step seems like it could cause discontinuities at the mosaic'd edges, but it really depends on the datasets and where the edges are.

  • I can't run your script: it fails with "ImportError: No module named '_gdal'" . I've tried several solutions found here to no avail. Nov 5, 2021 at 18:16
  • Nevermaind that: For some reason it didn't work inside a conda environment, but it does work from the command line running python3 Nov 5, 2021 at 18:21
  • Maybe the shebang line is choosing the wrong python in your environment. You can edit it to use the python in your preferred conda environment directly, and then the script wouldn't depend on activating that environment in its parent shell.
    – Dave X
    Nov 6, 2021 at 11:53
  • 1
    I ended up incorporating your script into mine with good results. After that normalization is accomplished with gdal_translate, and averaging with gdal_calc as you suggested. Maptiles are generated with gdal2tiles from a VRT made of all averaged rasters. Thanks! Nov 9, 2021 at 14:31
  • If github.com/OSGeo/gdal/blob/master/swig/python/gdal-utils/… didn't fault on differing size, extents, or projection for the different layers, and instead warped to match, I wouldn't need my script. But it for my own use case it was easier for me to preprocess the data. I was surprised at how easily gdal made it to make it general.
    – Dave X
    Nov 9, 2021 at 14:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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