Comparison between two images with different projections using GDAL

Question

I am pretty new to GDAL, and I want to compare two data sets with different projections, specifically a MODIS L2 product and a simulation--both having their longitude/latitude and environmental measurement/simulation as Numpy arrays. These two images do not completely overlap, rather they intersect.

My plan to compare them is to project both data into the same projection and then just compare them pixel-by-pixel, so I would need to generate a GDAL object for both of them, since they're both Numpy arrays, sharing the same geotransform and projection.

I compute the geotransform for both data as such:

nx, ny = data.shape[1], data.shape[0]
x, y = (lon.max() - lon.min()) / float(nx), (lat.max() - lat.min()) / float(ny)
gt = (lon.min(), xRes, 0, lat.max(), 0, -yRes)

And then generate a GDAL object with this transform:

gdal_object = gdal_driver.Create('gdalObj', nx, ny, 1, gdal.GDT_Float32)
gdal_object.SetGeoTransform(gt)

The projection of the GDAL object will then be based on:

srs = osr.SpatialReference()
srs.ImportfromEPSG(4326)
gdal_object.SetProjection(srs.ExportToWkt())
gdal_object.GetRasterBand(1).WriteArray(data)
gdal_object.FlushCache()

I then resample both these GDAL objects onto the same resolution:

resampled_data = gdal.Warp('gdalObject', gdal_object, xRes, yRes, resampleAlg = 'bilinear')

My question is then, is this approach on comparing two different data sets with different projections correct or efficient?

Answer 1

I would advise you to use more up-to-date and pythonic libraries such as rasterio or rioxarray to simplify your code instead of using the GDAL bindings.

What you should do is:

• determine the overlapping area, i.e. with geopandas
• convert the area into a rasterio Window for each raster (from lat/lon to x/y)
• open the two rasters using their respective window
• reproject one raster onto the other to align their grids
• then you can compare pixel per pixel

Once you've aligned your two rasters on the same grid, it is completely correct to compare your rasters. Just pay attention to your resampling method and what you try to achieve, i.e. nearest will keep you values scientifically acceptable when bilinear will smooth your results.

However, since you are simulating something, I would advise you to simulate directly on MODIS grid, it will simplify the process.