From 001445de82ecf89648a716944574305c71c14b85 Mon Sep 17 00:00:00 2001 From: cuill Date: Wed, 28 Jul 2021 09:56:09 -0400 Subject: [PATCH] Merged NWM hindcast inventory class into nwm --- .idea/.gitignore | 3 + .../inspectionProfiles/profiles_settings.xml | 6 + .idea/modules.xml | 8 ++ .idea/pyschism.iml | 12 ++ .idea/vcs.xml | 6 + pyschism/forcing/hycom/hycom.py | 1 - pyschism/forcing/hydrology/nwm.py | 129 ++++++++++++++++++ pyschism/mesh/vgrid.py | 4 +- pyschism/outputs/outputs.py | 38 +++++- pyschism/outputs/plot.py | 2 +- 10 files changed, 204 insertions(+), 5 deletions(-) create mode 100644 .idea/.gitignore create mode 100644 .idea/inspectionProfiles/profiles_settings.xml create mode 100644 .idea/modules.xml create mode 100644 .idea/pyschism.iml create mode 100644 .idea/vcs.xml diff --git a/.idea/.gitignore b/.idea/.gitignore new file mode 100644 index 00000000..26d33521 --- /dev/null +++ b/.idea/.gitignore @@ -0,0 +1,3 @@ +# Default ignored files +/shelf/ +/workspace.xml diff --git a/.idea/inspectionProfiles/profiles_settings.xml b/.idea/inspectionProfiles/profiles_settings.xml new file mode 100644 index 00000000..105ce2da --- /dev/null +++ b/.idea/inspectionProfiles/profiles_settings.xml @@ -0,0 +1,6 @@ + + + + \ No newline at end of file diff --git a/.idea/modules.xml b/.idea/modules.xml new file mode 100644 index 00000000..7053ba80 --- /dev/null +++ b/.idea/modules.xml @@ -0,0 +1,8 @@ + + + + + + + + \ No newline at end of file diff --git a/.idea/pyschism.iml b/.idea/pyschism.iml new file mode 100644 index 00000000..8b8c3954 --- /dev/null +++ b/.idea/pyschism.iml @@ -0,0 +1,12 @@ + + + + + + + + + + \ No newline at end of file diff --git a/.idea/vcs.xml b/.idea/vcs.xml new file mode 100644 index 00000000..94a25f7f --- /dev/null +++ b/.idea/vcs.xml @@ -0,0 +1,6 @@ + + + + + + \ No newline at end of file diff --git a/pyschism/forcing/hycom/hycom.py b/pyschism/forcing/hycom/hycom.py index 687a276e..7d05a51e 100644 --- a/pyschism/forcing/hycom/hycom.py +++ b/pyschism/forcing/hycom/hycom.py @@ -14,7 +14,6 @@ from matplotlib.transforms import Bbox from pyschism.mesh import Hgrid -from pyschism.dates import localize_datetime, nearest_cycle_date, pivot_time logger = logging.getLogger(__name__) diff --git a/pyschism/forcing/hydrology/nwm.py b/pyschism/forcing/hydrology/nwm.py index f169adae..6f370a8d 100644 --- a/pyschism/forcing/hydrology/nwm.py +++ b/pyschism/forcing/hydrology/nwm.py @@ -14,6 +14,7 @@ import boto3 from botocore import UNSIGNED from botocore.config import Config +import pandas as pd import geopandas as gpd from netCDF4 import Dataset import numpy as np @@ -110,6 +111,7 @@ def __init__(self, hgrid): 'No National Water model intersections found on the mesh.') intersection = gpd.GeoDataFrame(data, crs=hgrid.crs) del data + print(intersection) # 2) Generate element centroid KDTree centroids = [] @@ -422,6 +424,133 @@ def requested_product(self): }[self.product] +class AWSHindcastInventory: + + def __init__( + self, + start_date: datetime = None, + rnday: Union[int, float, timedelta] = timedelta(days=5.), + product='CHRTOUT_DOMAIN1.comp', + verbose=False, + fallback=True, + ): + """This will download the National Water Model retro data. + A 26-year (January 1993 through December 2018) retrospective + simulation using version 2.0 of the NWM. + + NetCDF files are saved to the system's temporary directory. + """ + self.start_date = dates.nearest_cycle() if start_date is None \ + else dates.nearest_cycle(dates.localize_datetime(start_date)) + # self.start_date = self.start_date.replace(tzinfo=None) + self.rnday = rnday if isinstance(rnday, timedelta) \ + else timedelta(days=rnday) + self.product = product + self.fallback = fallback + self._files = {_: None for _ in np.arange( + self.start_date, + self.start_date + self.rnday + self.output_interval, + self.output_interval + ).astype(datetime)} + + paginator=self.s3.get_paginator('list_objects_v2') + pages=paginator.paginate(Bucket=self.bucket, + Prefix=f'full_physics/{self.start_date.year}') + + self.data=[] + for page in pages: + for obj in page['Contents']: + self.data.append(obj) + + self.file_metadata = list(sorted([ + _['Key'] for _ in self.data if 'CHRTOUT_DOMAIN1.comp' in _['Key'] + ])) + + timevector=np.arange(datetime(self.start_date.year, 1, 1), + datetime(self.start_date.year+1,1, 1), + np.timedelta64(1, 'h'), + dtype='datetime64') + + timefile= {pd.to_datetime(str(timevector[i])): self.file_metadata[i] + for i in range(len(timevector))} + + for requested_time, _ in self._files.items(): + #print(f'Requesting NWM data for time {requested_time}') + key=timefile.get(requested_time) + print(f'Requesting NWM data for time {requested_time}, {key}') + logger.info(f'Requesting NWM data for time {requested_time}') + self._files[requested_time] = self.request_data(key, requested_time) + + def request_data(self, key, request_time): + + #print(key) + filename = pathlib.Path(self.tmpdir.name) / key + filename.parent.mkdir(parents=True, exist_ok=True) + + with open(filename, 'wb') as f: + #logger.info(f'Downloading file {key}, ') + print(f'Downloading file {key}, ') + self.s3.download_fileobj(self.bucket, key, f) + return filename + + def get_nc_pairing_indexes(self, pairings: NWMElementPairings): + nc_feature_id = Dataset(self.files[0])['feature_id'][:] + + def get_aggregated_features(features): + aggregated_features = [] + for source_feats in features: + aggregated_features.extend(list(source_feats)) + in_file = np.where( + np.in1d(nc_feature_id, aggregated_features, + assume_unique=True))[0] + in_file_2 = [] + sidx = 0 + for source_feats in features: + eidx = sidx + len(source_feats) + in_file_2.append(in_file[sidx:eidx].tolist()) + sidx = eidx + return in_file_2 + + sources = get_aggregated_features(pairings.sources.values()) + sinks = get_aggregated_features(pairings.sinks.values()) + return sources, sinks + + @property + def bucket(self): + return 'noaa-nwm-retro-v2.0-pds' + + @property + def nearest_cycle(self) -> datetime: + return dates.nearest_cycle(self.start_date) + + @property + def output_interval(self) -> timedelta: + return { + 'CHRTOUT_DOMAIN1.comp': timedelta(hours=1) + }[self.product] + + @property + def s3(self): + try: + return self._s3 + except AttributeError: + self._s3 = boto3.client( + 's3', config=Config(signature_version=UNSIGNED)) + return self._s3 + + @property + def tmpdir(self): + try: + return self._tmpdir + except AttributeError: + self._tmpdir = tempfile.TemporaryDirectory() + return self._tmpdir + + @property + def files(self): + return sorted(list(pathlib.Path(self.tmpdir.name).glob('**/*.comp'))) + + class NationalWaterModel(Hydrology): def __init__(self, aggregation_radius=None): diff --git a/pyschism/mesh/vgrid.py b/pyschism/mesh/vgrid.py index 5f756d63..f7b3cf49 100644 --- a/pyschism/mesh/vgrid.py +++ b/pyschism/mesh/vgrid.py @@ -14,7 +14,7 @@ class Vgrid: def __init__(self): """Represents a SCHISM vertical grid.""" - #self._vgrid = self._get_2D_string() + self._vgrid = self._get_2D_string() pass @classmethod @@ -63,7 +63,7 @@ def read_vgrid(fname): lines=lines[2:] kbp=np.array([int(i.split()[1])-1 for i in lines]) NP=len(kbp) - print(NP) + #print(NP) sigma=-np.ones([NP,nvrt]) for i, line in enumerate(lines): sigma[i,kbp[i]:]=np.array(line.strip().split()[2:]).astype('float') diff --git a/pyschism/outputs/outputs.py b/pyschism/outputs/outputs.py index cda65c83..42e01573 100644 --- a/pyschism/outputs/outputs.py +++ b/pyschism/outputs/outputs.py @@ -28,7 +28,6 @@ ) from pyschism.mesh.base import Gr3 - def get_stack_id_by_datetime(dt: datetime, flattened_timevector, stacks) -> str: output_stack = None for stack_id, stack_data in stacks.items(): @@ -188,6 +187,43 @@ def animate(index): return anim + def get_zcor_interp_coefficient(self, z0, zcor, level): + ''' + z0: water elevation (elev) + zcor: schout outputs zcor[np,nvrt] for each time record. + level: fixed depth (< 0, e.g., 4.5 m below the surface, level=[-4.5]) where values is interplated at. + ''' + NP=zcor.shape[0] + nvrt=zcor.shape[1] + k1=np.full((NP), np.nan) + coeff=np.full((NP), np.nan) + zinter=np.ones(NP)*level+z0 + + #calculate kbp + + + #surface + idxs=zinter>=zcor[:,-1] + k1[idxs]=nvrt-2 + coeff[idxs]=1.0 + + #bottom + idxs=zinter=zcor[:,k])*(zinter