Search code examples
pythongisrastergdalosgeo

How to reproject and area match a netCDF and shapefile

I am new to GDAL and just getting my feet wet.

I'm trying to compare rasters stored in netCDFs to shapefiles. The shapefiles are sub-sections of the area covered by the netCDFs and the datasets use slightly different projections. I convert the shapefile dataset to the netCDF projection. The netCDF file contains the raster array and 1d arrays for lat, lon, x, y.

Right now, my code uses gdal.RasterizeLayer to rasterize a shapefile into a tiff, and then gdal.ReprojectImage to reproject that into a new tiff. My issue is that I cannot figure out how to determine the extents of the second tiff - which I need to select the subsection of the netCDF data.

Here is the relevant sections of my code:

#Extract projection information
obs_driver = ogr.GetDriverByName('ESRI Shapefile')  
obs_dataset = obs_driver.Open(obsfiles[0])
layer = obs_dataset.GetLayer()
obs_proj = layer.GetSpatialRef()
mod_proj = obs_proj.SetProjParm('parameter',90)   #Hardcode param difference
xmin,xmax,ymin,ymax = layer.GetExtent()   #Extents pre-reproject

Rasterizing 

tiff1 = gdal.GetDriverByName('GTiff').Create('temp1.tif', 1000, 1000, 1, gdal.GDT_Byte)
tiff1.SetGeoTransform((xmin, pixelWidth, 0, ymin, 0, pixelHeight))
gdal.RasterizeLayer(tiff1,[1],layer,options = ['ATTRIBUTE=attribute'])

Re-projecting

dst1 = gdal.GetDriverByName('GTiff').Create('temp3.tif', 1000, 1000, 1, gdal.GDT_Byte)
dst1.SetGeoTransform((xmin, pixelWidth, 0, ymin, 0, pixelHeight))
dst1.SetProjection(mod_proj.ExportToWkt())
gdal.ReprojectImage(gdal.Open('./temp1.tif'), dst1, obs_proj.ExportToWkt(), mod_proj.ExportToWkt(), gdalconst.GRA_Bilinear)

And convert raster to array for point-by-point comparison

import matplotlib.pyplot as plt
obs = plt.imread('temp3.tif')

So now I need to get the extents of this array (in the new projection) so I can match it up with the right subsection of the netCDF array, and interpolate it to match.

EDIT:Now I am thinking that I need to individually transform the extents and use that to redefine the GeoTransform for the projection conversion. Looking into it.

Solution

  • Figured it out - you need to transform the extents to the new projection type before converting the data.

    #Calculate new x,y positions of extents
tx = osr.CoordinateTransformation(obs_proj,mod_proj)
  #Projection corrected extent of area in question
corr_ext = [0,0,0,0]  #[min,max,ymin,ymax]
(corr_ext[0],corr_ext[2],ignore_z) = tx.TransformPoint(ext[0],ext[2]) #Ignore z componenet - data is 2D
(corr_ext[1],corr_ext[3],ignore_z) = tx.TransformPoint(ext[1],ext[3])

    And then use those extents to set the GeoTransform of the object containing the reprojected data. (The GeoTransform holds the spatial location information for the raster).