From edbae6460ba188b6d87117222f6c0891bcbc56ef Mon Sep 17 00:00:00 2001 From: Zach Date: Sun, 10 Jan 2021 13:31:25 -0500 Subject: [PATCH 1/6] add vectorfield draft --- packetraven/model.py | 322 +++++++++++++++++++++++++++++++++++++++++++ setup.py | 2 + 2 files changed, 324 insertions(+) create mode 100644 packetraven/model.py diff --git a/packetraven/model.py b/packetraven/model.py new file mode 100644 index 00000000..427397d1 --- /dev/null +++ b/packetraven/model.py @@ -0,0 +1,322 @@ +from datetime import datetime, timedelta + +from matplotlib import pyplot, quiver +import numpy +import pyproj +from pyproj import CRS, Transformer +import xarray as xarray + +WGS84 = CRS.from_epsg(4326) +WEB_MERCATOR = CRS.from_epsg(3857) + + +class VectorField: + """ + Vector field of (u, v) values. + """ + + def __init__(self, time_deltas: [timedelta], projection: CRS = None): + """ + Build vector field of (u, v) values. + + :param time_deltas: list of time deltas + :param projection: native projection of field + """ + + self.time_deltas = [numpy.timedelta64(time_delta) for time_delta in time_deltas] + self.projection = projection if projection is not None else WGS84 + + time_delta = numpy.nanmean(self.time_deltas).item() + if type(time_delta) is int: + time_delta *= 1e-9 + elif type(time_delta) is timedelta: + time_delta = int(time_delta / timedelta(seconds=1)) + + self.delta_t = timedelta(seconds=time_delta) + + def u(self, point: numpy.array, time: datetime) -> float: + """ + u velocity in m/s at coordinates + + :param point: coordinates in linear units + :param time: time + :return: u value at coordinate in m/s + """ + + pass + + def v(self, point: numpy.array, time: datetime) -> float: + """ + v velocity in m/s at coordinates + + :param point: lon / lat coordinates + :param time: time + :return: v value at coordinate in m/s + """ + + pass + + def velocity(self, point: numpy.array, time: datetime) -> float: + """ + absolute velocity in m/s at coordinate + + :param point: lon / lat coordinates + :param time: time + :return: magnitude of uv vector in m/s + """ + + return math.sqrt(self.u(point, time) ** 2 + self.v(point, time) ** 2) + + def direction(self, point: numpy.array, time: datetime) -> float: + """ + angle of uv vector + + :param point: lon / lat coordinates + :param time: time + :return: radians from east of uv vector + """ + + return math.atan2(self.u(point, time), self.v(point, time)) + + def plot(self, time: datetime, axis: pyplot.Axes = None, **kwargs) -> quiver.Quiver: + """ + Plot vector field at the given time. + + :param time: time at which to plot + :param axis: pyplot axis on which to plot + :return: quiver plot + """ + + pass + + def __getitem__(self, position: (numpy.array, datetime)) -> numpy.array: + """ + velocity vector (u, v) in m/s at coordinates + + :param point: coordinates in linear units + :param time: time + :return: (u, v) vector + """ + + point, time = position + vector = numpy.array([self.u(point, time), self.v(point, time)]) + vector[numpy.isnan(vector)] = 0 + return vector + + def __repr__(self) -> str: + return f'{self.__class__.__name__}({self.delta_t})' + + +class RankineVortex(VectorField): + """ + Time-invariant circular vector field moving in solid-body rotation (like a turntable). + """ + + def __init__( + self, + center: (float, float), + radius: float, + period: timedelta, + time_deltas: numpy.array, + ): + """ + Construct a 2-dimensional solid rotating disk surrounded by inverse falloff of tangential velocity. + + :param center: tuple of geographic coordinates + :param radius: radius of central solid rotation in meters + :param period: rotational period + :param time_deltas: time differences + """ + + transformer = Transformer.from_crs(WGS84, WEB_MERCATOR) + self.center = numpy.array(transformer.transform(*center)) + self.radius = radius + self.angular_velocity = 2 * numpy.pi / (period / timedelta(seconds=1)) + + super().__init__(time_deltas) + + def u(self, point: numpy.array, time: datetime) -> float: + return -self.velocity(point, time) * numpy.cos(numpy.atan2(*(point - self.center))) + + def v(self, point: numpy.array, time: datetime) -> float: + return self.velocity(point, time) * numpy.sin(numpy.atan2(*(point - self.center))) + + def velocity(self, point: numpy.array, time: datetime) -> float: + radial_distance = numpy.sqrt(numpy.sum((point - self.center) ** 2)) + + if radial_distance <= self.radius: + return self.angular_velocity * radial_distance + else: + return self.angular_velocity * self.radius ** 2 / radial_distance + + def plot(self, axis: pyplot.Axes = None, **kwargs) -> quiver.Quiver: + if axis is None: + axis = pyplot.axes(projection=cartopy.crs.PlateCarree()) + + points = [] + radii = numpy.linspace(1, self.radius * 2, 20) + + for radius in radii: + num_points = 50 + points.extend( + [ + ( + numpy.cos(2 * numpy.pi / num_points * point_index) * radius + + self.center[0], + numpy.sin(2 * numpy.pi / num_points * point_index) * radius + + self.center[1], + ) + for point_index in range(0, num_points + 1) + ] + ) + + vectors = [self[point, datetime.now()] for point in points] + points = list( + zip(*pyproj.transform(WEB_MERCATOR, WGS84, *zip(*points))) + ) + + quiver_plot = axis.quiver(*zip(*points), *zip(*vectors), units='width', **kwargs) + axis.quiverkey( + quiver_plot, 0.9, 0.9, 1, r'$1 \frac{m}{s}$', labelpos='E', coordinates='figure' + ) + + return quiver_plot + + +class VectorDataset(VectorField): + """ + Vector field with time component using xarray observation. + """ + + def __init__( + self, + dataset: xarray.Dataset, + u_name: str = 'u', + v_name: str = 'v', + x_name: str = 'lon', + y_name: str = 'lat', + t_name: str = 'time', + coordinate_system: pyproj.Proj = None, + ): + """ + Create new velocity field from given observation. + + :param dataset: xarray observation containing velocity data (u, v) + :param u_name: name of u variable + :param v_name: name of v variable + :param x_name: name of x coordinate + :param y_name: name of y coordinate + :param t_name: name of time coordinate + :param coordinate_system: coordinate system of observation + """ + + self.coordinate_system = ( + coordinate_system if coordinate_system is not None else utilities.WGS84 + ) + + variables_to_rename = { + u_name: 'u', + v_name: 'v', + x_name: 'x', + y_name: 'y', + t_name: 'time', + } + self.dataset = dataset.rename(variables_to_rename) + + x, y = pyproj.transform( + self.coordinate_system, + utilities.WEB_MERCATOR, + *numpy.meshgrid(self.dataset['x'].values, self.dataset['y'].values), + ) + + self.dataset['x'] = x[0, :] + self.dataset['y'] = y[:, 0] + + self.delta_x = numpy.mean(numpy.diff(self.dataset['x'])) + self.delta_y = numpy.mean(numpy.diff(self.dataset['y'])) + + super().__init__(numpy.diff(self.dataset['time'].values)) + + def _interpolate( + self, variable: str, point: numpy.array, time: datetime + ) -> xarray.DataArray: + transformed_point = pyproj.transform( + utilities.WEB_MERCATOR, self.coordinate_system, point[0], point[1] + ) + + x_name = f'{variable}_x' + y_name = f'{variable}_y' + + x_range = slice( + self.dataset[x_name] + .sel({x_name: numpy.min(transformed_point[0]) - 1}, method='bfill') + .values.item(), + self.dataset[x_name] + .sel({x_name: numpy.max(transformed_point[0]) + 1}, method='ffill') + .values.item(), + ) + y_range = slice( + self.dataset[y_name] + .sel({y_name: numpy.min(transformed_point[1]) - 1}, method='bfill') + .values.item(), + self.dataset[y_name] + .sel({y_name: numpy.max(transformed_point[1]) + 1}, method='ffill') + .values.item(), + ) + time_range = slice( + self.dataset['time'].sel(time=time, method='bfill').values, + self.dataset['time'].sel(time=time, method='ffill').values, + ) + + if time_range.start == time_range.stop: + time_range = time_range.start + + cell = self.dataset[variable].sel( + {'time': time_range, x_name: x_range, y_name: y_range} + ) + + if len(transformed_point.shape) > 1: + cell = cell.interp({'time': time}) if 'time' in cell.dims else cell + return xarray.concat( + [ + cell.interp({x_name: location[0], y_name: location[1]}) + for location in transformed_point.T + ], + dim='point', + ) + else: + cell = cell.interp({x_name: transformed_point[0], y_name: transformed_point[1]}) + return cell.interp({'time': time}) if 'time' in cell.dims else cell + + def u(self, point: numpy.array, time: datetime) -> float: + return self._interpolate('u', point, time).values + + def v(self, point: numpy.array, time: datetime) -> float: + return self._interpolate('v', point, time).values + + def plot(self, time: datetime, axis: pyplot.Axes = None, **kwargs) -> quiver.Quiver: + if time is None: + time = self.dataset['time'].values[0] + + if axis is None: + axis = pyplot.axes(projection=cartopy.crs.PlateCarree()) + + lon, lat = pyproj.transform( + utilities.WEB_MERCATOR, + utilities.WGS84, + *numpy.meshgrid(self.dataset['x'].values, self.dataset['y'].values), + ) + + quiver_plot = axis.quiver( + lon, + lat, + self.dataset['u'].sel(time=time, method='nearest'), + self.dataset['v'].sel(time=time, method='nearest'), + units='width', + **kwargs, + ) + axis.quiverkey( + quiver_plot, 0.9, 0.9, 1, r'$1 \frac{m}{s}$', labelpos='E', coordinates='figure' + ) + + return quiver_plot diff --git a/setup.py b/setup.py index 453742af..e7d2d7dd 100644 --- a/setup.py +++ b/setup.py @@ -38,6 +38,7 @@ install_requires=[ 'aprslib', 'haversine', + 'netcdf4', 'numpy==1.19.3', 'pyserial', 'geojson', @@ -49,6 +50,7 @@ 'shapely', 'sshtunnel', 'tablecrow', + 'xarray', ], extras_require={'testing': ['flake8', 'pytest', 'pytest-cov', 'pytz'], 'development': ['oitnb']}, entry_points={'console_scripts': ['packetraven=client.cli:main']}, From eef4e0b3ed7e94aa73155427bab993a93e94a558 Mon Sep 17 00:00:00 2001 From: Zach Date: Sat, 30 Jan 2021 20:56:12 -0500 Subject: [PATCH 2/6] workflow updates --- .github/workflows/build.yml | 20 +++++++++----------- .github/workflows/tests.yml | 10 ++++++---- 2 files changed, 15 insertions(+), 15 deletions(-) diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index 05f4a347..3e38caf5 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -7,12 +7,12 @@ on: jobs: build_wheels: - name: build wheel on ${{ matrix.os }} + name: Python ${{ matrix.python-version }} wheel on ${{ matrix.os }} runs-on: ${{ matrix.os }} strategy: matrix: - os: [ubuntu-latest] - python-version: [3.8] + os: [ ubuntu-latest ] + python-version: [ 3.9 ] steps: - name: Checkout repository uses: actions/checkout@v2 @@ -24,8 +24,8 @@ jobs: uses: actions/cache@v2 with: path: ~/.cache/pip - key: ${{ matrix.os }}-pip-${{ matrix.python-version }}-${{ hashFiles('setup.py') }} - restore-keys: ${{ matrix.os }}-pip-${{ matrix.python-version }}- + key: ${{ runner.os }}-pip-${{ matrix.python-version }}-${{ hashFiles('setup.py') }} + restore-keys: ${{ runner.os }}-pip-${{ matrix.python-version }}- - name: Install dependencies run: pip install --upgrade pip setuptools wheel - name: Build wheel @@ -42,8 +42,6 @@ jobs: uses: actions/checkout@v2 - name: Install Python uses: actions/setup-python@v2 - with: - python-version: '3.8' - name: Package source run: python setup.py sdist - name: Save source package @@ -51,16 +49,16 @@ jobs: with: path: ./dist/*.tar.gz upload_pypi: - name: publish wheel to PyPI - needs: [build_wheels, build_sdist] + name: publish to PyPI + needs: [ build_wheels, build_sdist ] runs-on: ubuntu-latest steps: - - name: Retrieve wheel + - name: Retrieve wheel(s) and source uses: actions/download-artifact@v2 with: name: artifact path: dist - - name: Upload wheel + - name: Upload wheel(s) and source uses: pypa/gh-action-pypi-publish@master with: user: __token__ diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index 6c9142f2..be262b95 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -4,12 +4,12 @@ on: [ push ] jobs: tests: - name: test with Python ${{ matrix.python-version }} on ${{ matrix.os }} + name: Python ${{ matrix.python-version }} on ${{ matrix.os }} runs-on: ${{ matrix.os }} strategy: matrix: - os: [ ubuntu-latest ] - python-version: [ 3.8 ] + os: [ ubuntu-latest, windows-latest ] + python-version: [ 3.8, 3.9 ] services: postgres: image: postgis/postgis:10-2.5 @@ -32,9 +32,11 @@ jobs: path: ~/.cache/pip key: ${{ runner.os }}-pip-${{ matrix.python-version }}-${{ hashFiles('setup.py') }} restore-keys: ${{ runner.os }}-pip-${{ matrix.python-version }}- + - name: Update pip + run: python -m pip install --upgrade pip - name: Install dependencies run: | - pip install --upgrade pip + pip install wheel pip install -e .[testing] - name: Lint with flake8 run: | From 4360bc739891264c4e5c0fc32bdaddb7e81d6446 Mon Sep 17 00:00:00 2001 From: Zach Date: Sun, 10 Jan 2021 13:31:25 -0500 Subject: [PATCH 3/6] add vectorfield draft --- packetraven/model.py | 320 +++++++++++++++++++++++++++++++++++++++++++ setup.py | 2 + 2 files changed, 322 insertions(+) diff --git a/packetraven/model.py b/packetraven/model.py index 7883661d..b6e40fd5 100644 --- a/packetraven/model.py +++ b/packetraven/model.py @@ -1,4 +1,11 @@ +from matplotlib import pyplot, quiver import numpy +import pyproj +from pyproj import CRS, Transformer +import xarray as xarray + +WGS84 = CRS.from_epsg(4326) +WEB_MERCATOR = CRS.from_epsg(3857) # `dh/dt` based on historical flight data FREEFALL_DESCENT_RATE = lambda altitude: -5.8e-08 * altitude ** 2 - 6.001 @@ -8,3 +15,316 @@ # integration of `(1/(dh/dt)) dh` based on historical flight data # TODO make this model better with ML FREEFALL_SECONDS_TO_GROUND = lambda altitude: 1695.02 * numpy.arctan(9.8311e-5 * altitude) +from datetime import datetime, timedelta + + +class VectorField: + """ + Vector field of (u, v) values. + """ + + def __init__(self, time_deltas: [timedelta], projection: CRS = None): + """ + Build vector field of (u, v) values. + + :param time_deltas: list of time deltas + :param projection: native projection of field + """ + + self.time_deltas = [numpy.timedelta64(time_delta) for time_delta in time_deltas] + self.projection = projection if projection is not None else WGS84 + + time_delta = numpy.nanmean(self.time_deltas).item() + if type(time_delta) is int: + time_delta *= 1e-9 + elif type(time_delta) is timedelta: + time_delta = int(time_delta / timedelta(seconds=1)) + + self.delta_t = timedelta(seconds=time_delta) + + def u(self, point: numpy.array, time: datetime) -> float: + """ + u velocity in m/s at coordinates + + :param point: coordinates in linear units + :param time: time + :return: u value at coordinate in m/s + """ + + pass + + def v(self, point: numpy.array, time: datetime) -> float: + """ + v velocity in m/s at coordinates + + :param point: lon / lat coordinates + :param time: time + :return: v value at coordinate in m/s + """ + + pass + + def velocity(self, point: numpy.array, time: datetime) -> float: + """ + absolute velocity in m/s at coordinate + + :param point: lon / lat coordinates + :param time: time + :return: magnitude of uv vector in m/s + """ + + return math.sqrt(self.u(point, time) ** 2 + self.v(point, time) ** 2) + + def direction(self, point: numpy.array, time: datetime) -> float: + """ + angle of uv vector + + :param point: lon / lat coordinates + :param time: time + :return: radians from east of uv vector + """ + + return math.atan2(self.u(point, time), self.v(point, time)) + + def plot(self, time: datetime, axis: pyplot.Axes = None, **kwargs) -> quiver.Quiver: + """ + Plot vector field at the given time. + + :param time: time at which to plot + :param axis: pyplot axis on which to plot + :return: quiver plot + """ + + pass + + def __getitem__(self, position: (numpy.array, datetime)) -> numpy.array: + """ + velocity vector (u, v) in m/s at coordinates + + :param point: coordinates in linear units + :param time: time + :return: (u, v) vector + """ + + point, time = position + vector = numpy.array([self.u(point, time), self.v(point, time)]) + vector[numpy.isnan(vector)] = 0 + return vector + + def __repr__(self) -> str: + return f'{self.__class__.__name__}({self.delta_t})' + + +class RankineVortex(VectorField): + """ + Time-invariant circular vector field moving in solid-body rotation (like a turntable). + """ + + def __init__( + self, + center: (float, float), + radius: float, + period: timedelta, + time_deltas: numpy.array, + ): + """ + Construct a 2-dimensional solid rotating disk surrounded by inverse falloff of tangential velocity. + + :param center: tuple of geographic coordinates + :param radius: radius of central solid rotation in meters + :param period: rotational period + :param time_deltas: time differences + """ + + transformer = Transformer.from_crs(WGS84, WEB_MERCATOR) + self.center = numpy.array(transformer.transform(*center)) + self.radius = radius + self.angular_velocity = 2 * numpy.pi / (period / timedelta(seconds=1)) + + super().__init__(time_deltas) + + def u(self, point: numpy.array, time: datetime) -> float: + return -self.velocity(point, time) * numpy.cos(numpy.atan2(*(point - self.center))) + + def v(self, point: numpy.array, time: datetime) -> float: + return self.velocity(point, time) * numpy.sin(numpy.atan2(*(point - self.center))) + + def velocity(self, point: numpy.array, time: datetime) -> float: + radial_distance = numpy.sqrt(numpy.sum((point - self.center) ** 2)) + + if radial_distance <= self.radius: + return self.angular_velocity * radial_distance + else: + return self.angular_velocity * self.radius ** 2 / radial_distance + + def plot(self, axis: pyplot.Axes = None, **kwargs) -> quiver.Quiver: + if axis is None: + axis = pyplot.axes(projection=cartopy.crs.PlateCarree()) + + points = [] + radii = numpy.linspace(1, self.radius * 2, 20) + + for radius in radii: + num_points = 50 + points.extend( + [ + ( + numpy.cos(2 * numpy.pi / num_points * point_index) * radius + + self.center[0], + numpy.sin(2 * numpy.pi / num_points * point_index) * radius + + self.center[1], + ) + for point_index in range(0, num_points + 1) + ] + ) + + vectors = [self[point, datetime.now()] for point in points] + points = list( + zip(*pyproj.transform(WEB_MERCATOR, WGS84, *zip(*points))) + ) + + quiver_plot = axis.quiver(*zip(*points), *zip(*vectors), units='width', **kwargs) + axis.quiverkey( + quiver_plot, 0.9, 0.9, 1, r'$1 \frac{m}{s}$', labelpos='E', coordinates='figure' + ) + + return quiver_plot + + +class VectorDataset(VectorField): + """ + Vector field with time component using xarray observation. + """ + + def __init__( + self, + dataset: xarray.Dataset, + u_name: str = 'u', + v_name: str = 'v', + x_name: str = 'lon', + y_name: str = 'lat', + t_name: str = 'time', + coordinate_system: pyproj.Proj = None, + ): + """ + Create new velocity field from given observation. + + :param dataset: xarray observation containing velocity data (u, v) + :param u_name: name of u variable + :param v_name: name of v variable + :param x_name: name of x coordinate + :param y_name: name of y coordinate + :param t_name: name of time coordinate + :param coordinate_system: coordinate system of observation + """ + + self.coordinate_system = ( + coordinate_system if coordinate_system is not None else utilities.WGS84 + ) + + variables_to_rename = { + u_name: 'u', + v_name: 'v', + x_name: 'x', + y_name: 'y', + t_name: 'time', + } + self.dataset = dataset.rename(variables_to_rename) + + x, y = pyproj.transform( + self.coordinate_system, + utilities.WEB_MERCATOR, + *numpy.meshgrid(self.dataset['x'].values, self.dataset['y'].values), + ) + + self.dataset['x'] = x[0, :] + self.dataset['y'] = y[:, 0] + + self.delta_x = numpy.mean(numpy.diff(self.dataset['x'])) + self.delta_y = numpy.mean(numpy.diff(self.dataset['y'])) + + super().__init__(numpy.diff(self.dataset['time'].values)) + + def _interpolate( + self, variable: str, point: numpy.array, time: datetime + ) -> xarray.DataArray: + transformed_point = pyproj.transform( + utilities.WEB_MERCATOR, self.coordinate_system, point[0], point[1] + ) + + x_name = f'{variable}_x' + y_name = f'{variable}_y' + + x_range = slice( + self.dataset[x_name] + .sel({x_name: numpy.min(transformed_point[0]) - 1}, method='bfill') + .values.item(), + self.dataset[x_name] + .sel({x_name: numpy.max(transformed_point[0]) + 1}, method='ffill') + .values.item(), + ) + y_range = slice( + self.dataset[y_name] + .sel({y_name: numpy.min(transformed_point[1]) - 1}, method='bfill') + .values.item(), + self.dataset[y_name] + .sel({y_name: numpy.max(transformed_point[1]) + 1}, method='ffill') + .values.item(), + ) + time_range = slice( + self.dataset['time'].sel(time=time, method='bfill').values, + self.dataset['time'].sel(time=time, method='ffill').values, + ) + + if time_range.start == time_range.stop: + time_range = time_range.start + + cell = self.dataset[variable].sel( + {'time': time_range, x_name: x_range, y_name: y_range} + ) + + if len(transformed_point.shape) > 1: + cell = cell.interp({'time': time}) if 'time' in cell.dims else cell + return xarray.concat( + [ + cell.interp({x_name: location[0], y_name: location[1]}) + for location in transformed_point.T + ], + dim='point', + ) + else: + cell = cell.interp({x_name: transformed_point[0], y_name: transformed_point[1]}) + return cell.interp({'time': time}) if 'time' in cell.dims else cell + + def u(self, point: numpy.array, time: datetime) -> float: + return self._interpolate('u', point, time).values + + def v(self, point: numpy.array, time: datetime) -> float: + return self._interpolate('v', point, time).values + + def plot(self, time: datetime, axis: pyplot.Axes = None, **kwargs) -> quiver.Quiver: + if time is None: + time = self.dataset['time'].values[0] + + if axis is None: + axis = pyplot.axes(projection=cartopy.crs.PlateCarree()) + + lon, lat = pyproj.transform( + utilities.WEB_MERCATOR, + utilities.WGS84, + *numpy.meshgrid(self.dataset['x'].values, self.dataset['y'].values), + ) + + quiver_plot = axis.quiver( + lon, + lat, + self.dataset['u'].sel(time=time, method='nearest'), + self.dataset['v'].sel(time=time, method='nearest'), + units='width', + **kwargs, + ) + axis.quiverkey( + quiver_plot, 0.9, 0.9, 1, r'$1 \frac{m}{s}$', labelpos='E', coordinates='figure' + ) + + return quiver_plot diff --git a/setup.py b/setup.py index 540271f5..21cd25aa 100644 --- a/setup.py +++ b/setup.py @@ -80,6 +80,7 @@ install_requires=[ 'aprslib', 'haversine', + 'netcdf4', 'numpy>=1.20.0', 'pandas', 'pyserial', @@ -92,6 +93,7 @@ 'shapely', 'sshtunnel', 'tablecrow>=1.3.3', + 'xarray', ], extras_require={ 'testing': ['flake8', 'pytest', 'pytest-cov', 'pytz'], From cb6cd6bbb8b02df675ee517117cd55b7e5290f6f Mon Sep 17 00:00:00 2001 From: Zach Date: Tue, 13 Apr 2021 19:24:27 -0400 Subject: [PATCH 4/6] add GFS class --- packetraven/model.py | 64 ++++++++++++++++++-------------------------- 1 file changed, 26 insertions(+), 38 deletions(-) diff --git a/packetraven/model.py b/packetraven/model.py index b6e40fd5..e388bc60 100644 --- a/packetraven/model.py +++ b/packetraven/model.py @@ -203,6 +203,7 @@ def __init__( v_name: str = 'v', x_name: str = 'lon', y_name: str = 'lat', + z_name: str = 'z', t_name: str = 'time', coordinate_system: pyproj.Proj = None, ): @@ -214,12 +215,15 @@ def __init__( :param v_name: name of v variable :param x_name: name of x coordinate :param y_name: name of y coordinate + :param z_name: name of z coordinate :param t_name: name of time coordinate :param coordinate_system: coordinate system of observation """ self.coordinate_system = ( - coordinate_system if coordinate_system is not None else utilities.WGS84 + coordinate_system + if coordinate_system is not None + else WGS84 ) variables_to_rename = { @@ -229,11 +233,13 @@ def __init__( y_name: 'y', t_name: 'time', } + if z_name in dataset: + variables_to_rename[z_name] = 'z' self.dataset = dataset.rename(variables_to_rename) x, y = pyproj.transform( self.coordinate_system, - utilities.WEB_MERCATOR, + WEB_MERCATOR, *numpy.meshgrid(self.dataset['x'].values, self.dataset['y'].values), ) @@ -249,27 +255,19 @@ def _interpolate( self, variable: str, point: numpy.array, time: datetime ) -> xarray.DataArray: transformed_point = pyproj.transform( - utilities.WEB_MERCATOR, self.coordinate_system, point[0], point[1] + WEB_MERCATOR, self.coordinate_system, point[0], point[1] ) x_name = f'{variable}_x' y_name = f'{variable}_y' x_range = slice( - self.dataset[x_name] - .sel({x_name: numpy.min(transformed_point[0]) - 1}, method='bfill') - .values.item(), - self.dataset[x_name] - .sel({x_name: numpy.max(transformed_point[0]) + 1}, method='ffill') - .values.item(), + self.dataset[x_name].sel({x_name: numpy.min(transformed_point[0]) - 1}, method='bfill').values.item(), + self.dataset[x_name].sel({x_name: numpy.max(transformed_point[0]) + 1}, method='ffill').values.item(), ) y_range = slice( - self.dataset[y_name] - .sel({y_name: numpy.min(transformed_point[1]) - 1}, method='bfill') - .values.item(), - self.dataset[y_name] - .sel({y_name: numpy.max(transformed_point[1]) + 1}, method='ffill') - .values.item(), + self.dataset[y_name].sel({y_name: numpy.min(transformed_point[1]) - 1}, method='bfill').values.item(), + self.dataset[y_name].sel({y_name: numpy.max(transformed_point[1]) + 1}, method='ffill').values.item(), ) time_range = slice( self.dataset['time'].sel(time=time, method='bfill').values, @@ -280,7 +278,11 @@ def _interpolate( time_range = time_range.start cell = self.dataset[variable].sel( - {'time': time_range, x_name: x_range, y_name: y_range} + { + 'time': time_range, + x_name: x_range, + y_name: y_range, + } ) if len(transformed_point.shape) > 1: @@ -302,29 +304,15 @@ def u(self, point: numpy.array, time: datetime) -> float: def v(self, point: numpy.array, time: datetime) -> float: return self._interpolate('v', point, time).values - def plot(self, time: datetime, axis: pyplot.Axes = None, **kwargs) -> quiver.Quiver: - if time is None: - time = self.dataset['time'].values[0] - if axis is None: - axis = pyplot.axes(projection=cartopy.crs.PlateCarree()) +class VectorGFS(VectorDataset): + opendap_url = 'http://nomads.ncep.noaa.gov:80/dods/gfs_0p25_1hr/gfs20210413/gfs_0p25_1hr_12z' - lon, lat = pyproj.transform( - utilities.WEB_MERCATOR, - utilities.WGS84, - *numpy.meshgrid(self.dataset['x'].values, self.dataset['y'].values), - ) + def __init__(self): + dataset = xarray.open_dataset(self.opendap_url) - quiver_plot = axis.quiver( - lon, - lat, - self.dataset['u'].sel(time=time, method='nearest'), - self.dataset['v'].sel(time=time, method='nearest'), - units='width', - **kwargs, + super().__init__( + dataset=dataset, + uname='ugrdprs', + vname='vgrdprs', ) - axis.quiverkey( - quiver_plot, 0.9, 0.9, 1, r'$1 \frac{m}{s}$', labelpos='E', coordinates='figure' - ) - - return quiver_plot From 22b74ee7edefb8426cec281b0e023ebe9f2a4a74 Mon Sep 17 00:00:00 2001 From: Zach Date: Fri, 16 Apr 2021 09:43:04 -0400 Subject: [PATCH 5/6] gfs class --- packetraven/model.py | 161 +++++++++++++++++++++++-------------------- 1 file changed, 87 insertions(+), 74 deletions(-) diff --git a/packetraven/model.py b/packetraven/model.py index 4bf55183..bbd94301 100644 --- a/packetraven/model.py +++ b/packetraven/model.py @@ -1,5 +1,4 @@ -from datetime import datetime, timedelta - +from dateutil.parser import parse as parse_date from matplotlib import pyplot, quiver import numpy import pyproj @@ -25,16 +24,20 @@ class VectorField: Vector field of (u, v) values. """ - def __init__(self, time_deltas: [timedelta], projection: CRS = None): + def __init__(self, time_deltas: [timedelta], crs: CRS = None): """ Build vector field of (u, v) values. :param time_deltas: list of time deltas - :param projection: native projection of field + :param crs: native CRS of field """ - self.time_deltas = [numpy.timedelta64(time_delta) for time_delta in time_deltas] - self.projection = projection if projection is not None else WGS84 + if crs is None: + crs = WGS84 + if not isinstance(time_deltas, numpy.ndarray) or time_deltas.dtype != numpy.timedelta64: + time_deltas = numpy.array(time_deltas, dtype='timedelta64[s]') + self.time_deltas = time_deltas + self.crs = crs time_delta = numpy.nanmean(self.time_deltas).item() if type(time_delta) is int: @@ -75,7 +78,7 @@ def velocity(self, point: numpy.array, time: datetime) -> float: :return: magnitude of uv vector in m/s """ - return math.sqrt(self.u(point, time) ** 2 + self.v(point, time) ** 2) + return numpy.sqrt(self.u(point, time) ** 2 + self.v(point, time) ** 2) def direction(self, point: numpy.array, time: datetime) -> float: """ @@ -86,7 +89,7 @@ def direction(self, point: numpy.array, time: datetime) -> float: :return: radians from east of uv vector """ - return math.atan2(self.u(point, time), self.v(point, time)) + return numpy.arctan2(self.u(point, time), self.v(point, time)) def plot(self, time: datetime, axis: pyplot.Axes = None, **kwargs) -> quiver.Quiver: """ @@ -114,7 +117,7 @@ def __getitem__(self, position: (numpy.array, datetime)) -> numpy.array: return vector def __repr__(self) -> str: - return f'{self.__class__.__name__}({self.delta_t})' + return f'{self.__class__.__name__}({repr(self.delta_t)})' class RankineVortex(VectorField): @@ -123,7 +126,7 @@ class RankineVortex(VectorField): """ def __init__( - self, + self, center: (float, float), radius: float, period: timedelta, @@ -146,10 +149,10 @@ def __init__( super().__init__(time_deltas) def u(self, point: numpy.array, time: datetime) -> float: - return -self.velocity(point, time) * numpy.cos(numpy.atan2(*(point - self.center))) + return -self.velocity(point, time) * numpy.cos(numpy.arctan2(*(point - self.center))) def v(self, point: numpy.array, time: datetime) -> float: - return self.velocity(point, time) * numpy.sin(numpy.atan2(*(point - self.center))) + return self.velocity(point, time) * numpy.sin(numpy.arctan2(*(point - self.center))) def velocity(self, point: numpy.array, time: datetime) -> float: radial_distance = numpy.sqrt(numpy.sum((point - self.center) ** 2)) @@ -161,7 +164,7 @@ def velocity(self, point: numpy.array, time: datetime) -> float: def plot(self, axis: pyplot.Axes = None, **kwargs) -> quiver.Quiver: if axis is None: - axis = pyplot.axes(projection=cartopy.crs.PlateCarree()) + axis = pyplot.axes() points = [] radii = numpy.linspace(1, self.radius * 2, 20) @@ -201,13 +204,13 @@ class VectorDataset(VectorField): def __init__( self, dataset: xarray.Dataset, - u_name: str = 'u', - v_name: str = 'v', - x_name: str = 'lon', - y_name: str = 'lat', - z_name: str = 'z', - t_name: str = 'time', - coordinate_system: pyproj.Proj = None, + u_name: str = None, + v_name: str = None, + x_name: str = None, + y_name: str = None, + z_name: str = None, + t_name: str = None, + coordinate_system: CRS = None, ): """ Create new velocity field from given observation. @@ -222,54 +225,57 @@ def __init__( :param coordinate_system: coordinate system of observation """ - self.coordinate_system = ( - coordinate_system - if coordinate_system is not None - else WGS84 - ) - - variables_to_rename = { - u_name: 'u', - v_name: 'v', - x_name: 'x', - y_name: 'y', - t_name: 'time', - } - if z_name in dataset: - variables_to_rename[z_name] = 'z' - self.dataset = dataset.rename(variables_to_rename) - - x, y = pyproj.transform( - self.coordinate_system, - WEB_MERCATOR, - *numpy.meshgrid(self.dataset['x'].values, self.dataset['y'].values), - ) - - self.dataset['x'] = x[0, :] - self.dataset['y'] = y[:, 0] - - self.delta_x = numpy.mean(numpy.diff(self.dataset['x'])) - self.delta_y = numpy.mean(numpy.diff(self.dataset['y'])) - - super().__init__(numpy.diff(self.dataset['time'].values)) - - def _interpolate( - self, variable: str, point: numpy.array, time: datetime - ) -> xarray.DataArray: - transformed_point = pyproj.transform( - WEB_MERCATOR, self.coordinate_system, point[0], point[1] - ) - - x_name = f'{variable}_x' - y_name = f'{variable}_y' + if coordinate_system is None: + coordinate_system = WGS84 + + self.coordinate_system = coordinate_system + + if u_name is None: + u_name = 'u' + if v_name is None: + v_name = 'v' + if x_name is None: + x_name = 'lon' + if y_name is None: + y_name = 'lat' + if z_name is None: + z_name = 'z' + if t_name is None: + t_name = 'time' + + self.u_name = u_name + self.v_name = v_name + self.x_name = x_name + self.y_name = y_name + self.z_name = z_name + self.t_name = t_name + + self.dataset = dataset + + self.delta_x = numpy.mean(numpy.diff(self.dataset[self.x_name])) + self.delta_y = numpy.mean(numpy.diff(self.dataset[self.x_name])) + + super().__init__(numpy.diff(self.dataset[self.t_name].values)) + + def _interpolate(self, variable: str, point: numpy.array, time: datetime) -> xarray.DataArray: + if not isinstance(point, numpy.ndarray): + point = numpy.array(point) + if isinstance(time, str): + time = parse_date(time) + elif isinstance(time, timedelta): + time += self.dataset[self.t_name][0] x_range = slice( - self.dataset[x_name].sel({x_name: numpy.min(transformed_point[0]) - 1}, method='bfill').values.item(), - self.dataset[x_name].sel({x_name: numpy.max(transformed_point[0]) + 1}, method='ffill').values.item(), + self.dataset[self.x_name].sel({self.x_name: numpy.min(point[0]) - 1}, method='bfill').values.item(), + self.dataset[self.x_name].sel({self.x_name: numpy.max(point[0]) + 1}, method='ffill').values.item(), ) y_range = slice( - self.dataset[y_name].sel({y_name: numpy.min(transformed_point[1]) - 1}, method='bfill').values.item(), - self.dataset[y_name].sel({y_name: numpy.max(transformed_point[1]) + 1}, method='ffill').values.item(), + self.dataset[self.y_name].sel({ + self.y_name: numpy.min(point[1]) - 1, + }, method='bfill').values.item(), + self.dataset[self.y_name].sel({ + self.y_name: numpy.max(point[1]) + 1, + }, method='ffill').values.item(), ) time_range = slice( self.dataset['time'].sel(time=time, method='bfill').values, @@ -282,29 +288,35 @@ def _interpolate( cell = self.dataset[variable].sel( { 'time': time_range, - x_name: x_range, - y_name: y_range, + self.x_name: x_range, + self.y_name: y_range, } ) - if len(transformed_point.shape) > 1: + if len(point.shape) > 1: cell = cell.interp({'time': time}) if 'time' in cell.dims else cell return xarray.concat( [ - cell.interp({x_name: location[0], y_name: location[1]}) - for location in transformed_point.T + cell.interp({ + self.x_name: location[0], + self.y_name: location[1], + }) + for location in point.T ], dim='point', ) else: - cell = cell.interp({x_name: transformed_point[0], y_name: transformed_point[1]}) + cell = cell.interp({ + self.x_name: point[0], + self.y_name: point[1], + }) return cell.interp({'time': time}) if 'time' in cell.dims else cell def u(self, point: numpy.array, time: datetime) -> float: - return self._interpolate('u', point, time).values + return self._interpolate(self.u_name, point, time).values def v(self, point: numpy.array, time: datetime) -> float: - return self._interpolate('v', point, time).values + return self._interpolate(self.v_name, point, time).values class VectorGFS(VectorDataset): @@ -315,6 +327,7 @@ def __init__(self): super().__init__( dataset=dataset, - uname='ugrdprs', - vname='vgrdprs', + u_name='ugrdprs', + v_name='vgrdprs', + z_name='lev', ) From e3c32b3ffb836ab717145e59eec2ff15195ab3c9 Mon Sep 17 00:00:00 2001 From: Lint Action Date: Sun, 13 Jun 2021 15:21:25 +0000 Subject: [PATCH 6/6] Fix code style issues with oitnb --- packetraven/model.py | 58 +++++++++++++++++++++----------------------- 1 file changed, 27 insertions(+), 31 deletions(-) diff --git a/packetraven/model.py b/packetraven/model.py index bbd94301..3f6677d9 100644 --- a/packetraven/model.py +++ b/packetraven/model.py @@ -34,7 +34,10 @@ def __init__(self, time_deltas: [timedelta], crs: CRS = None): if crs is None: crs = WGS84 - if not isinstance(time_deltas, numpy.ndarray) or time_deltas.dtype != numpy.timedelta64: + if ( + not isinstance(time_deltas, numpy.ndarray) + or time_deltas.dtype != numpy.timedelta64 + ): time_deltas = numpy.array(time_deltas, dtype='timedelta64[s]') self.time_deltas = time_deltas self.crs = crs @@ -184,9 +187,7 @@ def plot(self, axis: pyplot.Axes = None, **kwargs) -> quiver.Quiver: ) vectors = [self[point, datetime.now()] for point in points] - points = list( - zip(*pyproj.transform(WEB_MERCATOR, WGS84, *zip(*points))) - ) + points = list(zip(*pyproj.transform(WEB_MERCATOR, WGS84, *zip(*points)))) quiver_plot = axis.quiver(*zip(*points), *zip(*vectors), units='width', **kwargs) axis.quiverkey( @@ -257,7 +258,9 @@ def __init__( super().__init__(numpy.diff(self.dataset[self.t_name].values)) - def _interpolate(self, variable: str, point: numpy.array, time: datetime) -> xarray.DataArray: + def _interpolate( + self, variable: str, point: numpy.array, time: datetime + ) -> xarray.DataArray: if not isinstance(point, numpy.ndarray): point = numpy.array(point) if isinstance(time, str): @@ -266,16 +269,20 @@ def _interpolate(self, variable: str, point: numpy.array, time: datetime) -> xar time += self.dataset[self.t_name][0] x_range = slice( - self.dataset[self.x_name].sel({self.x_name: numpy.min(point[0]) - 1}, method='bfill').values.item(), - self.dataset[self.x_name].sel({self.x_name: numpy.max(point[0]) + 1}, method='ffill').values.item(), + self.dataset[self.x_name] + .sel({self.x_name: numpy.min(point[0]) - 1}, method='bfill') + .values.item(), + self.dataset[self.x_name] + .sel({self.x_name: numpy.max(point[0]) + 1}, method='ffill') + .values.item(), ) y_range = slice( - self.dataset[self.y_name].sel({ - self.y_name: numpy.min(point[1]) - 1, - }, method='bfill').values.item(), - self.dataset[self.y_name].sel({ - self.y_name: numpy.max(point[1]) + 1, - }, method='ffill').values.item(), + self.dataset[self.y_name] + .sel({self.y_name: numpy.min(point[1]) - 1,}, method='bfill') + .values.item(), + self.dataset[self.y_name] + .sel({self.y_name: numpy.max(point[1]) + 1,}, method='ffill') + .values.item(), ) time_range = slice( self.dataset['time'].sel(time=time, method='bfill').values, @@ -286,30 +293,20 @@ def _interpolate(self, variable: str, point: numpy.array, time: datetime) -> xar time_range = time_range.start cell = self.dataset[variable].sel( - { - 'time': time_range, - self.x_name: x_range, - self.y_name: y_range, - } + {'time': time_range, self.x_name: x_range, self.y_name: y_range,} ) if len(point.shape) > 1: cell = cell.interp({'time': time}) if 'time' in cell.dims else cell return xarray.concat( [ - cell.interp({ - self.x_name: location[0], - self.y_name: location[1], - }) + cell.interp({self.x_name: location[0], self.y_name: location[1],}) for location in point.T ], dim='point', ) else: - cell = cell.interp({ - self.x_name: point[0], - self.y_name: point[1], - }) + cell = cell.interp({self.x_name: point[0], self.y_name: point[1],}) return cell.interp({'time': time}) if 'time' in cell.dims else cell def u(self, point: numpy.array, time: datetime) -> float: @@ -320,14 +317,13 @@ def v(self, point: numpy.array, time: datetime) -> float: class VectorGFS(VectorDataset): - opendap_url = 'http://nomads.ncep.noaa.gov:80/dods/gfs_0p25_1hr/gfs20210413/gfs_0p25_1hr_12z' + opendap_url = ( + 'http://nomads.ncep.noaa.gov:80/dods/gfs_0p25_1hr/gfs20210413/gfs_0p25_1hr_12z' + ) def __init__(self): dataset = xarray.open_dataset(self.opendap_url) super().__init__( - dataset=dataset, - u_name='ugrdprs', - v_name='vgrdprs', - z_name='lev', + dataset=dataset, u_name='ugrdprs', v_name='vgrdprs', z_name='lev', )