Source code for pyremap.polar

# This software is open source software available under the BSD-3 license.
#
# Copyright (c) 2019 Triad National Security, LLC. All rights reserved.
# Copyright (c) 2019 Lawrence Livermore National Security, LLC. All rights
# reserved.
# Copyright (c) 2019 UT-Battelle, LLC. All rights reserved.
#
# Additional copyright and license information can be found in the LICENSE file
# distributed with this code, or at
# https://raw.githubusercontent.com/MPAS-Dev/MPAS-Analysis/master/LICENSE
import pyproj
import xarray
import numpy

from pyremap.descriptor import ProjectionGridDescriptor


[docs]def get_arctic_stereographic_projection(): """ Get a projection for an Arctic stereographic comparison grid Returns ------- projection : ``pyproj.Proj`` object The projection """ # Authors # ------- # Milena Veneziani projection = pyproj.Proj('+proj=stere +lat_ts=75.0 +lat_0=90 +lon_0=0.0 ' '+k_0=1.0 +x_0=0.0 +y_0=0.0 +ellps=WGS84') return projection
[docs]def get_antarctic_stereographic_projection(): """ Get a projection for an Antarctic steregraphic grid """ projection = pyproj.Proj('+proj=stere +lat_ts=-71.0 +lat_0=-90 +lon_0=0.0 ' '+k_0=1.0 +x_0=0.0 +y_0=0.0 +ellps=WGS84') return projection
[docs]def get_polar_descriptor_from_file(fileName, projection='antarctic'): """ Get a descriptor of a polar stereographic grid used for remapping Parameters ---------- fileName : str A file containing x and y coordinates for the grid projection : {'arctic', 'antarctic', pyproj.Proj} The projection to use. 'arctic' and 'antarctic' are polar stereographic projections with reference latitude at +/- 71 degrees Returns ------- descriptor : ``ProjectionGridDescriptor`` object A descriptor of the polar grid """ dsIn = xarray.open_dataset(fileName) x = dsIn.x.values y = dsIn.y.values dx = int((x[1]-x[0])/1000.) Lx = int((x[-1] - x[0])/1000.) Ly = int((y[-1] - y[0])/1000.) meshName = '{}x{}km_{}km_Antarctic_stereo'.format(Lx, Ly, dx) projection = _get_projection(projection) descriptor = ProjectionGridDescriptor.create(projection, x, y, meshName) return descriptor
[docs]def get_polar_descriptor(Lx, Ly, dx, dy, projection='antarctic'): """ Get a descriptor of a polar stereographic grid used for remapping Parameters ---------- Lx, Ly : double Size of the domain in x and y in km dx, dy : double Resolution of the grid in km projection : {'arctic', 'antarctic', pyproj.Proj} The projection to use. 'arctic' and 'antarctic' are polar stereographic projections with reference latitude at +/- 71 degrees Returns ------- descriptor : ``ProjectionGridDescriptor`` object A descriptor of the Antarctic grid """ upperProj = projection[0].upper() + projection[1:] meshName = '{}x{}km_{}km_{}_stereo'.format(Lx, Ly, dx, upperProj) xMax = 0.5 * Lx * 1e3 nx = int(Lx / dx) + 1 x = numpy.linspace(-xMax, xMax, nx) yMax = 0.5 * Ly * 1e3 ny = int(Ly / dy) + 1 y = numpy.linspace(-yMax, yMax, ny) projection = _get_projection(projection) descriptor = ProjectionGridDescriptor.create(projection, x, y, meshName) return descriptor
[docs]def to_polar(points): projection = get_antarctic_stereographic_projection() latLonProjection = pyproj.Proj(proj='latlong', datum='WGS84') transformer = pyproj.Transformer.from_proj(latLonProjection, projection) x, y = transformer.transform(points[:, 0], points[:, 1], radians=False) points[:, 0] = x points[:, 1] = y return points
[docs]def from_polar(points): projection = get_antarctic_stereographic_projection() latLonProjection = pyproj.Proj(proj='latlong', datum='WGS84') transformer = pyproj.Transformer.from_proj(projection, latLonProjection) lon, lat = transformer.transform(points[:, 0], points[:, 1], radians=False) points[:, 0] = lon points[:, 1] = lat return points
def _get_projection(projection): if isinstance(projection, str): if projection == 'arctic': projection = get_arctic_stereographic_projection() elif projection == 'antarctic': projection = get_antarctic_stereographic_projection() else: raise ValueError('Bad projection name {}'.format(projection)) return projection