I have two NetCDF4 files with different grids and dimensions (but both approximately cover the same area, see images below). One has information on land-use, the other on topography.
I now want to add this topography variable to the land-use file by assigning each elevation value in the topography file to the correct gridcell in the land-use file based on the nearest lonlat coordinate.
How would I do this in Python?
If I understood your question correctly, you want to simply resample the data from one grid to another.
I recommend using the python package pyresample.
Here's the code snippet:
from netCDF4 import Dataset
import numpy as np
from pyresample import geometry, kd_tree
LU_file = Dataset('LUindexHCLIM3.nc', mode='a')
topo_file = Dataset('topography.nc', mode='r')
# Information from file_LU
nature_cover = np.array(LU_file.variables['Main_Nature_Cover'])
lon_LU = np.array(LU_file.variables['lon'])
lat_LU = np.array(LU_file.variables['lat'])
# Information from file_topo
height = np.array(topo_file.variables['HGT_M'])
# Use pyresample for interpolation
geometry_topo = geometry.SwathDefinition(lons=lon_topo, lats=lat_topo)
geometry_LU = geometry.SwathDefinition(lons=lon_LU, lats=lat_LU)
interpolated_height = kd_tree.resample_nearest(geometry_topo, height, geometry_LU, radius_of_influence=500000, fill_value=None)
interpolated_height = interpolated_height.filled(np.nan)
# Add the interpolated data to the first netCDF file
height_LU = LU_file.createVariable('HGT_M', 'f4', ('y', 'x'))
height_LU[:] = interpolated_height
LU_file.close()
topo_file.close()
Of course, you can play around with the different interpolation schemes to get a better result.
