Source code for pyart.core.grid

"""
An class for holding gridded Radar data.

"""

import numpy as np
from netCDF4 import num2date

try:
    import xarray

    _XARRAY_AVAILABLE = True
except ImportError:
    _XARRAY_AVAILABLE = False

try:
    import pyproj

    _PYPROJ_AVAILABLE = True
except ImportError:
    _PYPROJ_AVAILABLE = False

from ..config import get_metadata
from ..exceptions import MissingOptionalDependency
from ..lazydict import LazyLoadDict
from .transforms import cartesian_to_geographic, cartesian_vectors_to_geographic


[docs]class Grid: """ A class for storing rectilinear gridded radar data in Cartesian coordinate. Refer to the attribute section for information on the parameters. To create a Grid object using legacy parameters present in Py-ART version 1.5 and before, use :py:func:`from_legacy_parameters`, grid = Grid.from_legacy_parameters(fields, axes, metadata). Attributes ---------- time : dict Time of the grid. fields : dict of dicts Moments from radars or other variables. metadata : dict Metadata describing the grid. origin_longitude, origin_latitude, origin_altitude : dict Geographic coordinate of the origin of the grid. x, y, z : dict, 1D Distance from the grid origin for each Cartesian coordinate axis in a one dimensional array. Defines the spacing along the three grid axes which is repeated throughout the grid, making a rectilinear grid. nx, ny, nz : int Number of grid points along the given Cartesian dimension. projection : dic or str Projection parameters defining the map projection used to transform from Cartesian to geographic coordinates. None will use the default dictionary with the 'proj' key set to 'pyart_aeqd' indicating that the native Py-ART azimuthal equidistant projection is used. Other values should specify a valid pyproj.Proj projparams dictionary or string. The special key '_include_lon_0_lat_0' is removed when interpreting this dictionary. If this key is present and set to True, which is required when proj='pyart_aeqd', then the radar longitude and latitude will be added to the dictionary as 'lon_0' and 'lat_0'. Use the :py:func:`get_projparams` method to retrieve a copy of this attribute dictionary with this special key evaluated. radar_longitude, radar_latitude, radar_altitude : dict or None, optional Geographic location of the radars which make up the grid. radar_time : dict or None, optional Start of collection for the radar which make up the grid. radar_name : dict or None, optional Names of the radars which make up the grid. nradar : int Number of radars whose data was used to make the grid. projection_proj : Proj pyproj.Proj instance for the projection specified by the projection attribute. If the 'pyart_aeqd' projection is specified accessing this attribute will raise a ValueError. point_x, point_y, point_z : LazyLoadDict The Cartesian locations of all grid points from the origin in the three Cartesian coordinates. The three dimensional data arrays contained these attributes are calculated from the x, y, and z attributes. If these attributes are changed use :py:func: `init_point_x_y_z` to reset the attributes. point_longitude, point_latitude : LazyLoadDict Geographic location of each grid point. The projection parameter(s) defined in the `projection` attribute are used to perform an inverse map projection from the Cartesian grid point locations relative to the grid origin. If these attributes are changed use :py:func:`init_point_longitude_latitude` to reset the attributes. point_altitude : LazyLoadDict The altitude of each grid point as calculated from the altitude of the grid origin and the Cartesian z location of each grid point. If this attribute is changed use :py:func:`init_point_altitude` to reset the attribute. """ def __init__( self, time, fields, metadata, origin_latitude, origin_longitude, origin_altitude, x, y, z, projection=None, radar_latitude=None, radar_longitude=None, radar_altitude=None, radar_time=None, radar_name=None, ): """Initalize object.""" self.time = time self.fields = fields self.metadata = metadata self.origin_latitude = origin_latitude self.origin_longitude = origin_longitude self.origin_altitude = origin_altitude self.x = x self.y = y self.z = z self.nx = len(x["data"]) self.ny = len(y["data"]) self.nz = len(z["data"]) if projection is None: self.projection = {"proj": "pyart_aeqd", "_include_lon_0_lat_0": True} else: self.projection = projection self.radar_latitude = radar_latitude self.radar_longitude = radar_longitude self.radar_altitude = radar_altitude self.radar_time = radar_time self.radar_name = radar_name self.nradar = self._find_and_check_nradar() # initialize attributes with Lazy load dictionaries self.init_point_x_y_z() self.init_point_longitude_latitude() self.init_point_altitude() return def __getstate__(self): """Return object's state which can be pickled.""" state = self.__dict__.copy() # copy the objects state # Remove unpicklable entries (those which are lazily loaded del state["point_x"] del state["point_y"] del state["point_z"] del state["point_latitude"] del state["point_longitude"] del state["point_altitude"] return state def __setstate__(self, state): """Restore unpicklable entries from pickled object.""" self.__dict__.update(state) self.init_point_x_y_z() self.init_point_longitude_latitude() self.init_point_altitude() @property def projection_proj(self): # Proj instance as specified by the projection attribute. # Raises a ValueError if the pyart_aeqd projection is specified. if not _PYPROJ_AVAILABLE: raise MissingOptionalDependency( "PyProj is required to create a Proj instance but it " + "is not installed" ) projparams = self.get_projparams() # Check if projparams is dictionary and check for pyart_aeqd if isinstance(projparams, dict): if projparams["proj"] == "pyart_aeqd": raise ValueError( "Proj instance can not be made for the pyart_aeqd projection" ) # Get proj instance from a proj str or dict proj = pyproj.Proj(projparams) return proj
[docs] def get_projparams(self): """Return a projparam dict or str from the projection attribute.""" if isinstance(self.projection, dict): projparams = self.projection.copy() if projparams.pop("_include_lon_0_lat_0", False): projparams["lon_0"] = self.origin_longitude["data"][0] projparams["lat_0"] = self.origin_latitude["data"][0] return projparams else: return self.projection
def _find_and_check_nradar(self): """ Return the number of radars which were used to create the grid. Examine the radar attributes to determine the number of radars which were used to create the grid. If the size of the radar attributes are inconsistent a ValueError is raised by this method. """ nradar_set = False nradar = 0 if self.radar_latitude is not None: nradar = len(self.radar_latitude["data"]) nradar_set = True if self.radar_longitude is not None: if nradar_set and len(self.radar_longitude["data"]) != nradar: raise ValueError("Inconsistent length of radar_ arguments.") nradar = len(self.radar_longitude["data"]) nradar_set = True if self.radar_altitude is not None: if nradar_set and len(self.radar_altitude["data"]) != nradar: raise ValueError("Inconsistent length of radar_ arguments.") nradar = len(self.radar_altitude["data"]) nradar_set = True if self.radar_time is not None: if nradar_set and len(self.radar_time["data"]) != nradar: raise ValueError("Inconsistent length of radar_ arguments.") nradar = len(self.radar_time["data"]) nradar_set = True if self.radar_name is not None: if nradar_set and len(self.radar_name["data"]) != nradar: raise ValueError("Inconsistent length of radar_ arguments.") nradar = len(self.radar_name["data"]) nradar_set = True return nradar # Attribute init/reset methods
[docs] def init_point_x_y_z(self): """Initialize or reset the point_{x, y, z} attributes.""" self.point_x = LazyLoadDict(get_metadata("point_x")) self.point_x.set_lazy("data", _point_data_factory(self, "x")) self.point_y = LazyLoadDict(get_metadata("point_y")) self.point_y.set_lazy("data", _point_data_factory(self, "y")) self.point_z = LazyLoadDict(get_metadata("point_z")) self.point_z.set_lazy("data", _point_data_factory(self, "z"))
[docs] def init_point_longitude_latitude(self): """ Initialize or reset the point_{longitude, latitudes} attributes. """ point_longitude = LazyLoadDict(get_metadata("point_longitude")) point_longitude.set_lazy("data", _point_lon_lat_data_factory(self, 0)) self.point_longitude = point_longitude point_latitude = LazyLoadDict(get_metadata("point_latitude")) point_latitude.set_lazy("data", _point_lon_lat_data_factory(self, 1)) self.point_latitude = point_latitude
[docs] def init_point_altitude(self): """Initialize the point_altitude attribute.""" point_altitude = LazyLoadDict(get_metadata("point_altitude")) point_altitude.set_lazy("data", _point_altitude_data_factory(self)) self.point_altitude = point_altitude
[docs] def write( self, filename, format="NETCDF4", arm_time_variables=False, arm_alt_lat_lon_variables=False, ): """ Write the the Grid object to a NetCDF file. Parameters ---------- filename : str Filename to save to. format : str, optional NetCDF format, one of 'NETCDF4', 'NETCDF4_CLASSIC', 'NETCDF3_CLASSIC' or 'NETCDF3_64BIT'. arm_time_variables : bool, optional True to write the ARM standard time variables base_time and time_offset. False will not write these variables. arm_alt_lat_lon_variables : bool, optional True to write the ARM standard alt, lat, lon variables. False will not write these variables. """ # delayed import to avoid circular import from ..io.grid_io import write_grid write_grid( filename, self, format=format, arm_time_variables=arm_time_variables, arm_alt_lat_lon_variables=arm_alt_lat_lon_variables, )
[docs] def to_xarray(self): """ Convert the Grid object to an xarray format. Attributes ---------- time : dict Time of the grid. fields : dict of dicts Moments from radars or other variables. longitude, latitude : dict, 2D Arrays of latitude and longitude for the grid height level. x, y, z : dict, 1D Distance from the grid origin for each Cartesian coordinate axis in a one dimensional array. """ if not _XARRAY_AVAILABLE: raise MissingOptionalDependency( "Xarray is required to use Grid.to_xarray but is not installed!" ) def _process_radar_name(radar_name): """Process radar_name to handle different formats.""" if radar_name.dtype.kind == "S" and radar_name.ndim > 1: # Join each row of bytes into a single byte string return np.array( [b"".join(row) for row in radar_name], dtype=f"|S{radar_name.shape[1]}", ) return radar_name lon, lat = self.get_point_longitude_latitude() z = self.z["data"] y = self.y["data"] x = self.x["data"] time = np.array([num2date(self.time["data"][0], units=self.time["units"])]) ds = xarray.Dataset() for field, field_info in self.fields.items(): field_data = field_info["data"] data = xarray.DataArray( np.ma.expand_dims(field_data, 0), dims=("time", "z", "y", "x"), coords={ "time": time, "z": z, "lat": (["y", "x"], lat), "lon": (["y", "x"], lon), "y": y, "x": x, }, ) data.attrs.update({k: v for k, v in field_info.items() if k != "data"}) ds[field] = data ds.lon.attrs = { "long_name": "longitude of grid cell center", "units": "degree_E", "standard_name": "Longitude", } ds.lat.attrs = { "long_name": "latitude of grid cell center", "units": "degree_N", "standard_name": "Latitude", } for attr in [ds.z, ds.lat, ds.lon]: attr.encoding["_FillValue"] = None from ..io.grid_io import _make_coordinatesystem_dict ds.coords["ProjectionCoordinateSystem"] = xarray.DataArray( data=np.array(1, dtype="int32"), attrs=_make_coordinatesystem_dict(self), ) projection = self.projection.copy() if "_include_lon_0_lat_0" in projection: projection["_include_lon_0_lat_0"] = str( projection["_include_lon_0_lat_0"] ).lower() ds.coords["projection"] = xarray.DataArray( data=np.array(1, dtype="int32"), attrs=projection, ) # Handle origin and radar attributes with appropriate dimensions for attr_name in [ "origin_latitude", "origin_longitude", "origin_altitude", ]: if hasattr(self, attr_name): attr_data = getattr(self, attr_name) if attr_data is not None: attr_value = np.ma.expand_dims(attr_data["data"][0], 0) ds.coords[attr_name] = xarray.DataArray( attr_value, dims=("time",), attrs=get_metadata(attr_name) ) # Radar-specific attributes that should have the nradar dimension for attr_name in [ "radar_altitude", "radar_latitude", "radar_longitude", "radar_time", ]: if hasattr(self, attr_name): attr_data = getattr(self, attr_name) if attr_data is not None: ds.coords[attr_name] = xarray.DataArray( attr_data["data"], dims=("nradar",), attrs=get_metadata(attr_name), ) if "radar_time" in ds.variables: ds.radar_time.attrs.pop("calendar") # Handle radar_name and ensure it has the correct dimension if self.radar_name is not None: radar_name = _process_radar_name(self.radar_name["data"]) ds.coords["radar_name"] = xarray.DataArray( radar_name, dims=("nradar",), attrs=get_metadata("radar_name") ) else: radar_name = np.array(["ExampleRadar"], dtype="U") ds.coords["radar_name"] = xarray.DataArray( radar_name, dims=("nradar",), attrs=get_metadata("radar_name") ) # Add radar_name to attributes ds.attrs["radar_name"] = ( radar_name.tolist() if radar_name.size > 1 else radar_name.item() ) ds.attrs["nradar"] = radar_name.size ds.attrs.update(self.metadata) for key in ds.attrs: try: ds.attrs[key] = ds.attrs[key].decode("utf-8") except AttributeError: pass ds.close() return ds
[docs] def add_field(self, field_name, field_dict, replace_existing=False): """ Add a field to the object. Parameters ---------- field_name : str Name of the field to the fields dictionary. field_dict : dict Dictionary containing field data and metadata. replace_existing : bool, optional True to replace the existing field with key field_name if it exists, overwriting the existing data. If False, a ValueError is raised if field_name already exists. """ # checks to make sure input field dictionary is valid if "data" not in field_dict: raise KeyError('Field dictionary must contain a "data" key') if field_name in self.fields and replace_existing is False: raise ValueError(f"A field named {field_name} already exists") if field_dict["data"].shape != (self.nz, self.ny, self.nx): raise ValueError("Field has invalid shape") self.fields[field_name] = field_dict
[docs] def get_point_longitude_latitude(self, level=0, edges=False): """ Return arrays of longitude and latitude for a given grid height level. Parameters ---------- level : int, optional Grid height level at which to determine latitudes and longitudes. This is not currently used as all height level have the same layout. edges : bool, optional True to calculate the latitude and longitudes of the edges by interpolating between Cartesian coordinates points and extrapolating at the boundaries. False to calculate the locations at the centers. Returns ------- longitude, latitude : 2D array Arrays containing the latitude and longitudes, in degrees, of the grid points or edges between grid points for the given height. """ x = self.x["data"] y = self.y["data"] projparams = self.get_projparams() return cartesian_vectors_to_geographic(x, y, projparams, edges=edges)
def _point_data_factory(grid, coordinate): """Return a function which returns the locations of all points.""" def _point_data(): """The function which returns the locations of all points.""" reg_x = grid.x["data"] reg_y = grid.y["data"] reg_z = grid.z["data"] if coordinate == "x": return np.tile(reg_x, (len(reg_z), len(reg_y), 1)).swapaxes(2, 2) elif coordinate == "y": return np.tile(reg_y, (len(reg_z), len(reg_x), 1)).swapaxes(1, 2) else: assert coordinate == "z" return np.tile(reg_z, (len(reg_x), len(reg_y), 1)).swapaxes(0, 2) return _point_data def _point_lon_lat_data_factory(grid, coordinate): """Return a function which returns the geographic locations of points.""" def _point_lon_lat_data(): """The function which returns the geographic point locations.""" x = grid.point_x["data"] y = grid.point_y["data"] projparams = grid.get_projparams() geographic_coords = cartesian_to_geographic(x, y, projparams) # Set point_latitude['data'] when point_longitude['data'] is evaluated # and vice-versa. This ensures that both attributes contain data from # the same map projection and that the map projection only needs to be # evaluated once. if coordinate == 0: grid.point_latitude["data"] = geographic_coords[1] else: grid.point_longitude["data"] = geographic_coords[0] return geographic_coords[coordinate] return _point_lon_lat_data def _point_altitude_data_factory(grid): """Return a function which returns the point altitudes.""" def _point_altitude_data(): """The function which returns the point altitudes.""" return grid.origin_altitude["data"][0] + grid.point_z["data"] return _point_altitude_data