Matching two rasters with different projections and resolution

Question

I have one elevation file (SRTM) and one satellite imagery file (sentinel-2), both are in GeoTIFF format. The elevation file has a much larger area covered. Here is the gdalinfo of the elavation file:

> Driver: GTiff/GeoTIFF Files: C:\DEV\3D\Assets\Heightmaps\nv\full.tif
> Size is 10801, 7201 Coordinate System is: GEOGCS["WGS 84",
>     DATUM["WGS_1984",
>         SPHEROID["WGS 84",6378137,298.257223563,
>             AUTHORITY["EPSG","7030"]],
>         AUTHORITY["EPSG","6326"]],
>     PRIMEM["Greenwich",0],
>     UNIT["degree",0.0174532925199433],
>     AUTHORITY["EPSG","4326"]] Origin = (-2.000138888888889,36.000138888888884) 
> Pixel Size =(0.000277777777778,-0.000277777777778) 
> Metadata:   AREA_OR_POINT=Area
> Image Structure Metadata:   INTERLEAVE=BAND Corner Coordinates: 
> Upper Left  (  -2.0001389,  36.0001389) (  2d 0' 0.50"W, 36d 0' 0.50"N)
> Lower Left  (  -2.0001389,  33.9998611) (  2d 0' 0.50"W,33d59'59.50"N) 
> Upper Right (   1.0001389,  36.0001389) (  1d 0'0.50"E, 36d 0' 0.50"N) 
> Lower Right (   1.0001389,  33.9998611) (  1d 0' 0.50"E,33d59'59.50"N) 
> Center(  -0.5000000,  35.0000000) ( 0d30' 0.00"W, 35d 0' 0.00"N) 
> Band 1 Block=10801x1 Type=Int16,ColorInterp=Gray

and the satellite image:

>   Driver: GTiff/GeoTIFF Files:
> C:\DEV\3D\Assets\Imagery\new\sub_big_merge.tif Size is 15017, 8131
> Coordinate System is: PROJCS["WGS 84 / UTM zone 30N",
>     GEOGCS["WGS 84",
>         DATUM["WGS_1984",
>             SPHEROID["WGS 84",6378137,298.257223563,
>                 AUTHORITY["EPSG","7030"]],
>             AUTHORITY["EPSG","6326"]],
>         PRIMEM["Greenwich",0,
>             AUTHORITY["EPSG","8901"]],
>         UNIT["degree",0.0174532925199433,
>             AUTHORITY["EPSG","9122"]],
>         AUTHORITY["EPSG","4326"]],
>     PROJECTION["Transverse_Mercator"],
>     PARAMETER["latitude_of_origin",0],
>     PARAMETER["central_meridian",-3],
>     PARAMETER["scale_factor",0.9996],
>     PARAMETER["false_easting",500000],
>     PARAMETER["false_northing",0],
>     UNIT["metre",1,
>         AUTHORITY["EPSG","9001"]],
>     AXIS["Easting",EAST],
>     AXIS["Northing",NORTH],
>     AUTHORITY["EPSG","32630"]] 
> Origin = (614090.000000000000000,3987110.000000000000000) 
> Pixel Size =(10.000000000000000,-10.000000000000000) 
> Metadata:  
> AREA_OR_POINT=Area Image Structure Metadata:   INTERLEAVE=PIXEL Corner
> Coordinates: 
> Upper Left  (  614090.000, 3987110.000) (1d44' 1.72"W,36d 1'18.55"N) 
> Lower Left  (  614090.000, 3905800.000) (1d44'43.19"W, 35d17'19.93"N) 
> Upper Right (  764260.000, 3987110.000) (0d 4' 6.98"W, 35d59'33.60"N) 
> Lower Right (  764260.000, 3905800.000) (0d 5'42.84"W, 35d15'37.76"N) 
> Center      (  689175.000, 3946455.000) (0d54'37.94"W, 35d38'37.76"N) 
> Band 1 Block=15017x1 Type=Byte,
> ColorInterp=Red Band 2 Block=15017x1 Type=Byte, ColorInterp=Green Band
> 3 Block=15017x1 Type=Byte, ColorInterp=Blue

Now, how do I get a subset of the elevation file that exactly matches the same area of the satellite image with gdal_translate or any other tool (since as you can see the two files have differents projections, resolutions, units ... etc)?

In the meantime I have tried these two "solutions"

Use GDAL command line to copy projections

Use gdalsrsinfo to get the srs of the tiff that still has the projection:

gdalsrsinfo -o wkt tiffwithsrs.tiff

Then copy the output and use gdal_translate to apply it to a new tiff:

gdal_translate -a_srs '...' tiffwithoutsrs.tif newfixedtif.tif

just paste your projection after the -a_srs option wich did nothing .

and this which completely broke my file.

Answer 1

First of all, I'd move the two raster layers in the same working directory in order to skip the long paths in command line.

We can transform the elevation layer to the same CRS of the satellite image, i.e. EPSG:32630:

gdalwarp -t_srs EPSG:32630 full.tif full_32630.tif

Then we can compute the common area between the two rasters:

gdaltindex elevation_tile.shp full_32630.tif
gdaltindex sat_image_tile.shp sub_big_merge.tif

We can write an OGR VRT layer called e.g. tiles.vrt like this:

<OGRVRTDataSource>
    <OGRVRTLayer name="elevation_tile">
        <SrcDataSource>elevation_tile.shp</SrcDataSource>
    </OGRVRTLayer>
    <OGRVRTLayer name="sat_image_tile">
        <SrcDataSource>sat_image_tile.shp</SrcDataSource>
    </OGRVRTLayer>
</OGRVRTDataSource>

and calculate the intersection:

ogr2ogr intersection.shp tiles.vrt -dialect sqlite -sql "SELECT ST_Intersection(elevation_tile.geometry, sat_image_tile.geometry) AS geometry FROM elevation_tile, sat_image_tile"

Finally, we can crop the two raster layers with intersection.shp and set the target resolution as 30 m (because 0.000277777777778° =~ 30m is the lowest resolution, it should be the target one):

gdalwarp -cutline intersection.shp -tr 30 30 full_32630.tif elevation_cropped_30m.tif
gdalwarp -cutline intersection.shp -tr 30 30 sub_big_merge.tif sat_image_cropped_30m.tif

Result: elevation_cropped_30m.tif and sat_image_cropped_30m.tif spatial extent and resolution should match as requested.

Answer 2

If you can use Python and its GDAL bindings, you can have a look at my reproject_image_to_raster method. Basically, it takes one image as the master (in this case, your S2), and it would reproject the slave to match extension and output resolution. I think you only need access to the gdal (e.g. you need an import gdal).

Answer 3

To reproject a raster to match another raster's extents and resolution I use python:

#!/usr/bin/env python
# https://github.com/drf5n/drf5n-public/blob/master/gdalwarp2match.py
from osgeo import gdal, gdalconst
import argparse
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)")
args = parser.parse_args()
print(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
GRA_Bilinear Bilinear give diamond-shaped artifacts
GRA_CubicSpline Cubic-Spline smooths them
GRA_NearestNeighbour nearest neighbor????
interpolation=gdalconst.GRA_NearestNeighbour
#gdal.ReprojectImage(src, dst, src_proj, match_proj, gdalconst.GRA_Bilinear)
#gdal.ReprojectImage(src, dst, src_proj, match_proj, gdalconst.GRA_CubicSpline)
gdal.ReprojectImage(src, dst, src_proj, match_proj, interpolation)
del dst # Flush

Answer 4

The following python function will do three things (you explicitly asked for 2, but it might be helpful for others). The function will:

• take a leader and a follower raster datasets and pass:
  • the crs
  • the domain
  • the resolution

of the leader dataset to the follower dataset

the follower_dataset and the leader_dataset should be rio.open(path to .tifs).

The out_name is the path to your output .tif (follower dataset with the crs, domain and resolution of the leader dataset).

import rasterio as rio
from sklearn.impute import SimpleImputer
def reshape_raster(follower_dataset,leader_dataset, out_name):
dst_crs = follower_dataset.crs

out_width = leader_dataset.shape[1]
out_height = leader_dataset.shape[0]
imputer = SimpleImputer(fill_value=np.nan, strategy='mean')

with follower_dataset as src:
    profile = src.profile.copy()
    transform, width, height = calculate_default_transform(
            src.crs, 
            dst_crs, 
            src.width, 
            src.height,
            left = leader_dataset.bounds.left,
            bottom = leader_dataset.bounds.bottom,
            right = leader_dataset.bounds.right,
            top = leader_dataset.bounds.top,
            dst_width=out_width, 
            dst_height=out_height)

    profile.update({
            'crs': dst_crs,
            'transform': transform,
            'width': width,
            'height': height,
            'compress':'lzw'
            })

    with rio.open(out_name, 'w', **profile) as dst:
        data = src.read(1)
        data = imputer.fit_transform(data)
        reproject(
                source=rio.band(src, 1),
                destination=rio.band(dst, 1),
                src_transform=src.transform,
                src_crs = dst_crs,
                dst_transform=transform,
                dst_crs=dst_crs,
                resampling=Resampling.nearest
                )
    dst.FlushCache()
    dst.close()

    out=rio.open(out_name)
return out