Comparison between two images with different projections using GDAL

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?

• As per the Tour, we have a One question per Question policy. Please Edit your Question to focus on one question. Commented Jun 13, 2023 at 18:53
• Where do the numpy arrays come from? Asking because you can perhaps use the sources of the arrays directly in gdal_warp, instead of opening as arrays, converting to "gdal_objects" and then warping. Commented Jun 14, 2023 at 8:08
• The MODIS L2 will be from an HDF file while the simulation will be from a NetCDF4 file. This method is the one I am most familiar with. Commented Jun 14, 2023 at 9:20

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.

• I am not familiar with GeoPandas and Rasterio. But I will look into this. Commented Jun 14, 2023 at 9:37
• I did some reading on GeoPandas. Am I correct in thinking that I should convert the raster files first into shapefiles (or polygons) to determine the overlapping area? Commented Jun 15, 2023 at 9:59
• You can derive the bounds of the raster from a rasterio.dataset and convert it into a Polygon (shapely) then in a GeoDataFrame. Pay attention that if you want to have the footprint and not the extent of the raster (i.e. if you have nodata), you should use a vectorization technique of your raster Commented Jun 15, 2023 at 10:03
• I did some tests using your suggestions and basically did the following: (1) My first hurdle is that both the MODIS data and the simulation data that I have have no projection and coordinate reference systems, and simply converting them to GeoTIFF does not result into the correct projection. So what I did was "re-project" both of them into the same projection and CRS using `pyresample`. (2) From there, I open them both with Rasterio and did what was suggested here gis.stackexchange.com/a/447159/208986, which results into two Numpy arrays with the same size. Commented Jun 29, 2023 at 7:44
• My only concern is then I technically resampled both of my data twice, one in the `pyresample` and once again in the Rasterio. I wonder if there is some way that I can work with the original MODIS and simulation data, i.e. the data without the projection and CRS. Commented Jun 29, 2023 at 7:46