Source code for pyjams.gridcellarea

#!/usr/bin/env python
"""
Area of grid cells on Earth

This module was written by Matthias Cuntz while at Department of
Computational Hydrosystems, Helmholtz Centre for Environmental
Research - UFZ, Leipzig, Germany, and continued while at Institut
National de Recherche pour l'Agriculture, l'Alimentation et
l'Environnement (INRAE), Nancy, France.

:copyright: Copyright 2009-2022 Matthias Cuntz, see AUTHORS.rst for details.
:license: MIT License, see LICENSE for details.

.. moduleauthor:: Matthias Cuntz

The following functions are provided

.. autosummary::
   gridcellarea

History
    * Written Jul 2009 by Matthias Cuntz - mc (at) macu (dot) de
    * Ported to Python 3, Feb 2013, Matthias Cuntz
    * Assert correct input, Apr 2014, Matthias Cuntz
    * Clip lat because -90 and 90 give negative area, Oct 2018, Matthias Cuntz
    * Corrected bug in checking ascending or descending lats,
      Apr 2022, Matthias Cuntz
    * Removed assertion of nlat/nlon > 2, Apr 2022, Matthias Cuntz
    * Keyword radius of Earth and changed default from 6371000 to 6371009
      as in the rest of pyjams, Apr 2022, Matthias Cuntz
    * Rename cellarea to gridcellarea, Apr 2022, Matthias Cuntz

"""
import numpy as np


__all__ = ['gridcellarea']


[docs]def gridcellarea(lat, lon, globe=False, rearth=6371009.): """ Area of grid cells on a spherical Earth in square metre Parameters ---------- lat : array_like Latitudes of grid cell centres in degrees N lon : array_like Longitudes of grid cell centres in degrees E globe : bool, optional Assumes that latitudes span the globe if True, i.e. they are bounded by 90 and -90 degrees latitude rearth : float, optional Radius of the spherical Earth (default: 6371009.) Returns ------- numpy array with area in m^2 Notes ----- This is a rather rough routine with lots of possible improvements, notably allowing irregular grids, grid boundaries, etc. Examples -------- Gaussian latitudes >>> lat = np.array([ 12.98898858, 9.27785325, 5.56671363]) >>> lon = np.array([ 0., 3.75, 7.5]) >>> print(gridcellarea(lat, lon)[0,:]) [1.67639557e+11 1.67639557e+11 1.67639557e+11] >>> print(gridcellarea(lat, lon)[1,:]) [1.69790907e+11 1.69790907e+11 1.69790907e+11] >>> print(gridcellarea(lat, lon)[2,:]) [1.71230373e+11 1.71230373e+11 1.71230373e+11] """ # assert numpy lati = np.array(lat) loni = np.array(lon) nlat = lati.size nlon = loni.size # assert -90 < lat < 90 assert np.abs(lat).max() <= 90., ( 'probably swapped lat and lon in call:' ' def gridcellarea(lat, lon, globe=False):') # lat size in degrees # still + or - dlat = lati - np.roll(lati, 1) if globe: if dlat[0] > 0: # descending lats l0 = 90. else: # ascending lats l0 = -90. dlat[0] = (lati[0] - l0) + 0.5 * (lati[1] - lati[0]) dlat[-1] = (-l0 - lati[-1]) - 0.5*(lati[-2] - lati[-1]) else: dlat[0] = lati[1] - lati[0] dlat[-1] = lati[-1] - lati[-2] # lon size in degrees # check if meridian in lon range -> shift to -180,180 # e.g. 358 359 0 1 if np.any(np.abs(np.diff(loni)) > 360. / nlon): loni = np.where(loni > 180., loni - 360., loni) # check if -180,180 longitude in lon range -> shift to 0,360 # e.g. 179 180 -179 -178 if np.any(np.abs(np.diff(loni)) > 360. / nlon): loni = np.where(loni < 0., loni + 360., loni) dlon = np.abs(loni - np.roll(loni, 1)) dlon[0] = np.abs(loni[1] - loni[0]) # Northern latitudes of grid cell edges n_lat = lati[0] + dlat[0]/2. + np.cumsum(dlat) # Area of grid cells in m^2 with lat/lon in degree d2r = np.pi / 180. # degree to radian area = np.empty([nlat, nlon]) dlat = np.abs(dlat[:]) # -90, 90 give negative np.cos(lat*d2r) lati = np.clip(lati, -89.99999, 89.99999) for i in range(nlon): area[:, i] = (2. * d2r * rearth**2 * dlon[i] * np.sin(0.5 * dlat * d2r) * np.cos(lati * d2r)) return area
if __name__ == '__main__': import doctest doctest.testmod(optionflags=doctest.NORMALIZE_WHITESPACE)