A recurrent problem, when analysing GEOS-Chem outputs and/or comparing it 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 don't provide the altitude a.s.l of the base levels (i.e., the model topography).
This notebooks shows how to load global topographic elevation data (any raster format supported by the GDAL library) and then resample it on a given GEOS-Chem grid. The original, high resolution elevation data is downsampled to a GC grid using the regridding methods available in the Iris package.
NOTE: Both data loading and regridding operations performed here are not optimized for very high resolution datasets in terms of memory consumption and computing efficiency!
Author: B. Bovy (Nov. 2014)
In [1]:
import numpy as np
import gdal, osr
import iris
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import iris.plot as iplt
import pygchem
plt.style.use('ggplot')
%matplotlib inline
In [2]:
iris.__version__
Out[2]:
In [3]:
pygchem.__version__
Out[3]:
The example below uses the USGS/NGA's GMTED2010 global elevation data (Mean Statistic, 30 arc-seconds version). It is a high quality, multi-resolution product for global and continental applications. More info here: http://topotools.cr.usgs.gov/gmted_viewer/. The dataset can be downloaded in the ESRI's ArcGIS/Info binary grid format.
NOTE: The elevation raster must have lat/lon coordinates (no projection)! It is also assumed that the geographic coordinate system of the elevation raster is equivalent to the coordinate system used by the GEOS-Chem grids (GMAO Met Fields). No datum/ellipsoid/spheroid transformation is performed here.
In [4]:
dem_file = "/home/bovy/GISData/Elevation/GMTED2010/mn30_grd"
dem_name = "GMTED2010_mn30"
dem_units = "m"
In [5]:
!gdalinfo $dem_file -norat -noct -nogcp
NOTE: Since all data is loaded into RAM, it can consume a lot of memory when using high resolution data!
In [6]:
dem_rst = gdal.Open(dem_file)
dem_data = dem_rst.ReadAsArray().astype('d')
In [7]:
dem_cs = osr.SpatialReference()
dem_cs.ImportFromWkt(dem_rst.GetProjection())
if not dem_cs.IsGeographic():
raise ValueError("The raster's coordinate system must be geographic!")
GeogCS object from the raster's coord system.
In [8]:
cube_cs = iris.coord_systems.GeogCS(
semi_major_axis=dem_cs.GetSemiMajor(),
semi_minor_axis=dem_cs.GetSemiMinor(),
longitude_of_prime_meridian=dem_cs.GetProjParm('PRIMEM')
)
print(cube_cs)
In [9]:
dem_gt = dem_rst.GetGeoTransform()
nrows, ncols = dem_data.shape
y_ind = np.arange(0, nrows + 1)
x_ind = np.arange(0, ncols + 1)
dem_lon_left = dem_gt[0] + (x_ind * dem_gt[1])
dem_lat_top = dem_gt[3] + (y_ind * dem_gt[5])
dem_lon_bounds = np.column_stack((dem_lon_left[:-1], dem_lon_left[1:]))
dem_lat_bounds = np.column_stack((dem_lat_top[:-1], dem_lat_top[1:]))
dem_lon_centers = (dem_lon_left[0:-1] + dem_lon_left[1:]) / 2.
dem_lat_centers = (dem_lat_top[0:-1] + dem_lat_top[1:]) / 2.
cube for the global topography (elevation + coordinates)
In [10]:
clon = iris.coords.DimCoord(
dem_lon_centers,
bounds=dem_lon_bounds,
standard_name='longitude',
units='degrees_east',
coord_system=cube_cs
)
clat = iris.coords.DimCoord(
dem_lat_centers,
bounds=dem_lat_bounds,
standard_name='latitude',
units='degrees_north',
coord_system=cube_cs
)
dem_cube = iris.cube.Cube(
dem_data,
var_name=dem_name,
units=dem_units,
dim_coords_and_dims=((clat, 0), (clon, 1))
)
print(dem_cube)
In [11]:
lon_extent = (-12, 30)
lat_extent = (35, 65)
lat_subset = iris.Constraint(
latitude=lambda cell: lat_extent[0] <= cell <= lat_extent[1]
)
lon_subset = iris.Constraint(
longitude=lambda cell: lon_extent[0] <= cell <= lon_extent[1]
)
dem_subset = dem_cube.extract(lat_subset & lon_subset)
plt.figure(figsize=(10, 6))
ax = plt.axes(projection=ccrs.PlateCarree())
topo_map = ax.imshow(dem_subset.data,
origin='upper',
extent=lon_extent + lat_extent,
transform=ccrs.PlateCarree(),
interpolation='nearest',
cmap=plt.cm.terrain,
vmin=-1000, vmax=4500)
ax.coastlines(resolution='50m')
plt.colorbar(topo_map,
ax=ax, orientation='horizontal',
pad=0.05, shrink=0.65,
label='m')
plt.title(dem_cube.name())
plt.show()
The pygchem.grid module can be used to get one of the standard grids used in GEOS-Chem simulations.
In [12]:
import pygchem.grid as grid
Grid setup is available for the following grids:
In [13]:
grid.get_supported_models()
Out[13]:
In [14]:
gcgrid_model = 'GEOSFP'
gcgrid_resolution = 5, 4 # horizontal grid resolution (lon, lat) in decimal degrees
gcgrid_name = "{model}_{rlat}x{rlon}".format(
model=gcgrid_model,
rlat=gcgrid_resolution[1],
rlon=gcgrid_resolution[0]
)
gcgrid_name
Out[14]:
CTMGrid object from the grid name and resolution given here above
In [15]:
gcgrid = grid.grid_from_model(gcgrid_model,
resolution=gcgrid_resolution)
DimCoord objects) from the grid defined here above. The pygchem.tools.irisutils module provides handful routines for creating or handling such Iris objects.
In [16]:
import pygchem.tools.irisutil as irisutil
gclon, gclat = irisutil.coord_from_grid(
gcgrid, coord=('longitude', 'latitude')
)
In [17]:
gclon.coord_system = cube_cs
gclat.coord_system = cube_cs
cube for the defined grid
In [18]:
gc_cube = iris.cube.Cube(
np.empty(gclat.shape + gclon.shape),
var_name=gcgrid_name,
units=dem_units,
dim_coords_and_dims=((gclat, 0), (gclon, 1))
)
print(gc_cube)
See the Iris user guide about cube regridding and interpolation for more info. Here the area weighted conservative scheme is used to properly resample the high-resolution elevation data down to the low-resolution GEOS-Chem grid.
NOTE: the execution of the cell below may take a while, depending on the size of the elevation data (it may take > 1 hour for the elevation data used here!). This regridding implementation is not optimized for processing huge sized arrays!
In [19]:
regrid_scheme = iris.analysis.AreaWeighted(mdtol=0.5)
dem_regrid_cube = dem_cube.regrid(gc_cube, regrid_scheme)
In [20]:
dem_regrid_cube.var_name = "TOPO_{}".format(gc_cube.name())
dem_regrid_cube.standard_name = "height"
dem_regrid_cube.attributes['topo_data_src'] = dem_cube.name()
dem_regrid_cube.attributes['regrid_scheme'] = "area-weighted mean"
dem_regrid_cube.add_cell_method(
iris.coords.CellMethod("mean")
)
print(dem_regrid_cube)
In [21]:
iris.save(
dem_regrid_cube,
'/home/bovy/Grids/topo_{gc_name}.nc'.format(gc_name=gc_cube.name())
)
In [22]:
plt.figure(figsize=(9, 5))
ax = plt.axes(projection=ccrs.PlateCarree())
iplt.pcolormesh(dem_regrid_cube,
cmap=plt.cm.terrain,
vmin=-1000, vmax=4500)
ax.figure.set_facecolor('white')
ax.coastlines()
ax.set_global()
plt.colorbar(ax=ax, orientation='horizontal',
pad=0.05, shrink=0.7,
label='m')
plt.title(dem_regrid_cube.var_name)
plt.show()
In [ ]: