A recurrent problem, when analysing GEOS-Chem outputs and/or comparing them with other data sources, is to get the altitude values at the vertical levels of a GEOS-Chem grid. GEOS-Chem outputs include grid-box heights (in meters) but doesn't provide the altitude (a.s.l) of the base levels (i.e., the model topography).
The code below allows to load global elevation data (any raster format supported by the GDAL library) and then resample it to a given GEOS-Chem grid. The original, high resolution elevation data is downsampled to a GC grid using several regridding methods.
The example below uses the USGS's and NGA's GMTED2010 global elevation data. It is a high quality, multi-resolution product for global and continental applications. More info here: http://topotools.cr.usgs.gov/gmted_viewer/.
WARNING: The elevation raster must have lat/lon coordinates (no projection) ! It is also assumed that the geographic coordinate system of the elevation raster is roughly similar to the coordinate system used by the GEOS-Chem grids (GMAO Met Fields). No datum/ellipsoid/spheroid transformation is performed in the example below.
In [16]:
import numpy as np
import gdal, osr
from scipy.io import readsav
import iris
from iris.coords import DimCoord
from iris.experimental.regrid import regrid_area_weighted_rectilinear_src_and_grid
from iris.experimental.regrid_conservative import regrid_conservative_via_esmpy
import iris.quickplot as qiplot
import ESMF
import pygchem.grid as gcgrid
In [2]:
# GMTED2010 product at 30 arc-seconds, aggregation method = mean elevation, ArcGRID format.
elev_src = "/home/bovy/GISData/Elevation/GMTED2010/mn30_grd"
elev_name = "GMTED2010_mn30"
elev_units = "m"
gcgrid_model = 'GEOS57' # GEOS-Chem grid model (see the cell below for model names)
# Note: many models use the same horizontal grid
gcgrid_resolution = 2.5, 2 # horizontal grid resolution (longitude, latitude) in decimal degrees
output_dir = '/home/bovy/Grids/' # directory where regridded elevation data files will be saved
In [3]:
gcgrid.CTMGrid.models
Out[3]:
In [4]:
elev_rst = gdal.Open(elev_src)
elev_data = elev_rst.ReadAsArray()
In [5]:
elev_cs = osr.SpatialReference()
elev_cs.ImportFromWkt(elev_rst.GetProjection())
if not elev_cs.IsGeographic():
raise ValueError("The raster's coordinate system must be geographic!")
# create an iris's GeogCS object from the raster's coord system
cube_cs = iris.coord_systems.GeogCS(
semi_major_axis=elev_cs.GetSemiMajor(),
semi_minor_axis=elev_cs.GetSemiMinor(),
longitude_of_prime_meridian=elev_cs.GetProjParm('PRIMEM')
)
print elev_cs
In [6]:
elev_gt = elev_rst.GetGeoTransform()
cols = elev_rst.RasterXSize
rows = elev_rst.RasterYSize
y_ind = np.arange(0, rows + 1)
x_ind = np.arange(0, cols + 1)
elev_lon_left = elev_gt[0] + (x_ind * elev_gt[1])
elev_lat_top = elev_gt[3] + (y_ind * elev_gt[5])
elev_lon_bounds = np.column_stack((elev_lon_left[0:-1], elev_lon_left[1:]))
elev_lat_bounds = np.column_stack((elev_lat_top[0:-1], elev_lat_top[1:]))
elev_lon_centers = (elev_lon_left[0:-1] + elev_lon_left[1:]) / 2.
elev_lat_centers = (elev_lat_top[0:-1] + elev_lat_top[1:]) / 2.
In [7]:
clon = DimCoord(elev_lon_centers,
bounds=elev_lon_bounds,
standard_name='longitude',
units='degrees_east',
coord_system=cube_cs)
clat = DimCoord(elev_lat_centers,
bounds=elev_lat_bounds,
standard_name='latitude',
units='degrees_north',
coord_system=cube_cs)
elev_cube = iris.cube.Cube(elev_data.astype('d'),
var_name=elev_name,
units=elev_units,
dim_coords_and_dims=((clat, 0), (clon, 1)))
print elev_cube
In [8]:
gc = gcgrid.CTMGrid.from_model(gcgrid_model, resolution=gcgrid_resolution)
gc_name = "{model}_{rlat}x{rlon}".format(model=gcgrid_model,
rlat=gcgrid_resolution[1],
rlon=gcgrid_resolution[0])
gc_lon_centers, gc_lat_centers = gc.lonlat_centers
elon, elat = gc.lonlat_edges
gc_lon_bounds = np.column_stack((elon[0:-1], elon[1:]))
gc_lat_bounds = np.column_stack((elat[0:-1], elat[1:]))
In [9]:
clon = DimCoord(gc_lon_centers,
bounds=gc_lon_bounds,
standard_name='longitude',
units='degrees_east',
coord_system=cube_cs)
clat = DimCoord(gc_lat_centers,
bounds=gc_lat_bounds,
standard_name='latitude',
units='degrees_north',
coord_system=cube_cs)
dummy_data = np.empty((gc_lat_centers.size, gc_lon_centers.size))
gc_cube = iris.cube.Cube(dummy_data,
var_name=gc_name,
units=elev_units,
dim_coords_and_dims=((clat, 0), (clon, 1)))
print gc_cube
In [13]:
gc_elev_cube_awm = regrid_area_weighted_rectilinear_src_and_grid(elev_cube, gc_cube)
gc_elev_cube_awm.var_name = "DEM_{}".format(gc_cube.name())
gc_elev_cube_awm.attributes['elevation_data_src'] = elev_cube.name()
gc_elev_cube_awm.attributes['regrid_method'] = "area-weighted mean"
print gc_elev_cube_awm
In [39]:
gc_elev_cube_conservative = regrid_conservative_via_esmpy(elev_cube, gc_cube)
gc_elev_cube_conservative.var_name = "DEM_{}".format(gc_cube.name())
gc_elev_cube_conservative.attributes['elevation_data_src'] = elev_cube.name()
gc_elev_cube_conservative.attributes['regrid_method'] = "conservative"
print gc_elev_cube_conservative
In [10]:
gc_elev_cube_bilinear = iris.analysis.interpolate.regrid(elev_cube, gc_cube, mode='bilinear')
gc_elev_cube_bilinear.var_name = "DEM_{}".format(gc_cube.name())
gc_elev_cube_bilinear.attributes['elevation_data_src'] = elev_cube.name()
gc_elev_cube_bilinear.attributes['regrid_method'] = "bilinear"
print gc_elev_cube_bilinear
In [14]:
iris.save(gc_elev_cube_awm, '/home/bovy/temp/dem_{gc_name}_awm.nc'.format(gc_name=gc_cube.name()))
iris.save(gc_elev_cube_bilinear, '/home/bovy/temp/dem_{gc_name}_bilinear.nc'.format(gc_name=gc_cube.name()))
In [11]:
lat_selection = iris.Constraint(latitude=lambda cell: 40. <= cell <= 55.)
lon_selection = iris.Constraint(longitude=lambda cell: 0. <= cell <= 20.)
elev_cube_subset = elev_cube.extract(lat_selection & lon_selection)
qiplot.pcolormesh(elev_cube_subset, vmax=3000, cmap=cm.gist_earth)
plt.gca().coastlines()
plt.show()
In [15]:
qiplot.pcolormesh(gc_elev_cube_awm, vmax=4000, cmap=cm.gist_earth)
plt.gca().coastlines()
plt.show()
In [12]:
qiplot.pcolormesh(gc_elev_cube_bilinear, vmax=4000, cmap=cm.gist_earth)
plt.gca().coastlines()
plt.show()
In [ ]: