Skip to content

Commit

Permalink
read fort.22 from file (#78)
Browse files Browse the repository at this point in the history 
* construct `fort.22` with `.join()`

* read `fort.22`

* add class method as boilerplate for inherited implementations
  • Loading branch information
@zacharyburnett
zacharyburnett authored Apr 19, 2021
1 parent 11b040b commit e1b33a2
Show file tree
Hide file tree
Showing 5 changed files with 373 additions and 41 deletions.
6 changes: 6 additions & 0 deletions adcircpy/forcing/winds/base.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,3 +12,9 @@ def __init__(self, nws: int, interval_seconds: int):
@abstractmethod
def write(self, directory: PathLike, overwrite: bool = False):
raise NotImplementedError

@classmethod
def from_fort22(cls, fort22: PathLike, nws: int = None) -> 'WindForcing':
raise NotImplementedError(
f'reading `fort.22` is not implemented for {cls}'
)
211 changes: 170 additions & 41 deletions adcircpy/forcing/winds/best_track.py
View file
Original file line number Diff line number Diff line change
@@ -1,3 +1,4 @@
from collections import Collection
from datetime import datetime, timedelta
from functools import wraps
import gzip
Expand All @@ -7,14 +8,16 @@
from os import PathLike
import pathlib
import time
from typing import Any
import urllib.request

from haversine import haversine
import matplotlib.pyplot as plt
from matplotlib.transforms import Bbox
import numpy as np
import pandas
from pandas import DataFrame, read_csv
from pyproj import Proj, CRS, Transformer
from pyproj import CRS, Proj, Transformer
from shapely import ops
from shapely.geometry import Point, Polygon
import utm
Expand All @@ -35,61 +38,82 @@ def __init__(self, storm_id, nws: int = 20,

def __str__(self):
record_number = self._generate_record_numbers()
fort22 = ''
fort22 = []
for i, (_, row) in enumerate(self.df.iterrows()):
fort22 += f'{row["basin"]:<2},{row["storm_number"]:>3},' \
f' {row["datetime"]:%Y%m%d%H},' \
f'{"":3},{row["record_type"]:>5},' \
f'{int((row["datetime"] - self.start_date) / timedelta(hours=1)):>4},'
line = []

line.extend([
f'{row["basin"]:<2}',
f'{row["storm_number"]:>3}',
f'{row["datetime"]:%Y%m%d%H}',
f'{"":3}',
f'{row["record_type"]:>5}',
f'{convert_value((row["datetime"] - self.start_date) / timedelta(hours=1), int):>4}',
])

if row["latitude"] >= 0:
fort22 += f'{int(row["latitude"] / .1):>4}N,'
line.append(f'{convert_value(row["latitude"] / .1, int):>4}N')
else:
fort22 += f'{int(row["latitude"] / -.1):>4}S,'
line.append(f'{convert_value(row["latitude"] / -.1, int):>4}S')
if row["longitude"] >= 0:
fort22 += f'{int(row["longitude"] / .1):>5}E,'
line.append(f'{convert_value(row["longitude"] / .1, int):>5}E')
else:
fort22 += f'{int(row["longitude"] / -.1):>5}W,'
fort22 += f'{int(row["max_sustained_wind_speed"]):>4},' \
f'{int(row["central_pressure"]):>5},' \
f'{row["development_level"]:>3},' \
f'{int(row["isotach"]):>4},' \
f'{row["quadrant"]:>4},' \
f'{int(row["radius_for_NEQ"]):>5},' \
f'{int(row["radius_for_SEQ"]):>5},' \
f'{int(row["radius_for_SWQ"]):>5},' \
f'{int(row["radius_for_NWQ"]):>5},'
line.append(
f'{convert_value(row["longitude"] / -.1, int):>5}W')

line.extend([
f'{convert_value(row["max_sustained_wind_speed"], int):>4}',
f'{convert_value(row["central_pressure"], int):>5}',
f'{row["development_level"]:>3}',
f'{convert_value(row["isotach"], int):>4}',
f'{row["quadrant"]:>4}',
f'{convert_value(row["radius_for_NEQ"], int):>5}',
f'{convert_value(row["radius_for_SEQ"], int):>5}',
f'{convert_value(row["radius_for_SWQ"], int):>5}',
f'{convert_value(row["radius_for_NWQ"], int):>5}',
])

if row["background_pressure"] is None:
row["background_pressure"] = \
self.df["background_pressure"].iloc[i - 1]
if (row["background_pressure"] <= row["central_pressure"]
and 1013 > row["central_pressure"]):
fort22 += f'{1013:>5},'
background_pressure = 1013
elif (row["background_pressure"] <= row["central_pressure"]
and 1013 <= row["central_pressure"]):
fort22 += f'{int(row["central_pressure"] + 1):>5},'
background_pressure = convert_value(
row["central_pressure"] + 1, int)
else:
fort22 += f'{int(row["background_pressure"]):>5},'
fort22 += f'{int(row["radius_of_last_closed_isobar"]):>5},' \
f'{int(row["radius_of_maximum_winds"]):>4},'
fort22 += f'{"":>5},' # gust
fort22 += f'{"":>4},' # eye
fort22 += f'{"":>4},' # subregion
fort22 += f'{"":>4},' # maxseas
fort22 += f'{"":>4},' # initials
fort22 += f'{row["direction"]:>3},' \
f'{row["speed"]:>4},' \
f'{row["name"]:^12},'
background_pressure = convert_value(row["background_pressure"],
int)
line.append(f'{background_pressure:>5}')

line.extend([
f'{convert_value(row["radius_of_last_closed_isobar"], int):>5}',
f'{convert_value(row["radius_of_maximum_winds"], int):>4}',
f'{"":>5}', # gust
f'{"":>4}', # eye
f'{"":>4}', # subregion
f'{"":>4}', # maxseas
f'{"":>4}', # initials
f'{row["direction"]:>3}',
f'{row["speed"]:>4}',
f'{row["name"]:^12}',
])

# from this point forwards it's all aswip
fort22 += f'{record_number[i]:>4},' \
f'\n'
return fort22
line.append(f'{record_number[i]:>4}')

fort22.append(','.join(line))

return '\n'.join(fort22)

def write(self, path: PathLike, overwrite: bool = False):
if not isinstance(path, pathlib.Path):
path = pathlib.Path(path)
if path.exists() and overwrite is False:
raise Exception(
'File exist, set overwrite=True to allow overwrite.')
'File exist, set overwrite=True to allow overwrite.')
with open(path, 'w') as f:
f.write(str(self))

Expand Down Expand Up @@ -404,14 +428,14 @@ def clip_to_bbox(self, bbox, bbox_crs):
ellps={
'GRS 1980': 'GRS80',
'WGS 84': 'WGS84'
}[df_crs.ellipsoid.name]
)
}[df_crs.ellipsoid.name]
)
transformer = Transformer.from_crs(
df_crs, utm_crs, always_xy=True)
df_crs, utm_crs, always_xy=True)
p = Point(*transformer.transform(lon, lat))
pol = p.buffer(radii)
transformer = Transformer.from_crs(
utm_crs, bbox_crs, always_xy=True)
utm_crs, bbox_crs, always_xy=True)
pol = ops.transform(transformer.transform, pol)
if _switch is True:
if not pol.intersects(bbox_pol):
Expand All @@ -431,7 +455,7 @@ def clip_to_bbox(self, bbox, bbox_crs):

if _found_start_date is False:
raise Exception(
f'No data within mesh bounding box for storm {self.storm_id}.')
f'No data within mesh bounding box for storm {self.storm_id}.')

def plot_track(self, axes=None, show=False, color='k', **kwargs):
kwargs.update({'color': color})
Expand Down Expand Up @@ -536,6 +560,111 @@ def _compute_velocity(data):
int(np.around(bearing, 0)))
return data

@classmethod
def from_fort22(
cls,
fort22: PathLike,
nws: int = None,
) -> 'BestTrackForcing':
if nws is None:
nws = 20

try:
with open(fort22) as fort22_file:
fort22 = fort22_file.readlines()
except:
fort22 = str(fort22).splitlines()

data = {
'basin': [],
'storm_number': [],
'datetime': [],
'record_type': [],
'latitude': [],
'longitude': [],
'max_sustained_wind_speed': [],
'central_pressure': [],
'development_level': [],
'isotach': [],
'quadrant': [],
'radius_for_NEQ': [],
'radius_for_SEQ': [],
'radius_for_SWQ': [],
'radius_for_NWQ': [],
'background_pressure': [],
'radius_of_last_closed_isobar': [],
'radius_of_maximum_winds': [],
'name': [],
'direction': [],
'speed': [],
}

for index, row in enumerate(fort22):
row = [value.strip() for value in row.split(',')]

row_data = {key: None for key in data}

row_data['basin'] = row[0]
row_data['storm_number'] = row[1]
row_data['datetime'] = datetime.strptime(row[2], '%Y%m%d%H')
row_data['record_type'] = row[4]

latitude = row[6]
if 'N' in latitude:
latitude = float(latitude[:-1]) * 0.1
elif 'S' in latitude:
latitude = float(latitude[:-1]) * -0.1
row_data['latitude'] = latitude

longitude = row[7]
if 'E' in longitude:
longitude = float(longitude[:-1]) * 0.1
elif 'W' in longitude:
longitude = float(longitude[:-1]) * -0.1
row_data['longitude'] = longitude

row_data['max_sustained_wind_speed'] = convert_value(row[8], int)
row_data['central_pressure'] = convert_value(row[9], int)
row_data['development_level'] = row[10]
row_data['isotach'] = convert_value(row[11], int)
row_data['quadrant'] = row[12]
row_data['radius_for_NEQ'] = convert_value(row[13], int)
row_data['radius_for_SEQ'] = convert_value(row[14], int)
row_data['radius_for_SWQ'] = convert_value(row[15], int)
row_data['radius_for_NWQ'] = convert_value(row[16], int)
row_data['background_pressure'] = convert_value(row[17], int)
row_data['radius_of_last_closed_isobar'] = convert_value(row[18],
int)
row_data['radius_of_maximum_winds'] = convert_value(row[19], int)
row_data['direction'] = row[25]
row_data['speed'] = row[26]
row_data['name'] = row[27]

for key, value in row_data.items():
if isinstance(data[key], Collection):
data[key].append(value)
elif data[key] is None:
data[key] = value

storm_id = f'{data["name"][0]}{data["datetime"][0]:%Y}'

instance = cls(
storm_id=storm_id,
nws=nws,
start_date=min(data['datetime']),
end_date=max(data['datetime']),
)

instance.__df = pandas.DataFrame(data=data)

return instance


def convert_value(value: Any, to_type: type) -> Any:
if value is not None and value != '':
value = to_type(value)
return value


def retry(ExceptionToCheck, tries=4, delay=3, backoff=2, logger=None):
"""Retry calling the decorated function using an exponential backoff.
Expand Down
Loading

0 comments on commit e1b33a2

Please sign in to comment.