diff --git a/pyschism/forcing/hycom/hycom2schism.py b/pyschism/forcing/hycom/hycom2schism.py index 768a38aa..bceec245 100644 --- a/pyschism/forcing/hycom/hycom2schism.py +++ b/pyschism/forcing/hycom/hycom2schism.py @@ -23,6 +23,20 @@ logger = logging.getLogger(__name__) +def convert_longitude(ds): + lon_name = 'lon' + ds['_lon_adjusted'] = xr.where( + ds[lon_name] > 180, + ds[lon_name] - 360, + ds[lon_name]) + ds = ( + ds.swap_dims({lon_name: '_lon_adjusted'}) + .sel(**{'_lon_adjusted': sorted(ds._lon_adjusted)}) + .drop(lon_name) + ) + ds = ds.rename({'_lon_adjusted': lon_name}) + return ds + def get_database(date, Bbox=None): if date >= datetime(2018, 12, 4): database = f'GLBy0.08/expt_93.0' @@ -415,7 +429,6 @@ def fetch_data(self, outdir: Union[str, os.PathLike], start_date, rnday, elev2D= #timeseries_uv[it,:,:,1]=vvel_int ds.close() - logger.info(f'Writing *th.nc takes {time()-t0} seconds') class Nudge: @@ -671,7 +684,11 @@ class DownloadHycom: def __init__(self, hgrid): - self.bbox = hgrid.bbox + xmin, xmax = hgrid.coords[:, 0].min(), hgrid.coords[:, 0].max() + ymin, ymax = hgrid.coords[:, 1].min(), hgrid.coords[:, 1].max() + xmin = xmin + 360. if xmin < 0 else xmin + xmax = xmax + 360. if xmax < 0 else xmax + self.bbox = Bbox.from_extents(xmin, ymin, xmax, ymax) def fetch_data(self, date): @@ -686,16 +703,13 @@ def fetch_data(self, date): foutname = f'SSH_{date.strftime("%Y%m%d")}.nc' logger.info(f'filename is {foutname}') ds = xr.open_dataset(url_ssh) - ds1 = ds.rename_dims({'lon':'xlon'}) - ds2 = ds1.rename_dims({'lat':'ylat'}) - ds3 = ds2.rename_vars({'lat':'ylat'}) - ds4 = ds3.rename_vars({'lon':'xlon'}) - ds4.to_netcdf(foutname, 'w', 'NETCDF3_CLASSIC', unlimited_dims='time') + ds = convert_longitude(ds) + ds = ds.rename_dims({'lon':'xlon'}) + ds = ds.rename_dims({'lat':'ylat'}) + ds = ds.rename_vars({'lat':'ylat'}) + ds = ds.rename_vars({'lon':'xlon'}) + ds.to_netcdf(foutname, 'w', 'NETCDF3_CLASSIC', unlimited_dims='time') ds.close() - ds1.close() - ds2.close() - ds3.close() - ds4.close() url_uv = f'https://tds.hycom.org/thredds/dodsC/{database}?lat[{lat_idx1}:1:{lat_idx2}],' + \ f'lon[{lon_idx1}:1:{lon_idx2}],depth[0:1:-1],time[{time_idx}],' + \ @@ -705,16 +719,13 @@ def fetch_data(self, date): foutname = f'UV_{date.strftime("%Y%m%d")}.nc' logger.info(f'filename is {foutname}') ds = xr.open_dataset(url_uv) - ds1 = ds.rename_dims({'lon':'xlon'}) - ds2 = ds1.rename_dims({'lat':'ylat'}) - ds3 = ds2.rename_vars({'lat':'ylat'}) - ds4 = ds3.rename_vars({'lon':'xlon'}) - ds4.to_netcdf(foutname, 'w', 'NETCDF3_CLASSIC', unlimited_dims='time') + ds = convert_longitude(ds) + ds = ds.rename_dims({'lon':'xlon'}) + ds = ds.rename_dims({'lat':'ylat'}) + ds = ds.rename_vars({'lat':'ylat'}) + ds = ds.rename_vars({'lon':'xlon'}) + ds.to_netcdf(foutname, 'w', 'NETCDF3_CLASSIC', unlimited_dims='time') ds.close() - ds1.close() - ds2.close() - ds3.close() - ds4.close() url_ts = f'https://tds.hycom.org/thredds/dodsC/{database}?lat[{lat_idx1}:1:{lat_idx2}],' + \ f'lon[{lon_idx1}:1:{lon_idx2}],depth[0:1:-1],time[{time_idx}],' + \ @@ -725,6 +736,7 @@ def fetch_data(self, date): logger.info(f'filename is {foutname}') ds = xr.open_dataset(url_ts) + temp = ds.water_temp.values salt = ds.salinity.values dep = ds.depth.values @@ -733,9 +745,9 @@ def fetch_data(self, date): ptemp = ConvertTemp(salt, temp, dep) #drop water_temp variable and add new temperature variable - ds1 = ds.drop('water_temp') - ds1['temperature']=(['time','depth','lat','lon'], ptemp) - ds1.temperature.attrs = { + ds = ds.drop('water_temp') + ds['temperature']=(['time','depth','lat','lon'], ptemp) + ds.temperature.attrs = { 'long_name': 'Sea water potential temperature', 'standard_name': 'sea_water_potential_temperature', 'units': 'degC' @@ -744,14 +756,10 @@ def fetch_data(self, date): #ds.assign(water_temp2=ptemp) #ds.assign.attrs = ds.water_temp.attrs - ds2 = ds1.rename_dims({'lon':'xlon'}) - ds3 = ds2.rename_dims({'lat':'ylat'}) - ds4 = ds3.rename_vars({'lat':'ylat'}) - ds5 = ds4.rename_vars({'lon':'xlon'}) - ds5.to_netcdf(foutname, 'w', unlimited_dims='time', encoding={'temperature':{'dtype': 'h', '_FillValue': -30000.,'scale_factor': 0.001, 'add_offset': 20., 'missing_value': -30000.}}) + ds = convert_longitude(ds) + ds = ds.rename_dims({'lon':'xlon'}) + ds = ds.rename_dims({'lat':'ylat'}) + ds = ds.rename_vars({'lat':'ylat'}) + ds = ds.rename_vars({'lon':'xlon'}) + ds.to_netcdf(foutname, 'w', unlimited_dims='time', encoding={'temperature':{'dtype': 'h', '_FillValue': -30000.,'scale_factor': 0.001, 'add_offset': 20., 'missing_value': -30000.}}) ds.close() - ds1.close() - ds2.close() - ds3.close() - ds4.close() - ds5.close() diff --git a/pyschism/forcing/nws/nws2/era5.py b/pyschism/forcing/nws/nws2/era5.py index f899059e..7bd77217 100644 --- a/pyschism/forcing/nws/nws2/era5.py +++ b/pyschism/forcing/nws/nws2/era5.py @@ -89,10 +89,10 @@ def xy_grid(self): lon.append(x) return np.meshgrid(np.array(lon), self.lat) -def put_sflux_fields(iday, date, timevector, ds, nx_grid, ny_grid, air, rad, prc, OUTDIR): +def put_sflux_fields(iday, date, timevector, ds, nx_grid, ny_grid, air, rad, prc, output_interval, OUTDIR): rt=pd.to_datetime(str(date)) idx=np.where(rt == timevector)[0].item() - times=[i/24 for i in np.arange(0, 25, 3)] + times=[i/24 for i in np.arange(0, 25, output_interval)] if air is True: with Dataset(OUTDIR / "sflux_air_1.{:04}.nc".format(iday+1), 'w', format='NETCDF3_CLASSIC') as dst: @@ -128,7 +128,7 @@ def put_sflux_fields(iday, date, timevector, ds, nx_grid, ny_grid, air, rad, prc dst['prmsl'].long_name = "Pressure reduced to MSL" dst['prmsl'].standard_name = "air_pressure_at_sea_level" dst['prmsl'].units = "Pa" - dst['prmsl'][:,:,:]=ds['msl'][idx:idx+26:3,::-1,:] + dst['prmsl'][:,:,:]=ds['msl'][idx:idx+25:output_interval, ::-1, :] # spfh dst.createVariable('spfh', 'f4', ('time', 'ny_grid', 'nx_grid')) @@ -137,8 +137,8 @@ def put_sflux_fields(iday, date, timevector, ds, nx_grid, ny_grid, air, rad, prc dst['spfh'].standard_name = "specific_humidity" dst['spfh'].units = "1" #convert dewpoint to specific humidity - d2m = ds['d2m'][idx:idx+26:3,::-1,:] - msl = ds['msl'][idx:idx+26:3,::-1,:] + d2m = ds['d2m'][idx:idx+25:output_interval, ::-1, :] + msl = ds['msl'][idx:idx+25:output_interval, ::-1, :] Td = d2m - 273.15 e1 = 6.112*np.exp((17.67*Td)/(Td + 243.5)) spfh = (0.622*e1)/(msl*0.01 - (0.378*e1)) @@ -149,7 +149,7 @@ def put_sflux_fields(iday, date, timevector, ds, nx_grid, ny_grid, air, rad, prc dst['stmp'].long_name = "Surface Air Temperature (2m AGL)" dst['stmp'].standard_name = "air_temperature" dst['stmp'].units = "K" - dst['stmp'][:,:,:]=ds['t2m'][idx:idx+26:3,::-1,:] + dst['stmp'][:,:,:]=ds['t2m'][idx:idx+25:output_interval, ::-1, :] # uwind dst.createVariable('uwind', 'f4', ('time', 'ny_grid', 'nx_grid')) @@ -157,7 +157,7 @@ def put_sflux_fields(iday, date, timevector, ds, nx_grid, ny_grid, air, rad, prc "(10m AGL)" dst['uwind'].standard_name = "eastward_wind" dst['uwind'].units = "m/s" - dst['uwind'][:,:,:]=ds['u10'][idx:idx+26:3,::-1,:] + dst['uwind'][:,:,:]=ds['u10'][idx:idx+25:output_interval, ::-1, :] # vwind dst.createVariable('vwind', 'f4', ('time', 'ny_grid', 'nx_grid')) @@ -165,7 +165,7 @@ def put_sflux_fields(iday, date, timevector, ds, nx_grid, ny_grid, air, rad, prc "(10m AGL)" dst['vwind'].standard_name = "northward_wind" dst['vwind'].units = "m/s" - dst['vwind'][:,:,:]=ds['v10'][idx:idx+26:3,::-1,:] + dst['vwind'][:,:,:]=ds['v10'][idx:idx+25:output_interval, ::-1, :] if prc is True: with Dataset(OUTDIR / "sflux_prc_1.{:04}.nc".format(iday+1), 'w', format='NETCDF3_CLASSIC') as dst: @@ -201,7 +201,7 @@ def put_sflux_fields(iday, date, timevector, ds, nx_grid, ny_grid, air, rad, prc dst['prate'].long_name = "Surface Precipitation Rate" dst['prate'].standard_name = "air_pressure_at_sea_level" dst['prate'].units = "kg/m^2/s" - dst['prate'][:,:,:]=ds['mtpr'][idx:idx+26:3,::-1,:] + dst['prate'][:,:,:]=ds['mtpr'][idx:idx+25:output_interval, ::-1, :] if rad is True: with Dataset(OUTDIR / "sflux_rad_1.{:04}.nc".format(iday+1), 'w', format='NETCDF3_CLASSIC') as dst: @@ -237,14 +237,14 @@ def put_sflux_fields(iday, date, timevector, ds, nx_grid, ny_grid, air, rad, prc dst['dlwrf'].long_name = "Downward Long Wave Radiation Flux" dst['dlwrf'].standard_name = "surface_downwelling_longwave_flux_in_air" dst['dlwrf'].units = "W/m^2" - dst['dlwrf'][:,:,:]=ds['msdwlwrf'][idx:idx+26:3,::-1,:] + dst['dlwrf'][:,:,:]=ds['msdwlwrf'][idx:idx+25:output_interval, ::-1, :] # dwrf dst.createVariable('dswrf', 'f4', ('time', 'ny_grid', 'nx_grid')) dst['dswrf'].long_name = "Downward Long Wave Radiation Flux" dst['dswrf'].standard_name = "surface_downwelling_shortwave_flux_in_air" dst['dswrf'].units = "W/m^2" - dst['dswrf'][:,:,:]=ds['msdwswrf'][idx:idx+26:3,::-1,:] + dst['dswrf'][:,:,:]=ds['msdwswrf'][idx:idx+25:output_interval, ::-1, :] class ERA5(SfluxDataset): @@ -273,6 +273,7 @@ def write( rad: bool = True, bbox: Bbox = None, overwrite: bool=False, + output_interval: int = 1, ): self.start_date=start_date self.rnday=rnday @@ -284,6 +285,10 @@ def write( np.timedelta64(1, 'D'), dtype='datetime64')} + #write sflux_inputs.txt + with open(f'{outdir}/sflux_inputs.txt', 'w') as f: + f.write('&sflux_inputs\n/\n') + logger.info('Start downloading ERA5') self.inventory = ERA5DataInventory( self.start_date, @@ -300,4 +305,5 @@ def write( times=nc4.num2date(time1,units=time1.units,only_use_cftime_datetimes=False) for iday, date in enumerate(dates): - put_sflux_fields(iday, date, times, ds, nx_grid, ny_grid, air=air, rad=rad, prc=prc, OUTDIR=outdir) + put_sflux_fields(iday, date, times, ds, nx_grid, ny_grid, air=air, rad=rad, prc=prc, output_interval=output_interval, OUTDIR=outdir) + diff --git a/pyschism/forcing/source_sink/nwm.py b/pyschism/forcing/source_sink/nwm.py index b745cf83..0369d3e1 100644 --- a/pyschism/forcing/source_sink/nwm.py +++ b/pyschism/forcing/source_sink/nwm.py @@ -113,6 +113,8 @@ def __init__(self, hgrid: Gr3, nwm_file=None, workers=-1): raise IOError( "No National Water model intersections found on the mesh.") intersection = gpd.GeoDataFrame(data, crs=hgrid.crs) + #TODO: add exporting intersection as an option + #intersection.to_file('intersections.shp') del data # 2) Generate element centroid KDTree @@ -370,12 +372,14 @@ def __new__( return AWSHindcastInventory.__new__(cls) # GOOGLEHindcastInventory -> January 2019 through 30 days earlier than today + # data before April 20, 2021 (including Apirl 20) is 3-hr interval, after that is hourly + # Missing data 20211231, 20220101, 20220102, 20220103 elif start_date >= dates.localize_datetime( datetime(2021, 1, 1, 0, 0) - ) and start_date + rnday < dates.nearest_zulu() - timedelta(days=30): + ) and start_date + rnday < dates.nearest_zulu() - timedelta(days=2): return GOOGLEHindcastInventory.__new__(cls) - elif start_date >= dates.nearest_zulu() - timedelta(days=30): + elif start_date >= dates.nearest_zulu() - timedelta(days=2): return AWSForecastInventory.__new__(cls) else: @@ -561,7 +565,7 @@ def __init__( self, start_date: datetime = None, rnday: Union[int, float, timedelta] = timedelta(days=5.0), - product='analysis_assim.channel_rt', + product='medium_range_mem1', verbose=False, fallback=True, cache=None, @@ -569,7 +573,6 @@ def __init__( 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 @@ -588,37 +591,54 @@ def __init__( dtype='datetime64') for requested_time, _ in self._files.items(): - #print(f'Requesting NWM data for time {requested_time}') - print(f'Requesting NWM data for time {requested_time}') + logger.info(f'Requesting NWM data for time {requested_time}') - self._files[requested_time] = self.request_data(requested_time) + if requested_time.hour == 0: + requested_time2 = requested_time - timedelta(days=1) + else: + requested_time2 = requested_time - def request_data(self, request_time): + self._files[requested_time] = self.request_data(requested_time, requested_time2) + + def request_data(self, request_time, request_time2): + + fname = self.tmpdir / f'nwm.t00z.{self.product[:12]}.channel_rt_1.' \ + f'{request_time.strftime("%Y%m%d%H")}.conus.nc' + #logger.info(f'fname is {fname}') - fname = self.tmpdir / f'{request_time.strftime("%Y%m%d%H")}00.CHRTOUT_DOMAIN1.comp' if fname.is_file(): cached_file = list(self.tmpdir.glob(f'**/{fname.name}')) if len(cached_file) == 1: fname = cached_file[0] logger.info(f'Use cached file {fname}') - else: - - fname = f'./analysis_assim/{request_time.strftime("%Y%m%d%H")}00.CHRTOUT_DOMAIN1.comp' - - logger.info(f'Downloading file {request_time}, ') - #nwm.20180917/forcing_medium_range/nwm.t12z.medium_range.forcing.f160.conus.nc - #https://storage.googleapis.com/national-water-model/nwm.20210409/analysis_assim/nwm.t00z.analysis_assim.channel_rt.tm01.conus.nc - url=f'https://storage.googleapis.com/national-water-model/nwm.{request_time.strftime("%Y%m%d")}' \ - f'/analysis_assim/nwm.t{request_time.strftime("%H")}z.analysis_assim.channel_rt.tm00.conus.nc' - logger.info(f'{url}') + else: + + fname = f'{self.start_date.strftime("%Y%m%d")}/nwm.t00z.' \ + f'{self.product[:12]}.channel_rt_1.{request_time.strftime("%Y%m%d%H")}.conus.nc' + + logger.info(f'Downloading file {request_time}, ') + + it = request_time.strftime("%H") + if it == '00': + logger.info(f'Requesting data at 00Z from yesterday!') + it = str(int(it) + 24) + it = it.zfill(3) + + url = f'https://storage.googleapis.com/national-water-model/nwm.{request_time2.strftime("%Y%m%d")}' \ + f'/{self.product}/nwm.t00z.{self.product[:12]}.channel_rt_1.f{it}.conus.nc' + logger.info(f'{url}') + try: wget.download(url, fname) + except: + logger.info(f'No data for {request_time}!') + return fname @property def output_interval(self) -> timedelta: return { - 'analysis_assim.channel_rt': timedelta(hours=1) + 'medium_range_mem1': timedelta(hours=3) }[self.product] @property @@ -694,6 +714,7 @@ def __init__( ) for it, (requested_time, _) in enumerate(self._files.items()): + if it == 0 and requested_time.hour == 0: logger.info(f'set it=0 {requested_time}') yesterday = (self.start_date - timedelta(days=1)).strftime("%Y%m%d") diff --git a/pyschism/mesh/vgrid.py b/pyschism/mesh/vgrid.py index aa371489..d7dd7dfa 100644 --- a/pyschism/mesh/vgrid.py +++ b/pyschism/mesh/vgrid.py @@ -151,13 +151,24 @@ def open(cls, path): nvrt = int(lines[1].strip().split()[0]) - kbp = np.array([int(i.split()[1])-1 for i in lines[2:]]) - - sigma = -np.ones((len(kbp), nvrt)) - - for i, line in enumerate(lines[2:]): - sigma[i, kbp[i]:] = np.array( - line.strip().split()[2:]).astype('float') + sline = np.array(lines[2].split()).astype('float') + if sline.min() < 0: + #old version + kbp = np.array([int(i.split()[1])-1 for i in lines[2:]]) + sigma = -np.ones((len(kbp), nvrt)) + + for i, line in enumerate(lines[2:]): + sigma[i, kbp[i]:] = np.array( + line.strip().split()[2:]).astype('float') + + else: + #new version + sline = sline.astype('int') + kbp = sline-1 + sigma = np.array([line.split()[1:] for line in lines[3:]]).T.astype('float') + #replace -9. with -1. + fpm = sigma<-1 + sigma[fpm] = -1 return cls(sigma) diff --git a/setup.cfg b/setup.cfg index 46243bfd..8514b55a 100644 --- a/setup.cfg +++ b/setup.cfg @@ -1,5 +1,6 @@ [metadata] name = pyschism +version = 1.0.0 author = "SCHISM development group" author_email = jreniel@gmail.com, lcui@vims.edu description = Python package for working with SCHISM input and output files. diff --git a/setup.py b/setup.py index 3850c702..f61b4bd3 100755 --- a/setup.py +++ b/setup.py @@ -9,18 +9,19 @@ subprocess.check_call([sys.executable, '-m', 'pip', 'install', 'wheel']) -try: - from dunamai import Version -except ImportError: - subprocess.check_call( - [sys.executable, '-m', 'pip', 'install', 'dunamai'] - ) - from dunamai import Version # type: ignore[import] +#try: +# from dunamai import Version +#except ImportError: +# subprocess.check_call( +# [sys.executable, '-m', 'pip', 'install', 'dunamai'] +# ) +# from dunamai import Version # type: ignore[import] -try: - version = Version.from_any_vcs().serialize() -except RuntimeError: - version = '0.0.0' +#try: +# print(Version.from_any_vcs().serialize()) +# version = Version.from_any_vcs().serialize() +#except RuntimeError: +# version = '0.0.0' class BuildSchism(setuptools.Command): @@ -81,7 +82,7 @@ def run(self): meta = conf['metadata'] setuptools.setup( name=meta['name'], - version=version, + version=meta['version'], author=meta['author'], author_email=meta['author_email'], description=meta['description'],