diff --git a/dev/compute_interpolation_weights/README.txt b/dev/compute_interpolation_weights/README.txt new file mode 100644 index 0000000..4b5bf7f --- /dev/null +++ b/dev/compute_interpolation_weights/README.txt @@ -0,0 +1,23 @@ +Scripts to create interpolation weights and distance to boundary files for interpolating forcing to RWPS mesh. +To run: + +$ cd RWPS/dev/compute_interpolation_weights/sorc +$ sh compute_interpolation_weights.sh oc_1500m_30km + +to generate interpolation files for mesh rwps.oc_1500m_30km.msh. Files will be writen to directory: +RWPS/fix/ + +Interpolation weights are created for: + +nbm oc domain +rrfs hi domain (used for wind) +rrfs pr domain (used for wind) +rrfs ak domain (used for wind) +rrfs na domain (used for wind) +rrfs conus domain (used for wind) +nbm ak domain (used for ice concentration) +rtofs glo domain (used for current without extrapolation) +rtofs glo domain (used for ice with extrapolation) +stofs domain (used for waterlevel and current) + +Currently only setup to work on wcoss2. diff --git a/dev/compute_interpolation_weights/ecf/compute_interpolation_weights.ecf b/dev/compute_interpolation_weights/ecf/compute_interpolation_weights.ecf new file mode 100644 index 0000000..af2bc7b --- /dev/null +++ b/dev/compute_interpolation_weights/ecf/compute_interpolation_weights.ecf @@ -0,0 +1,26 @@ +#!/bin/bash +#PBS -N ESMPy +#PBS -j oe +#PBS -S /bin/bash +#PBS -q dev +#PBS -A NWPS-DEV +#PBS -l walltime=04:00:00 +#PBS -l select=1:ncpus=16:mem=500GB +#PBS -l place=excl +#PBS -l debug=true + +module reset +module load PrgEnv-intel/8.5.0 +module load intel/19.1.3.304 +module load craype/2.7.17 +module load cray-mpich/8.1.19 +module load hdf5-C/1.14.0 +module load netcdf-C/4.9.2 +module load esmf-C/8.6.0 +module load ve/hafs/2.1 + +pip list -v + +$RWPSroot/dev/compute_interpolation_weights/jobs/COMPUTE_INTERPOLATION_WEIGHTS + + diff --git a/dev/compute_interpolation_weights/jobs/COMPUTE_INTERPOLATION_WEIGHTS b/dev/compute_interpolation_weights/jobs/COMPUTE_INTERPOLATION_WEIGHTS new file mode 100755 index 0000000..a44dd71 --- /dev/null +++ b/dev/compute_interpolation_weights/jobs/COMPUTE_INTERPOLATION_WEIGHTS @@ -0,0 +1,171 @@ +#!/bin/bash + +###################################################################### +# This job computes interpolation weights and distance to boundary for +# a RWPS mesh. Currently configured for: +# nbm oc domain (wind) +# nbm ak domain (ice) +# rrfs hi domain (wind) +# rrfs pr domain (wind) +# rrfs ak domain (wind) +# rrfs conus domain (wind) +# rrfs na domain (wind) +# rtofs glo domain (surface currents, ice) +# stofs v3 (waterlevel, currents) +# +# Different sets of interpolation weights are created for rtofs ice +# and surface currents so that ice concentration can be extrapolated +# by nearest node while currents are not extrapolated into coastal +# regions outside of rtofs marine domain. +# +# Files are created in local RWPS/fix directory. +###################################################################### + +cd $RWPSroot/dev/compute_interpolation_weights/ush + +# Sample input directory +SID="/lfs/h2/emc/couple/noscrub/keston.smith/SampleInput" +date="20260730" +cycl="00" +Nprocs=32 + +out=$RWPSroot/fix + +meshname="${mesh##*/}" +meshname="${meshname: 0: -4}" + + +###################################################################### +# Compute rrfs interpolation weights for interpolating wind. +# no extrapolation for rrfs forecasts. +###################################################################### +rrfs_hi=$SID/wind.$date.$cycl/rrfs.$date.$cycl.wind10m.hi.nc +rrfs_pr=$SID/wind.$date.$cycl/rrfs.$date.$cycl.wind10m.pr.nc +rrfs_ak=$SID/wind.$date.$cycl/rrfs.$date.$cycl.wind10m.ak.nc +rrfs_na=$SID/wind.$date.$cycl/rrfs.$date.$cycl.wind10m.na.nc +rrfs_conus=$SID/wind.$date.$cycl/rrfs.$date.$cycl.wind10m.conus.nc + + +rrfsIntrp="$out/InterpolationWeights.$meshname.rrfs" +rrfsDst="$out/DistToBndy.$meshname.rrfs" + +rm -rf $rrfsIntrp.hi.nc $rrfsDst.hi.nc __pycache__ +python ./compute_gridded_to_rwps_interp_weights.py $rrfs_hi $mesh $rrfsIntrp.hi.nc $rrfsDst.hi.nc -1 > GenWhts.rrfs.hi.out +rm -rf $rrfsIntrp.pr.nc $rrfsDst.pr.nc __pycache__ +python ./compute_gridded_to_rwps_interp_weights.py $rrfs_pr $mesh $rrfsIntrp.pr.nc $rrfsDst.pr.nc -1 > GenWhts.rrfs.pr.out +rm -rf $rrfsIntrp.ak.nc $rrfsDst.ak.nc __pycache__ +python ./compute_gridded_to_rwps_interp_weights.py $rrfs_ak $mesh $rrfsIntrp.ak.nc $rrfsDst.ak.nc -1 > GenWhts.rrfs.ak.out +rm -rf $rrfsIntrp.na.nc $rrfsDst.na.nc __pycache__ +python ./compute_gridded_to_rwps_interp_weights.py $rrfs_na $mesh $rrfsIntrp.na.nc $rrfsDst.na.nc -1 > GenWhts.rrfs.na.out +rm -rf $rrfsIntrp.conus.nc $rrfsDst.conus.nc __pycache__ +python ./compute_gridded_to_rwps_interp_weights.py $rrfs_conus $mesh $rrfsIntrp.conus.nc $rrfsDst.conus.nc -1 > GenWhts.rrfs.conus.out + + +###################################################################### +# Compute nbm interpolation weights on oc domain for interpolating wind. +# Extrapolate to cover minor edge differences on south and east boundaries. +###################################################################### +nbm_oc=$SID/wind.$date.$cycl/nbm.$date.$cycl.wind10m.oc.nc +nbmIntrp="$out/InterpolationWeights.$meshname.nbm" +nbmDst="$out/DistToBndy.$meshname.nbm" + +#extrapolate nbm oc wind to all points +rm $nbmIntrp.oc.nc $nbmDst.oc.nc +python ./compute_gridded_to_rwps_interp_weights.py $nbm_oc $mesh $nbmIntrp.oc.nc $nbmDst.oc.nc 1 > GenWhts.nbm.oc.out + + +###################################################################### +# Compute nbm interpolation weights on ak domain for interpolating ice. +# No extrapolation outside the ak region +###################################################################### +nbm_ak=$SID/ice.$date.$cycl/nbm.$date.$cycl.ice.ak.nc +#do not extrapolate for NBM .ak. domain +rm $nbmIntrp.ak.nc $nbmDst.ak.nc +python ./compute_gridded_to_rwps_interp_weights.py $nbm_ak $mesh $nbmIntrp.ak.nc $nbmDst.ak.nc -1 > GenWhts.nbm.ak.out + +rtofs_ice=$SID/ice.$date.$cycl/rtofs.ice.$date.nc +rtofsIntrp="$out/InterpolationWeights.$meshname.rtofs" +rtofsDst="$out/DistToBndy.$meshname.rtofs" + +###################################################################### +# Compute rtofs interpolation weights for interpolating ice. +# Do extrapolate rtofs .glo. for full ice coverage in near shore. +###################################################################### +rm $rtofsIntrp.ice.nc $rtofsDst.ice.nc +python ./compute_gridded_to_rwps_interp_weights.py $rtofs_ice $mesh $rtofsIntrp.ice.nc $rtofsDst.ice.nc 1 > GenWhts.rtofs.ice.out + + +###################################################################### +# Compute rtofs interpolation weights for interpolating currents. +# Do not extrapolate into coastal regions +###################################################################### +rtofs_cur=$SID/rtofs.$date/currents.$date.rtofs.nc +#do not extrapolate rtofs .glo. currents +rm $rtofsIntrp.current.nc $rtofsDst.current.nc +python ./compute_gridded_to_rwps_interp_weights.py $rtofs_cur $mesh $rtofsIntrp.current.nc $rtofsDst.current.nc -1 > GenWhts.rtofs.current.out + + +###################################################################### +# Compute stofs interpolation weights for stofs v3 output +# All dates after ~ August 18 2026 +###################################################################### +stofs_glo="$SID/stofs.v3/stofs.20260819.00/stofs_2d_glo.t00z.fields.cwl.nc" + +mkdir -p STOFSInterpWeights.$meshname + +# First construct parallel jobscript to compute weights using Nprocs division of destination domain(meshfile) +python compute_unstr_to_rwps_interp_weights.py $stofs_glo $mesh $Nprocs + +# Now run the parallel job script +qsub -W block=true jobcardcompute_unstr_to_rwps_interp_weightsPBS +wait + +# combine text file output of parallel jobs +cat STOFSInterpWeights.$meshname/Part.IntrpWghts.*.txt > InterpWeights.$meshname.stofs.txt + +# sort order of nodes and convert output weights to netcdf file +python convert_weights_to_netcdf.py $stofs_glo $mesh InterpWeights.$meshname.stofs.txt 1 +cp InterpWeights.$meshname.stofs.nc $out/InterpolationWeights.$meshname.stofs.nc +cp DistToBndy.$meshname.stofs.nc $out/DistToBndy.$meshname.stofs.nc + +# clean up temporary directory and files +rm -rf STOFSInterpWeights.$meshname +rm InterpWeights.$meshname.stofs.txt +rm *.out +rm ESMPy* + + +###################################################################### +# Compute stofs interpolation weights for old(v2) stofs output +# All dates before ~ August 18 2026 +# Note: stofs v3 and stofs v2 appear to use the same mesh but file sizes differ. +# There should not be any problem just useing the stofs v3 based weights for v2 output files +###################################################################### +stofs_glo="$SID/stofs.20260730.00/stofs_2d_glo.t00z.fields.cwl.nc" + +mkdir -p STOFSInterpWeights.$meshname + +# First construct parallel jobscript to compute weights using Nprocs division of destination domain(meshfile) +python compute_unstr_to_rwps_interp_weights.py $stofs_glo $mesh $Nprocs + +# Now run parallel job script +qsub -W block=true jobcardcompute_unstr_to_rwps_interp_weightsPBS +wait + +#combine text file output of parallel jobs +cat STOFSInterpWeights.$meshname/Part.IntrpWghts.*.txt > InterpWeights.$meshname.stofs.txt + +# convert output weights to netcdf file +python convert_weights_to_netcdf.py $stofs_glo $mesh InterpWeights.$meshname.stofs.txt 1 +cp InterpWeights.$meshname.stofs.nc $out/InterpolationWeights.$meshname.stofs.v2.nc +cp DistToBndy.$meshname.stofs.nc $out/DistToBndy.$meshname.stofs.v2.nc +mv DistToBndy.$meshname.stofs.nc DistToBndy.$meshname.stofs.v2.nc +mv InterpWeights.$meshname.stofs.nc InterpWeights.$meshname.stofs.v2.nc + +# clean up temporary directory and files +rm -rf STOFSInterpWeights.$meshname +rm InterpWeights.$meshname.stofs.txt +rm *.out +rm ESMPy* + + diff --git a/dev/compute_interpolation_weights/ush/compute_gridded_to_rwps_interp_weights.py b/dev/compute_interpolation_weights/ush/compute_gridded_to_rwps_interp_weights.py new file mode 100644 index 0000000..4852594 --- /dev/null +++ b/dev/compute_interpolation_weights/ush/compute_gridded_to_rwps_interp_weights.py @@ -0,0 +1,219 @@ +# This is a a set of routines that interpolate wind forecasts to an unstructured +# mesh from regular or curvilinear grids (using ESMPY). In addition to interpolation +# field the out file contains a spatially variable error variance estimate based on the +# distance to the forecast boundary. This error variance is used later to update +# dispirate forecasts in a bayesian manner to yield a spatially smooth estimate +# incorporating all forecasts. The understanding is that the more local forecast +# products are of higher accuracy than the coarser broader scale forecasts + +import numpy as np +import os + +import datetime +import netCDF4 as nc +import sys +import re + +import xarray as xr +import esmpy +import scipy.sparse as sp + +sys.path.append(os.path.abspath(os.path.join(os.path.dirname(__file__), '../../../ush/preprocess'))) + +import interp_utilities as iutil + +# Main program +AddExtrapolationSupport=True + +nargin = len(sys.argv) - 1 + +flin=sys.argv[1] +mshfl=sys.argv[2] +weights_file=sys.argv[3] +dist2bnd_file=sys.argv[4] +# Don't use nearest neighbor interpolation unless 6th positive integer argument present. +# This may be needed on boundary of RWPS mesh if node alignment is outside NBM OC domain +Extrapolate=False +if nargin > 4: + if int(sys.argv[5])>0: + print("using nearest neighbor to extrapolate wind field beyond geometric coverage") + Extrapolate=True + +xi, yi, ei, zi = iutil.loadWW3Mesh(mshfl) +nn=len(xi) +#shift coords to[129,370] +data = nc.Dataset(flin,"r") + +#read spaital dimensions and determine if input mesh is curvilinear or regular +x1=np.asarray(data["longitude"][:]) +y1=np.asarray(data["latitude"][:]) + +# rwps :[-231, 11] +# nbm :[129, 370] +# rrfs,pr :[ -75.5000 -62.5087] +# rrfs,hi :[ -161.5250 -153.8690] +# rrfs,ak :[150.2012 266.2886] +# rrfs,na :[67.5721 427.0000] +# rrfs,conus :[ 225.9045 299.0828] +# rtofs,glo :[74.1552,434.0146] +# stofs, glo :[-180,180] + +#Shift longitude coordinates to RWPS specs for various files +if "hi" in flin: + x1=x1+360. +if "pr" in flin: + x1=x1+360. + +if "rtofs" in flin: #remove bad geometry edges + x1=x1[1:-1,1:-1] + y1=y1[1:-1,1:-1] + x1=x1-360. + +if (len(x1.shape)==2 and len(y1.shape)==2): + IsCrvLn=True +elif (len(x1.shape)==1 and len(y1.shape)==1): + #represent regular grid as curvilinear grid + IsCrvLn=False + nx=len(x1) + ny=len(y1) + x1 = np.tile(x1,(ny,1)) + y1 = np.tile(y1,(nx,1)).T + IsCrvLn=True +else: + print("input file spatial dimension is not recognized. ending program") + sys.exit() + +####################################### +# === Create weights ===# +####################################### +meshslash=mshfl.rfind('/')+1 +dom=flin.split(".") +dom=dom[len(dom)-2] + +print("interpolation weights will be written to file = "+ weights_file) + +nx=x1.shape[0] +ny=x1.shape[1] +n1=nx*ny + +print("Computing weights and saving to file: "+ weights_file) + +#Use esmpy to construct bilinear interpolation weights +iutil.CurvilinearGridCreateInterpWeights(xi, yi, x1, y1, weights_file) + +####################################### +# === read weights and ===# +####################################### +with xr.open_dataset(weights_file) as ds_s: + # Standard sparse storage uses 'row', 'col', and 'S'==weights variables + row = ds_s['row'].values + col = ds_s['col'].values + weights = ds_s['S'].values + Nrows=ds_s.attrs.get('Nrows') + Ncols=ds_s.attrs.get('Ncols') +print("nn = "+str(nn)+": Nrows = "+str(Nrows)) +print("n1 = "+str(n1)+": Ncols = "+str(Ncols)) +if not ((nn==Nrows) and (n1==Ncols)): + print("Wrong matrix weights: number of rows from "+ mshfl +" = "+str(nn)+ + " but number of rows in "+ weights_file +" = "+str(Nrows)+ + ", number of spatial points in "+ flin +" = "+str(n1)+ + " but number of columns in "+ weights_file +" = "+str(Ncols) ) + print(" You probably need to remove file "+ weights_file +" and rerun to generate appropriate weights") +matrix = sp.coo_matrix((weights, (row-1, col-1)), shape=(nn,n1)).tocsr() +print("sparse interpolation matrix") +print(matrix) + +################################################################################## +# START: Extrapolate for nodes not covered by interpolator +################################################################################## +x1v=np.transpose(x1).reshape(n1) # vectorize src nodes, consistant with data to interpolate +y1v=np.transpose(y1).reshape(n1) + +if Extrapolate: + from scipy.interpolate import NearestNDInterpolator + srcp = np.array((x1v,y1v)).T + srcv = 1.+x1v**2 + y1v**2 #dummy input field + #dstv = matrix @ srcv.T + row_sum = matrix.sum(axis=1) + j0=np.where( row_sum==0 ) # destination nodes with no coverage from interpolation matrix + j0=np.array(j0[0]).tolist() + dstp = np.array((xi[j0],yi[j0])).T + interpolator = NearestNDInterpolator(srcp, srcv) + distances, j0src = interpolator.tree.query(dstp) + weightsExtrp=weights.tolist().append([1.0] * len(j0) ) + rowExtrp=np.concatenate( (row, np.array(j0)) ) + colExtrp=np.concatenate( (col, np.array(j0src+1)) ) + weightsExtrp=np.concatenate( (weights, np.array([1.0] * len(j0))) ) + os.replace(weights_file, weights_file[0:-3]+".NoExtrap.nc") + iutil.WriteInterpolationWeightsToNetCDF(weights_file,rowExtrp,colExtrp,weightsExtrp,len(xi),len(x1v)) +################################################################################## +# FINISHED: Extrapolate for nodes not covered by interpolator +################################################################################## + +################################################################################## +# START: Extrapolation support for NaN occurances in source field +################################################################################## +if AddExtrapolationSupport: + with nc.Dataset(weights_file, 'r+', format='NETCDF4') as ncadd: + ncadd.createDimension('nn_src' , len(x1v)) + ncadd.createDimension('nn_dst' , len(xi)) + + xsrc_var=ncadd.createVariable('x_src', 'f8', ('nn_src',)) + xsrc_var.long_name = 'interpolation source node longitude' + xsrc_var[:]=x1v[:] + + ysrc_var=ncadd.createVariable('y_src', 'f8', ('nn_src',)) + ysrc_var.long_name = 'interpolation source node latitude' + ysrc_var[:]=y1v[:] + + xdst_var=ncadd.createVariable('x_dst', 'f8', ('nn_dst',)) + xdst_var.long_name = 'interpolation destination node longitude' + xdst_var[:]=xi[:] + + ydst_var=ncadd.createVariable('y_dst', 'f8', ('nn_dst',)) + ydst_var.long_name = 'interpolation destination node latitude' + ydst_var[:]=yi[:] + +################################################################################## +# FINISHED: Extrapolation support for NaN occurances in source field +################################################################################## +with nc.Dataset(weights_file, 'r+', format='NETCDF4') as ncadd: + ncadd.setncattr("SrcFieldType", "gridded") + ncadd.setncattr("InputFile", flin) + ncadd.setncattr("MeshFile", mshfl) + + +################################################################################## +# START: Compute distance to boundary for each node in mesh: +################################################################################## +row_sum = matrix.sum(axis=1) +j0=np.where( row_sum==0 ) # destination nodes with no coverage from interpolation matrix +j0=np.array(j0[0]).tolist() +u0=np.ones(xi.shape) +nan=float("nan") +u0[j0]=nan + +if not j0: #all points are interpolation envelope + dist2bnd=0.*xi + np.inf +else: + dist2bnd=iutil.CalculateDistanceToInterpEnvelope(xi,yi,u0, 1.) + +if Extrapolate: + dist2bnd=0.*dist2bnd + np.inf #all points are inside boundary- No boundary with this type of extrapolation + print(np.mean(dist2bnd)) + + +with nc.Dataset(dist2bnd_file, 'w', format='NETCDF4') as ncout: + ncout.createDimension('node' , nn) + d_var=ncout.createVariable('dist2bnd', 'f4', ('node',)) + d_var.long_name = 'distance to boundary' + d_var.units = 'km' + d_var.standard_name = 'distance to boundary' + d_var[:]=dist2bnd[:] + + z_var=ncout.createVariable('depth', 'f4', ('node',)) + z_var.long_name = 'mesh depth' + z_var.units = 'm' + z_var.standard_name = 'depth' + z_var[:]=zi[:] + diff --git a/dev/compute_interpolation_weights/ush/compute_interpolation_weights.sh b/dev/compute_interpolation_weights/ush/compute_interpolation_weights.sh new file mode 100755 index 0000000..6c2afed --- /dev/null +++ b/dev/compute_interpolation_weights/ush/compute_interpolation_weights.sh @@ -0,0 +1,53 @@ +!/bin/bash + +# --------------------------------------------------------------------------- # +# Launch scripts to compute interpolation weights for an RWPS mesh. To run: # +# # +# $sh ./compute_interpolation_weights.sh oc_500m_10km # +# # +# to generate interpolation weights for unstructured mesh # +# rwps.oc_500m_10km.msh. # +# # +# Last Changed : 09-04-2026 # +# --------------------------------------------------------------------------- # + +echo 'setting paths...' + +export meshID=$1 + +export RWPSroot=$(cd "$(dirname "$(readlink -f -n "${BASH_SOURCE[0]}")")" && git rev-parse --show-toplevel) + +export mesh="$RWPSroot/fix/rwps.$meshID.msh" + + +export fix="$RWPSroot/fix" +export prep="$RWPSroot/PrepInputs" +export tmp="$prep/tmpfiles" +export frc="$prep/forcing" +export outdir=$prep + +mkdir -p $prep +mkdir -p $tmp +mkdir -p $frc + + +#machine dependend path to rtofs, nbm, rrfs, and stofs forecast files +export COMINrtofs="/lfs/h1/ops/prod/com/rtofs/v2.5/rtofs.$PDY/" +export COMINnbm="/lfs/h3/mdl/ptmp/mdl.nbm/blend/v5.2/blend.$PDY/$cyc/grib2" +export COMINrrfs="/lfs/h1/ops/prod/com/rrfs/v1.0/rrfs.$PDY/$cyc" +export COMINstofs="/lfs/h1/ops/prod/com/stofs/v3.1/stofs_2d_glo.$PDY" + +export COMINlocal=$tmp + +#machine dependend path to RWPS fix files +export RWPSfix=/lfs/h2/emc/couple/noscrub/keston.smith/RWPS + +# copy mesh to local fix directory +cp -p $RWPSfix/fix/$meshID/20260722/rwps.$meshID.msh $RWPSroot/fix/ + +meshname="${mesh##*/}" +export meshname="${meshname: 0: -4}" + +#Retrieve current and process for forecast cycle +qsub -V $RWPSroot/dev/compute_interpolation_weights/ecf/compute_interpolation_weights.ecf + diff --git a/dev/compute_interpolation_weights/ush/compute_unstr_to_rwps_interp_weights.py b/dev/compute_interpolation_weights/ush/compute_unstr_to_rwps_interp_weights.py new file mode 100644 index 0000000..dd000e8 --- /dev/null +++ b/dev/compute_interpolation_weights/ush/compute_unstr_to_rwps_interp_weights.py @@ -0,0 +1,138 @@ +import numpy as np +import netCDF4 as nc +import sys +from scipy.interpolate import RegularGridInterpolator +import os + +sys.path.append(os.path.abspath(os.path.join(os.path.dirname(__file__), '../../../ush/preprocess'))) +import interp_utilities as iutil + +nargin = len(sys.argv) - 1 + +flin=sys.argv[1] +mshfl=sys.argv[2] + +meshslash=mshfl.rfind('/')+1 +TmpOutDir="STOFSInterpWeights."+mshfl[meshslash:len(mshfl)-4] + +#if nargin==3 then write the parallel jobcard that carries out the interpolation weight calculation +if nargin ==3 : + Njobs=int(sys.argv[3]) + os.makedirs(TmpOutDir, exist_ok=True) + try: + os.remove(TmpOutDir+"/*.txt") + except: + print("directory "+TmpOutDir+" is alread empty") + iutil.WriteInterpJobscriptPBS("jobcardcompute_unstr_to_rwps_interp_weightsPBS",flin,mshfl,Njobs, 1) + iutil.WriteInterpJobscriptSLURM("jobcardcompute_unstr_to_rwps_interp_weightsSLURM",flin,mshfl,Njobs, 1) + print("Made parallel jobcards to create interpolation weights with "+str(Njobs)+" processes. Next step:") + print("sbatch jobcardcompute_unstr_to_rwps_interp_weightsSLURM") + sys.exit() + +#if nargin==4 then do the part of the mesh for jobID of Njobs +if nargin == 4: + jobID=int(sys.argv[3]) + Njobs=int(sys.argv[4]) + +#if nargin==6 then do the interpolation weight calculation for a specified window +else: + lonW=int(sys.argv[3]) + lonE=int(sys.argv[4]) + latS=int(sys.argv[5]) + latN=int(sys.argv[6]) + +meshslash=mshfl.rfind('/')+1 +if nargin < 5: + weights_file = TmpOutDir+"/Part.IntrpWghts."+str(jobID)+".txt" +else: + weights_file = TmpOutDir+"/Part.IntrpWghts.W"+str(lonW)+".E"+str(lonE)+".S"+str(latS)+".N"+str(latN)+"."+mshfl[meshslash:len(mshfl)-3]+"txt" + + +print("saving output to:"+weights_file) + +xi, yi, ei, zi = iutil.loadWW3Mesh(mshfl) +nni=len(xi) +j=np.where(xi>90.) # should be empty +xi[j]=xi[j]-360. +#if np.mean(xi)<0.: +# xi=xi+360. + +# divide domain from east west based on job id=0 ... Njobs +# and make North-South window contain the full target domain +# NOTE: for efficiency the sections of the domain should have approximately +# the same number of destination nodes and/or the same number of source elements + +if nargin < 5: +# balance node load for + xis=np.sort(xi) + mm=round(nni/Njobs) + xil=xis[range(0,nni,mm)] + lonW=xil[jobID] + if jobID==Njobs-1: + lonE=xis[-1]+1. + else: + lonE=xil[jobID+1] + if jobID==0: + lonW=xis[0]-1. + + latS=np.min(yi)-1. + latN=np.max(yi)+1. + +data = nc.Dataset(flin,"r") +#read spaital dimensions and determine if input mesh is curvilinear or regular +x=np.asarray(data["x"][:]) +j=np.where(x>90.) #elements broken here should not effect RWPS interpolation +x[j]=x[j]-360. + +y=np.asarray(data["y"][:]) +e=np.asarray(data["element"][:,:]) + +e0=e-1 +xc = np.mean(x[e0], axis=1) +yc = np.mean(y[e0], axis=1) + +nn=len(x) + +SearchWidth=.5 # search rectangle extended beyond window to this number (degrees lat lon) + +print(e.shape) +ne=e.shape[0] +print(ne) + +xwin=[lonW,lonE] +ywin=[latS,latN] +print("calculating interpolation weights for [W"+str(lonW)+": E "+str(lonE)+": S "+str(latS)+": N "+str(latN)+"]") +#Find target mesh nodes in window +jxUi=np.where( xi <= np.max(xwin))[0].tolist() +jxDi=np.where( xi >= np.min(xwin))[0].tolist() +jyUi=np.where( yi <= np.max(ywin))[0].tolist() +jyDi=np.where( yi >= np.min(ywin))[0].tolist() +ji=list( set(jxUi) & set(jxDi) & set(jyUi) & set(jyDi) ) + +#Find source elements in and near window +SearchWidth=.5 +jxU=np.where( xc < np.max(xwin)+SearchWidth )[0].tolist() +jxD=np.where( xc > np.min(xwin)-SearchWidth )[0].tolist() +jyU=np.where( yc < np.max(ywin)+SearchWidth )[0].tolist() +jyD=np.where( yc > np.min(ywin)-SearchWidth )[0].tolist() +je=list( set(jxU) & set(jxD) & set(jyU) & set(jyD) ) + +weights, nodes, elenum, Dist2EleCenter=iutil.compute_mesh_to_mesh_interp_weights(x, y, e[je,:], xi[ji], yi[ji]) + +#Move back to global index's +for k in range(len(elenum)): + elenum[k]=je[elenum[k]] + ji[k]=ji[k] + +#return to convention numbering from 1 ... +elenum = list(np.array(elenum) + 1) +ji = list(np.array(ji) + 1) +#nodes = np.array(nodes) + 1 # as writen using the correct 1 .. nn node numbering + +# File format +# ["taget node number" "source element number" "source n1" "source n2 "source n3" "distance from target node to element center" "weight1" "weight2" "weight3"] + +Fout=np.vstack((np.array(ji),elenum,nodes[:,0],nodes[:,1],nodes[:,2],Dist2EleCenter,weights[:,0],weights[:,1],weights[:,2])) +np.savetxt(weights_file,Fout.T, fmt='%d %d %d %d %d %.6f %.6f %.6f %.6f') + + diff --git a/ecf/jrwps_prep_current.ecf b/ecf/jrwps_prep_current.ecf new file mode 100644 index 0000000..c176f6e --- /dev/null +++ b/ecf/jrwps_prep_current.ecf @@ -0,0 +1,36 @@ +#!/bin/sh +#PBS -q dev +#PBS -l walltime=00:59:00 +#PBS -A GLWU-DEV +#PBS -N RWPS.PreProc.current +#PBS -j oe +#PBS -l select=1:ncpus=32:mem=64GB +#PBS -o RWPS.PreProc.out + +cd $HOMErwps/ush/preprocess + +module use ${HOMErwps}/modulefiles + +############################################################ +# Load modules +############################################################ +source $HOMErwps/versions/preproc.${MACHINE_ID}.ver + +#modules for esmf +module reset +module load PrgEnv-intel/${PrgEnv_intel_ver} +module load craype/${craype_ver} +module load intel/${intel_ver} +module load cray-mpich/${cray_mpich_ver} +module load hdf5-C/${hdf5_C_ver} +module load netcdf-C/${netcdf_C_ver} +module load esmf-C/${esmf_C_ver} +module load ve/hafs/${ve_hafs_ver} + +#modules for wgrib2 +module load intel-oneapi/${intel_oneapi_ver} +module load wgrib2/${wgrib2_ver} + +module list + +${HOMErwps}/jobs/JRWPS_PREP_CURRENT diff --git a/ecf/jrwps_prep_ice.ecf b/ecf/jrwps_prep_ice.ecf new file mode 100644 index 0000000..8a771be --- /dev/null +++ b/ecf/jrwps_prep_ice.ecf @@ -0,0 +1,36 @@ +#!/bin/sh +#PBS -q dev +#PBS -l walltime=00:59:00 +#PBS -A GLWU-DEV +#PBS -N RWPS.PreProc.ice +#PBS -j oe +#PBS -l select=1:ncpus=16:mem=64GB +#PBS -o RWPS.PreProc.out + +cd $HOMErwps/ush/preprocess + +module use ${HOMErwps}/modulefiles + +############################################################ +# Load modules +############################################################ +source $HOMErwps/versions/preproc.${MACHINE_ID}.ver + +#modules for esmf +module reset +module load PrgEnv-intel/${PrgEnv_intel_ver} +module load craype/${craype_ver} +module load intel/${intel_ver} +module load cray-mpich/${cray_mpich_ver} +module load hdf5-C/${hdf5_C_ver} +module load netcdf-C/${netcdf_C_ver} +module load esmf-C/${esmf_C_ver} +module load ve/hafs/${ve_hafs_ver} + +#modules for wgrib2 +module load intel-oneapi/${intel_oneapi_ver} +module load wgrib2/${wgrib2_ver} + +module list + +${HOMErwps}/jobs/JRWPS_PREP_ICE diff --git a/ecf/jrwps_prep_waterlevel.ecf b/ecf/jrwps_prep_waterlevel.ecf new file mode 100644 index 0000000..b4ee30d --- /dev/null +++ b/ecf/jrwps_prep_waterlevel.ecf @@ -0,0 +1,36 @@ +#!/bin/sh +#PBS -q dev +#PBS -l walltime=00:59:00 +#PBS -A GLWU-DEV +#PBS -N RWPS.PreProc.waterlevel +#PBS -j oe +#PBS -l select=1:ncpus=16:mem=64GB +#PBS -o RWPS.PreProc.out + +cd $HOMErwps/ush/preprocess + +module use ${HOMErwps}/modulefiles + +############################################################ +# Load modules +############################################################ +source $HOMErwps/versions/preproc.${MACHINE_ID}.ver + +#modules for esmf +module reset +module load PrgEnv-intel/${PrgEnv_intel_ver} +module load craype/${craype_ver} +module load intel/${intel_ver} +module load cray-mpich/${cray_mpich_ver} +module load hdf5-C/${hdf5_C_ver} +module load netcdf-C/${netcdf_C_ver} +module load esmf-C/${esmf_C_ver} +module load ve/hafs/${ve_hafs_ver} + +#modules for wgrib2 +module load intel-oneapi/${intel_oneapi_ver} +module load wgrib2/${wgrib2_ver} + +module list + +${HOMErwps}/jobs/JRWPS_PREP_WATERLEVEL diff --git a/ecf/jrwps_prep_wind.ecf b/ecf/jrwps_prep_wind.ecf new file mode 100644 index 0000000..02c7e4a --- /dev/null +++ b/ecf/jrwps_prep_wind.ecf @@ -0,0 +1,36 @@ +#!/bin/sh +#PBS -q dev +#PBS -l walltime=04:00:00 +#PBS -A GLWU-DEV +#PBS -N RWPS.PreProc.wind +#PBS -j oe +#PBS -l select=1:ncpus=32:mem=64GB +#PBS -o RWPS.PreProc.out + +cd $HOMErwps/ush/preprocess + +module use ${HOMErwps}/modulefiles + +############################################################ +# Load modules +############################################################ +source $HOMErwps/versions/preproc.${MACHINE_ID}.ver + +#modules for esmf +module reset +module load PrgEnv-intel/${PrgEnv_intel_ver} +module load craype/${craype_ver} +module load intel/${intel_ver} +module load cray-mpich/${cray_mpich_ver} +module load hdf5-C/${hdf5_C_ver} +module load netcdf-C/${netcdf_C_ver} +module load esmf-C/${esmf_C_ver} +module load ve/hafs/${ve_hafs_ver} + +#modules for wgrib2 +module load intel-oneapi/${intel_oneapi_ver} +module load wgrib2/${wgrib2_ver} + +module list + +${HOMErwps}/jobs/JRWPS_PREP_WIND diff --git a/jobs/JRWPS_PREP_CURRENT b/jobs/JRWPS_PREP_CURRENT new file mode 100755 index 0000000..1457d1a --- /dev/null +++ b/jobs/JRWPS_PREP_CURRENT @@ -0,0 +1,6 @@ +#!/bin/sh + +cd $HOMErwps/ush/preprocess +sh get_current.sh +sh process_current.sh + diff --git a/jobs/JRWPS_PREP_ICE b/jobs/JRWPS_PREP_ICE new file mode 100755 index 0000000..31b77cf --- /dev/null +++ b/jobs/JRWPS_PREP_ICE @@ -0,0 +1,7 @@ +#!/bin/sh + +cd $HOMErwps/ush/preprocess +sh get_ice.sh +sh process_ice.sh + + diff --git a/jobs/JRWPS_PREP_WATERLEVEL b/jobs/JRWPS_PREP_WATERLEVEL new file mode 100755 index 0000000..0e1f6ae --- /dev/null +++ b/jobs/JRWPS_PREP_WATERLEVEL @@ -0,0 +1,5 @@ +#!/bin/sh + +cd $HOMErwps/ush/preprocess +sh get_waterlevel.sh +sh process_waterlevel.sh diff --git a/jobs/JRWPS_PREP_WIND b/jobs/JRWPS_PREP_WIND new file mode 100755 index 0000000..dcdcb4c --- /dev/null +++ b/jobs/JRWPS_PREP_WIND @@ -0,0 +1,7 @@ +#!/bin/sh + +cd $HOMErwps/ush/preprocess +sh get_wind.sh +sh process_wind.sh + + diff --git a/sorc/get_externals.sh b/sorc/get_externals.sh index 47b172e..b48d8f6 100755 --- a/sorc/get_externals.sh +++ b/sorc/get_externals.sh @@ -7,7 +7,18 @@ # --------------------------------------------------------------------------- # echo 'Fetching externals...' -cp -p /lfs/h2/emc/couple/noscrub/keston.smith/meshes/RWPS.v0.msh ../fix/ + + +# copy mesh to local fix directory +cp -p $RWPSfix/fix/$meshname/20260722/rwps.$meshname.msh ../fix/ +# copy Interpoplation weights for nbm, rrfs, rtofs and stofs to local fix directory +cp -p $RWPSfix/fix/$meshname/20260722/InterpolationWeights*$meshname*.nc ../fix/ +# copy distance to boundary for nbm, rrfs, rtofs and stofs to local fix directory +cp -p $RWPSfix/fix/$meshname/20260722/DistToBndy*$meshname*.nc ../fix/ + + +#cp -p /lfs/h2/emc/couple/noscrub/keston.smith/meshes/RWPS.v0.msh ../fix/ + cp -p /lfs/h2/emc/couple/noscrub/saeideh.banihashemi/git/fv3gfs/fix/fix_glwu/mesh.glwu ../fix/mesh.rwps cp -p /lfs/h2/emc/couple/noscrub/saeideh.banihashemi/git/fv3gfs/fix/fix_glwu/grint_weights.grlc_2p5km ../fix/ cp -p /lfs/h2/emc/couple/noscrub/saeideh.banihashemi/git/fv3gfs/fix/fix_glwu/grint_weights.grlr ../fix/ diff --git a/ush/prepare_forcing.sh b/ush/prepare_forcing.sh new file mode 100755 index 0000000..e05d3aa --- /dev/null +++ b/ush/prepare_forcing.sh @@ -0,0 +1,76 @@ +#!/bin/bash + +# --------------------------------------------------------------------------- # +# # +# Launches scripts to retrieve and process wind, currents, ice, and # +# waterlevel for an rwps forecast. Command line inputs are date cyc and # +# meshID. For example, run: # +# # +# $sh prepare_forcing 20260904 00 oc_500m_10km # +# # +# to generate forcing for mesh rwps.oc_500m_10km.msh # +# # +# Last Changed : 09-04-2026 Sep 2026 # +# --------------------------------------------------------------------------- # + +echo 'setting paths...' + +export PDY=$1 +export cyc=$2 +export meshID=$3 + +readonly HOMErwps=$(cd "$(dirname "$(readlink -f -n "${BASH_SOURCE[0]}")")" && git rev-parse --show-toplevel) +cd "${HOMErwps}/sorc" || exit 1 + +source "${HOMErwps}/ush/detect_machine.sh" +source "${HOMErwps}/ush/module-setup.sh" +#source "${HOMErwps}/versions/build.ver" + +export MACHINE_ID +export HOMErwps + +if [[ -z "${MACHINE_ID}" ]]; then + echo "FATAL: Unable to determine target machine" + exit 1 +fi + + +export mesh="$HOMErwps/fix/rwps.$meshID.msh" +export fix="$HOMErwps/fix" +export prep="$HOMErwps/PrepInputs" +export tmp="$prep/tmpfiles" +export frc="$prep/forcing" +export outdir=$prep + +echo $HOMErwps +echo $fix $prep $tmp $outdir +mkdir -p $prep +mkdir -p $tmp +mkdir -p $frc + +#machine dependend path to rtofs, nbm, rrfs, and stofs forecast files +export COMINrtofs="/lfs/h1/ops/prod/com/rtofs/v2.5/rtofs.$PDY/" +export COMINnbm="/lfs/h3/mdl/ptmp/mdl.nbm/blend/v5.2/blend.$PDY/$cyc/grib2" +export COMINrrfs="/lfs/h1/ops/prod/com/rrfs/v1.0/rrfs.$PDY/$cyc" +export COMINstofs="/lfs/h1/ops/prod/com/stofs/v3.1/stofs_2d_glo.$PDY" + +export COMINlocal=$tmp + +#machine dependend path to RWPS fix files +export RWPSfix=/lfs/h2/emc/couple/noscrub/keston.smith/RWPS + +# copy mesh to local fix directory +cp -p $RWPSfix/fix/$meshID/20260722/rwps.$meshID.msh $HOMErwps/fix/ +# copy Interpoplation weights for nbm, rrfs, rtofs and stofs to local fix directory +cp -p $RWPSfix/fix/$meshID/20260722/InterpolationWeights*$meshID*.nc $HOMErwps/fix/ +# copy distance to boundary for nbm, rrfs, rtofs and stofs to local fix directory +cp -p $RWPSfix/fix/$meshID/20260722/DistToBndy*$meshID*.nc $HOMErwps/fix/ + +meshname="${mesh##*/}" +export meshname="${meshname: 0: -4}" + +#Retrieve current and process for forecast cycle +qsub -V $HOMErwps/ecf/jrwps_prep_current.ecf +qsub -V $HOMErwps/ecf/jrwps_prep_ice.ecf +qsub -V $HOMErwps/ecf/jrwps_prep_waterlevel.ecf +qsub -V $HOMErwps/ecf/jrwps_prep_wind.ecf diff --git a/ush/preprocess/README.txt b/ush/preprocess/README.txt new file mode 100644 index 0000000..0e3763d --- /dev/null +++ b/ush/preprocess/README.txt @@ -0,0 +1,6 @@ +Scripts for preprocessing RWPS forcing for a forecast cycle. To run: + +$cd RWPS/sorc +$../prepare_forcing.sh 20260825 00 oc_1500m_30km + +to prepare forcing for the RWPS forecast. Currently setup to work on wcoss2. diff --git a/ush/preprocess/add_err_var_to_file.py b/ush/preprocess/add_err_var_to_file.py new file mode 100644 index 0000000..0e8545d --- /dev/null +++ b/ush/preprocess/add_err_var_to_file.py @@ -0,0 +1,134 @@ +import numpy as np +import netCDF4 as nc +import sys +import interp_utilities as iutil +import os +#Add compute and add error variance field to existing file + +UseUnixTime=True +nargin = len(sys.argv) - 1 + +flin=sys.argv[1] +dist2bnd_file=sys.argv[2] + +VarParam0=sys.argv[3] +VarParam=VarParam0.split(":") + +data=nc.Dataset(dist2bnd_file,"r") +dist2bnd=np.array(data["dist2bnd"][:]) +zi=np.array(data["depth"][:]) +nn=len(dist2bnd) +data.close() + +#remove path and file +suffixp = flin.rfind(".") +dirp = flin.rfind("/") +filename=flin[dirp+1:suffixp] + + +VariableType="WaterLevel" +if "wind" in filename: + VariableType="Wind" +elif "uv" in filename: + VariableType="Current" +elif "vel" in filename: + VariableType="Current" +elif "ice" in filename: + VariableType="Ice" +if "water" in filename: + VariableType="WaterLevel" +if "level" in filename: + VariableType="WaterLevel" + +print("flin = "+flin+ ", VariableType = "+ VariableType) +if VariableType=="Current": + VarShallow=float(VarParam[0]) # variance (m/s)**2 for shallow regions + VarDeep=float(VarParam[1]) # variance (m/s)**2 for deep regions + BatShallow=float(VarParam[2]) # isobath (m) for shallow regions + BatDeep=float(VarParam[3]) # isobath (m) for deep regions + if "stofs" in flin: +# Variance = iutil.VarianceLinearDepth(zi,1.,100.,50.,250.) + Variance = iutil.VarianceLinearDepth (zi, VarShallow, VarDeep, BatShallow, BatDeep) + if "rtofs" in flin: #variance high in shallows and near boundary of coverage +# VarianceDepth = iutil.VarianceLinearDepth(zi,100.,1.,50.,250.) + VarianceDepth = iutil.VarianceLinearDepth(zi,VarShallow,VarDeep,BatShallow,BatDeep) + VarLambda= float(VarParam[4]) # lengthscale (km) for linear transition from bounadry variance(==VarShallow) to interior variance(==VarDeep) + VarianceBnd = iutil.VarianceLinearDistanceToBndy( dist2bnd, VarDeep,VarShallow, VarLambda) + Variance = np.maximum(VarianceDepth, VarianceBnd) + +if VariableType=="WaterLevel": + if (("stofs" in flin) or (True)): +# Variance = 1.+0*zi + VarInterior=float(VarParam[0]) # variance (m)**2 for stofs water level + Variance = VarInterior+0.*zi + +if VariableType=="Wind": + ##LocalFS = [ rwps_pr, rwps_hi, rwps_ak, rwps_conus, rwps_na] # file names + ##VarFS = [ 4. , 4. , 9. , 16. , 25. ] # (m m /s /s) + ##LambdaFS = [ 150. , 200. , 500. , 1000. , 1500. ] # (km) + VarInterior=float(VarParam[0]) # variance (m/s)**2 for interior of forecast + if "nbm" in flin: + Variance = VarInterior + np.zeros(nn) + if "rrfs" in flin: + VarBoundary = float(VarParam[1]) # variance (m/s)**2 for boundary of forecast + VarLambda = float(VarParam[2]) # lengthscale (km) for linear transition from bounadry variance to interior variance + Variance = iutil.VarianceLinearDistanceToBndy( dist2bnd, VarInterior, VarBoundary,VarLambda ) + + +if VariableType=="Ice": + ##LocalFS = [ rwps_pr, rwps_hi, rwps_ak, rwps_conus, rwps_na] # file names + ##VarFS = [ 4. , 4. , 9. , 16. , 25. ] # (m m /s /s) + ##LambdaFS = [ 150. , 200. , 500. , 1000. , 1500. ] # (km) + VarInterior=float(VarParam[0]) # variance (m/s)**2 for interior of forecast + if "rtofs" in flin: + Variance = VarInterior + np.zeros(nn) + if "nbm" in flin: + VarBoundary = float(VarParam[1]) # variance (m/s)**2 for boundary of forecast + VarLambda = float(VarParam[2]) # lengthscale (km) for linear transition from bounadry variance to interior variance + Variance = iutil.VarianceLinearDistanceToBndy( dist2bnd, VarInterior, VarBoundary,VarLambda ) + + + + + + +fltmp=flin+".tmp.nc" +try: + os.remove(fltmp) +except: + print("creating "+fltmp+" temporarily") + +data0 = nc.Dataset(flin,"r") + +with nc.Dataset(fltmp, "w", format="NETCDF4") as ncout: + # 1. Copy Global Attributes + ncout.setncatts({attr: data0.getncattr(attr) for attr in data0.ncattrs()}) + # 2. Copy Dimensions + for name, dimension in data0.dimensions.items(): + # If the dimension is unlimited, pass None to createDimension + dim_len = len(dimension) if not dimension.isunlimited() else None + ncout.createDimension(name, dim_len) + for name, src_var in data0.variables.items(): + dst_var = ncout.createVariable(name, src_var.datatype, src_var.dimensions) + dst_var.setncatts({attr: src_var.getncattr(attr) for attr in src_var.ncattrs()}) + dst_var[:] = src_var[:] + + time=np.asarray(data0["time"][:]) + nt=len(time) + ErrorVariance=np.zeros((nt,nn)) + for k in range(nt): + ErrorVariance[k,:]=Variance[:] + + if not 'node' in ncout.dimensions: + ncadd.createDimension('node' , nn) + if not 'time' in ncout.dimensions: + ncadd.createDimension('time' , nt) + if not 'ErrorVariance' in ncout.variables: + ErrorVariance_var=ncout.createVariable('ErrorVariance', 'f8', ('time','node')) + ErrorVariance_var.long_name = 'forecast error variance' + ErrorVariance_var.units = "(field units)**2" + ErrorVariance_var.standard_name = 'errror variance' + ErrorVariance_var[:]=ErrorVariance + + ncout.close +os.rename(fltmp, flin) diff --git a/ush/preprocess/add_mesh_geom_to_file.py b/ush/preprocess/add_mesh_geom_to_file.py new file mode 100644 index 0000000..c06deda --- /dev/null +++ b/ush/preprocess/add_mesh_geom_to_file.py @@ -0,0 +1,73 @@ +import numpy as np +import netCDF4 as nc +import sys +import interp_utilities as iutil +import os + +UseUnixTime=True +nargin = len(sys.argv) - 1 + +flin=sys.argv[1] +mshfl=sys.argv[2] + +meshslash=mshfl.rfind('/')+1 +meshname=mshfl[meshslash:len(mshfl)-3] + +xi, yi, ei, zi = iutil.loadWW3Mesh(mshfl) +nn=len(xi) +ne=ei.shape[0] + +fltmp=flin+".tmp.nc" +try: + os.remove(fltmp) +except: + print("creating "+fltmp+" temporarily") + +data0 = nc.Dataset(flin,"r") + +with nc.Dataset(fltmp, "w", format="NETCDF4") as ncout: + # 1. Copy Global Attributes + ncout.setncatts({attr: data0.getncattr(attr) for attr in data0.ncattrs()}) + # 2. Copy Dimensions + for name, dimension in data0.dimensions.items(): + # If the dimension is unlimited, pass None to createDimension + dim_len = len(dimension) if not dimension.isunlimited() else None + ncout.createDimension(name, dim_len) + for name, src_var in data0.variables.items(): + dst_var = ncout.createVariable(name, src_var.datatype, src_var.dimensions) + dst_var.setncatts({attr: src_var.getncattr(attr) for attr in src_var.ncattrs()}) + dst_var[:] = src_var[:] + + if not 'node' in data0.dimensions: + ncout.createDimension('node' , nn) + if not 'element' in data0.dimensions: + ncout.createDimension('element' , ne) + if not 'noel' in data0.dimensions: + ncout.createDimension('noel', 3) + ncout.meshname=meshname + ncout.mesh=mshfl + + if not 'longitude' in ncout.variables: + lon_var=ncout.createVariable('longitude', 'f8', ('node',)) + lon_var.units = 'degree_east' + lon_var.long_name = 'longitude' + lon_var.standard_name = 'longitude' + lon_var.axis = 'X' + lon_var[:]=xi[:] + + if not 'latitude' in ncout.variables: + lat_var=ncout.createVariable('latitude', 'f8', ('node',)) + lat_var.units = 'degree_north' + lat_var.long_name = 'latitude' + lat_var.standard_name = 'latitude' + lat_var.axis = 'Y' + lat_var[:]=yi[:] + + if not 'tri' in ncout.variables: + tri_var=ncout.createVariable('tri', 'i4', ('noel','element')) + tri_var.long_name = 'element list' + tri_var.standard_name = 'element list' + tri_var[:]=np.transpose(ei) + + ncout.close +os.rename(fltmp, flin) diff --git a/ush/preprocess/bayes_forecast_update.py b/ush/preprocess/bayes_forecast_update.py new file mode 100644 index 0000000..fcc04c7 --- /dev/null +++ b/ush/preprocess/bayes_forecast_update.py @@ -0,0 +1,109 @@ +import datetime +import numpy as np +import netCDF4 as nc +import sys +import re +import interp_utilities as iutil +flin0=sys.argv[1] +flin1=sys.argv[2] +flout=sys.argv[3] + +varname0=sys.argv[4] +varname=varname0.split(":") +nvar=len(varname) + +data0 = nc.Dataset(flin0,"r") +t=np.asarray(data0["time"][:]) +var=np.asarray(data0["ErrorVariance"][:,:]) + +x=np.asarray(data0["longitude"][:]) +y=np.asarray(data0["latitude"][:]) +e=np.asarray(data0["tri"][:,:]) + +nt=len(t) +nn=len(x) +ne=e.shape[1] +noel=e.shape[0] + +nvar=len(varname) +field=np.zeros((nvar,nt,nn)) + +nan=float("nan") +for jv in range(nvar): + tmp=np.asarray(data0[varname[jv]][:,:]) + if "_FillValue" in data0[varname[jv]].ncattrs(): + fill_value0=data0[varname[jv]]._FillValue + tmp[np.where(tmp==fill_value0)]=nan + field[jv,:,:]=tmp + +data1 = nc.Dataset(flin1,"r") +t1=np.asarray(data1["time"][:]) +var1=np.asarray(data1["ErrorVariance"][:,:]) +nt1,nn1=var1.shape +field1=np.zeros((nvar,nt1,nn1)) +for jv in range(nvar): + tmp=np.asarray(data1[varname[jv]][:,:]) + if "_FillValue" in data1[varname[jv]].ncattrs(): + fill_value1=data1[varname[jv]]._FillValue + tmp[np.where(tmp==fill_value1)]=nan + field1[jv,:,:]=tmp + +Inf=float('inf') +um1=field1[0,0,:] +ng1=np.where( um1**2>=0 ) # find points that are valid floats for field1 + +for k in range(nt): + #j=np.where(t[k]==t1) # find common merge point in data + j=np.where(np.abs(t[k]-t1)<120.) # find common merge point in data + j=j[0].tolist() + if len(j)==1: + print("j="+str(j)+" : k="+str(k)) + for jv in range(nvar): + field[jv,k,ng1] = ( field[jv,k,ng1] + + ( var[k,ng1] / ( var[k,ng1] + var1[j,ng1] ) ) + * (field1[jv,j,ng1]-field[jv,k,ng1]) + ) + #update error variance to posterior + var[k,ng1]=var[k,ng1] * ( var1[j,ng1] / ( var[k,ng1]+var1[j,ng1] ) ) + +with nc.Dataset(flout, 'w', format='NETCDF4') as ncout: + ncout.createDimension('level' , 1) + ncout.createDimension('node' , nn) + ncout.createDimension('element' , ne) + ncout.createDimension('time', nt) + ncout.createDimension('noel', 3) + + lon_var=ncout.createVariable('longitude', 'f8', ('node',)) + lon_var.units = 'degree_east' + lon_var.long_name = 'longitude' + lon_var.standard_name = 'longitude' + lon_var.axis = 'X' + lon_var[:]=x[:] + + lat_var=ncout.createVariable('latitude', 'f8', ('node',)) + lat_var.units = 'degree_north' + lat_var.long_name = 'latitude' + lat_var.standard_name = 'latitude' + lat_var.axis = 'Y' + lat_var[:]=y[:] + + time_var=ncout.createVariable('time', 'f4', ('time',)) + iutil.CopyAttributes(data0["time"], time_var) + time_var[:]=t[:] + + tri_var=ncout.createVariable('tri', 'i4', ('noel','element')) + tri_var.long_name = 'element list' + tri_var.standard_name = 'element list' + tri_var[:,:]=e + + #Copy attributes from old file to new + for jv in range(nvar): + f_var=ncout.createVariable(varname[jv], 'f4', ('time','node')) + iutil.CopyAttributes(data0[varname[jv]], f_var) + f_var[:,:]=field[jv,:,:] + + ErrVar0=data0["ErrorVariance"] + ErrVar_var=ncout.createVariable("ErrorVariance", 'f4', ('time','node')) + iutil.CopyAttributes(ErrVar0, ErrVar_var) + ErrVar_var[:,:]=var[:,:] + ncout.close diff --git a/ush/preprocess/convert_weights_to_netcdf.py b/ush/preprocess/convert_weights_to_netcdf.py new file mode 100644 index 0000000..7d1ee50 --- /dev/null +++ b/ush/preprocess/convert_weights_to_netcdf.py @@ -0,0 +1,134 @@ +import numpy as np +import netCDF4 as nc +import sys +import interp_utilities as iutil +import scipy.sparse as sp + +nargin = len(sys.argv) - 1 + +flin=sys.argv[1] +mshfl=sys.argv[2] +meshslash=mshfl.rfind('/')+1 + +AddExtrapolationSupport=True + +TextWeightFl=sys.argv[3] + +flout=TextWeightFl[0:len(TextWeightFl)-4]+".nc" + +xi, yi, ei, zi = iutil.loadWW3Mesh(mshfl) +nn_dst=len(xi) +data = nc.Dataset(flin,"r") +#read spaital dimensions and determine if input mesh is curvilinear or regular +x=np.asarray(data["x"][:]) +#convert to RWPS coordinates + +jEast=np.where(x>90) +x[jEast]=x[jEast]-360. + +nn_src=len(x) + +#Read in weights from cat of text output(rows are not in order) +F=np.loadtxt(TextWeightFl) + +node_dst=F[:,0] # row number (destination node) +ele_num_src=F[:,1] # Source element number (not used) +node_src=F[:,2:5] # Source node +DistToCenterKM=F[:,5] # distance of destination node from source element center (not used) +weights=F[:,6:9] # interpolation weights from source nodes to destination node + +node_dst=node_dst.astype(int) +node_src=node_src.astype(int) +ele_num_src=ele_num_src.astype(int) + +nnz=F.shape[0] +n_s=3*nnz +row=np.zeros(n_s,dtype=int) +col=np.zeros(n_s,dtype=int) +val=np.zeros(n_s) + +n=0 +for k in range(nnz): + for j in range(3): + row[n]=node_dst[k] + col[n]=node_src[k,j] + val[n]=weights[k,j] + n=n+1 +print("n="+str(n)) +print("n_s="+str(n_s)) +n_s=n + +with nc.Dataset(flout, 'w', format='NETCDF4') as ncout: + ncout.createDimension('n_s' , n_s) + ncout.setncattr("Nrows", nn_dst) + ncout.setncattr("Ncols", nn_src) + ncout.setncattr("SrcFieldType", "unstructured") + + r_var=ncout.createVariable('row', 'i4', ('n_s',)) + r_var.long_name = 'row index' + r_var[:]=row[:] + + c_var=ncout.createVariable('col', 'i4', ('n_s',)) + c_var.long_name = 'column index' + c_var[:]=col[:] + + s_var=ncout.createVariable('S', 'f8', ('n_s',)) + s_var.long_name = 'matrix value' + s_var[:]=val[:] + +#Consider adding (x,y) destination and (x,y) for source. This is usefull for extrapolation + if AddExtrapolationSupport: + x=np.asarray(data["x"][:]) + y=np.asarray(data["y"][:]) + ncout.createDimension('nn_src' , nn_src) + ncout.createDimension('nn_dst' , nn_dst) + + xsrc_var=ncout.createVariable('x_src', 'f8', ('nn_src',)) + xsrc_var.long_name = 'interpolation source node longitude' + xsrc_var[:]=x[:] + + ysrc_var=ncout.createVariable('y_src', 'f8', ('nn_src',)) + ysrc_var.long_name = 'interpolation source node latitude' + ysrc_var[:]=y[:] + + xdst_var=ncout.createVariable('x_dst', 'f8', ('nn_dst',)) + xdst_var.long_name = 'interpolation destination node longitude' + xdst_var[:]=xi[:] + + ydst_var=ncout.createVariable('y_dst', 'f8', ('nn_dst',)) + ydst_var.long_name = 'interpolation destination node latitude' + ydst_var[:]=yi[:] + + zdst_var=ncout.createVariable('depth_dst', 'f8', ('nn_dst',)) + zdst_var.long_name = 'interpolation destination node latitude' + zdst_var[:]=zi[:] + +Extrapolate=sys.argv[4] + +if Extrapolate: + dist2bnd=np.full(len(xi), np.inf) #all points are inside boundary- No boundary with this type of extrapolation +else: + matrix = sp.coo_matrix((weights, (row-1, col-1)), shape=(nn_src,nn_dst)).tocsr() + row_sum = matrix.sum(axis=1) + j0=np.where( row_sum==0 ) # destination nodes with no coverage from interpolation matrix + j0=np.array(j0[0]).tolist() + u0=np.ones(xi.shape) + nan=float("nan") + u0[j0]=nan + dist2bnd=iutil.CalculateDistanceToInterpEnvelope(xi,yi,u0, 1.) + +dist2bnd_file = "DistToBndy."+mshfl[meshslash:len(mshfl)-4]+".stofs.nc" +with nc.Dataset(dist2bnd_file, 'w', format='NETCDF4') as ncout: + ncout.createDimension('node' , nn_dst) + d_var=ncout.createVariable('dist2bnd', 'f4', ('node',)) + d_var.long_name = 'distance to boundary' + d_var.units = 'km' + d_var.standard_name = 'distance to boundary' + d_var[:]=dist2bnd[:] + + z_var=ncout.createVariable('depth', 'f4', ('node',)) + z_var.long_name = 'mesh depth' + z_var.units = 'm' + z_var.standard_name = 'depth' + z_var[:]=zi[:] + diff --git a/ush/preprocess/get_current.sh b/ush/preprocess/get_current.sh new file mode 100644 index 0000000..ac51b5d --- /dev/null +++ b/ush/preprocess/get_current.sh @@ -0,0 +1,8 @@ +#!/bin/bash + +# This script retrieves global RTOFS and STOFS as netcdf files +cd $HOMErwps/ush/preprocess + +rtofs/get_rtofs.sh $PDY & +stofs/get_stofs.sh $PDY $cyc current& +wait; diff --git a/ush/preprocess/get_ice.sh b/ush/preprocess/get_ice.sh new file mode 100644 index 0000000..c99dc27 --- /dev/null +++ b/ush/preprocess/get_ice.sh @@ -0,0 +1,9 @@ +#!/bin/bash + +# This script retrieves ice forecasts for global RTOFS and Alaska NBM as netcdf files + +cd $HOMErwps/ush/preprocess + +sh rtofs/get_rtofs_ice.sh $PDY $cyc & +sh nbm/get_nbm_ice.sh $PDY $cyc & +wait; diff --git a/ush/preprocess/get_waterlevel.sh b/ush/preprocess/get_waterlevel.sh new file mode 100644 index 0000000..26b099e --- /dev/null +++ b/ush/preprocess/get_waterlevel.sh @@ -0,0 +1,6 @@ +#!/bin/bash + +## This script retrieves global STOFS water level as netcdf file + +cd $HOMErwps/ush/preprocess +stofs/get_stofs.sh $PDY $cyc waterlevel diff --git a/ush/preprocess/get_wind.sh b/ush/preprocess/get_wind.sh new file mode 100644 index 0000000..056d3c5 --- /dev/null +++ b/ush/preprocess/get_wind.sh @@ -0,0 +1,48 @@ +#!/bin/bash + +# This script retrieves rrfs and nbm winbds and exports as +# netcdf files + +cd $HOMErwps/ush/preprocess + +echo "retrieving winds from rrfs and nbm for rwps wind" +( + echo "sh nbm/make_nbm_wind.sh $PDY $cyc oc > $tmp/nbm.oc.out" + nbm/make_nbm_wind.sh $PDY $cyc oc > $tmp/nbm.oc.out + echo "retrieved winds from nbm oc domain" + echo "Not retrieving other nbm domain winds" +)& + +( + echo "sh rrfs/make_rrfs_wind.sh $PDY $cyc na > $tmp/rrfs.na.out" + rrfs/make_rrfs_wind.sh $PDY $cyc na > $tmp/rrfs.na.out + echo "retrieved winds from rrfs na domain" +)& + +( + echo "sh rrfs/make_rrfs_wind.sh $PDY $cyc ak > $tmp/rrfs.ak.out" + rrfs/make_rrfs_wind.sh $PDY $cyc ak > $tmp/rrfs.ak.out + echo "retrieved winds from rrfs ak domain" +)& + +( + echo "sh rrfs/make_rrfs_wind.sh $PDY $cyc pr > rrfs.pr.out" + rrfs/make_rrfs_wind.sh $PDY $cyc pr > $tmp/rrfs.pr.out + echo "retrieved winds from rrfs pr domain" +)& + +( + echo "sh rrfs/make_rrfs_wind.sh $PDY $cyc hi > rrfs.hi.out" + rrfs/make_rrfs_wind.sh $PDY $cyc hi > $tmp/rrfs.hi.out + echo "retrieved winds from rrfs hi domain" +)& + +( + echo "sh rrfs/make_rrfs_wind.sh $PDY $cyc conus > rrfs.conus.out" + rrfs/make_rrfs_wind.sh $PDY $cyc conus > $tmp/rrfs.conus.out + echo "retrieved winds from rrfs conus domain" +)& +wait + +echo "finished retrieving winds from rrfs and nbm for rwps wind" + diff --git a/ush/preprocess/interp_time.py b/ush/preprocess/interp_time.py new file mode 100644 index 0000000..e64c4b9 --- /dev/null +++ b/ush/preprocess/interp_time.py @@ -0,0 +1,156 @@ +import os +import argparse +import numpy as np +import netCDF4 as nc +import sys +import math + +from scipy.interpolate import interp1d + +nargin = len(sys.argv) - 1 + +flin=sys.argv[1] +flinNewTimes=sys.argv[2] +flout=sys.argv[3] + +varname0=sys.argv[4] +varname=varname0.split(":") +nvar=len(varname) + +# If InterpSecondFile is True interpolate the values from the second file are interpolated to the common time points +# and extrapolate (via persistance) for all times beyond the range of time from the second file + +InterpSecondFile=False +if nargin>4: + InterpSecondFile=eval(sys.argv[5]) + print("Interpolating fields from second file: "+str(InterpSecondFile)+" and extrapolating via persistance in time") + +#might want to rewrite with an input file to derive full set of times +print("running interp_time.py: takes a forecast (arg1="+flin+") and interpolates to times ") +print("in another file (arg2="+flinNewTimes+"). The origonal and interpoated values are output to ") +print("a new file(arg3="+flout+") .") +print("interpolating for variables:") +for jv in range(nvar): + print(varname[jv]) + +data0 = nc.Dataset(flin,"r") +t=np.asarray(data0["time"][:]) + +data1 = nc.Dataset(flinNewTimes,"r") +t1=np.asarray(data1["time"][:]) + +#make sorted union of times in both forecasts +tf = np.union1d(t, t1) + +#eliminate times occuring outside the range of arg1 to avoid extrapolation +j=np.where( np.logical_and( tf>=np.min(t) , tf<=np.max(t) ) ) +print(j) +j=j[0].tolist() +tf=tf[j] + +print("Initial times:") +print(t) +print("New times:") +print(tf) + +nt=len(tf) + +x=np.asarray(data0["longitude"][:]) +y=np.asarray(data0["latitude"][:]) +e=np.asarray(data0["tri"][:,:]) +print("e") +print(e) +print(e.shape) +nn=len(x) +ne=e.shape[1] +noel=e.shape[0] + +InterpolatedVariables=np.zeros((nvar,nt,nn)) + +IsExtrap=False +if InterpSecondFile: + t=t1 + data=data1 + IsExtrap=True +else: + data=data0 + +ntt=len(t) +#unique_vals, indices = np.unique(arr, return_index=True) +tu,indx=np.unique(t,return_index=True) +for jv in range(nvar): + u=np.asarray(data[varname[jv]][:,:]) + fill_value0 = data[varname[jv]]._FillValue + nan=float('nan') + jb=np.where(u==fill_value0) + u[jb]=nan + +#set up interpolator for u + if IsExtrap: + fi = interp1d(tu, u[indx,:], axis=0, kind='linear',fill_value="extrapolate") #extrapolated values will be overwritten + else: + fi = interp1d(tu, u[indx,:], axis=0, kind='linear') + uf=fi(tf) + + #re-insert initial values at times that match initial time points + print("Interpolation compleate: now re-insert initial values at times that match initial time points") + print("to remove small interpolation artifacts") + for k in range(nt): + j=np.where(tf[k]==t) + j=j[0].tolist() + print("new time index: "+str(k)+" is same as old time index:"+str(j)) + if len(j)>0: + print("mapping exactly back to origonal values at time: "+str(tf[k])) + j=j[0] + uf[k,:]=u[j,:] + + if IsExtrap: + jExtrapEarly=np.where(tft[ntt-1]) + for jxtrp in range(len(jExtrapLate)): + uf[jExtrapLate[jxtrp],:]=u[ntt-1,:] + + InterpolatedVariables[jv,:,:]=uf[:,:] + +import interp_utilities as iutil +with nc.Dataset(flout, 'w', format='NETCDF4') as ncout: + + ncout.createDimension('level' , 1) + ncout.createDimension('node' , nn) + ncout.createDimension('element' , ne) + ncout.createDimension('time', nt) + ncout.createDimension('noel', 3) + + lon_var=ncout.createVariable('longitude', 'f8', ('node',)) + lon_var.units = 'degree_east' + lon_var.long_name = 'longitude' + lon_var.standard_name = 'longitude' + lon_var.axis = 'X' + lon_var[:]=x[:] + + lat_var=ncout.createVariable('latitude', 'f8', ('node',)) + lat_var.units = 'degree_north' + lat_var.long_name = 'latitude' + lat_var.standard_name = 'latitude' + lat_var.axis = 'Y' + lat_var[:]=y[:] + + time_var0=data0["time"] + time_var=ncout.createVariable('time', 'f8', ('time',)) + iutil.CopyAttributes(time_var0, time_var) + time_var[:]=tf[:] + + tri_var=ncout.createVariable('tri', 'i4', ('noel','element',)) + tri_var.long_name = 'element list' + tri_var.standard_name = 'element list' + tri_var[:,:]=e + + #Copy attributes from old file to new + for jv in range(nvar): + var=data0[varname[jv]] + f_var=ncout.createVariable(varname[jv], 'f4', ('time','node')) + iutil.CopyAttributes(var, f_var) + f_var[:,:]=InterpolatedVariables[jv,:,:] + ncout.close diff --git a/ush/preprocess/interp_utilities.py b/ush/preprocess/interp_utilities.py new file mode 100644 index 0000000..909a8f9 --- /dev/null +++ b/ush/preprocess/interp_utilities.py @@ -0,0 +1,764 @@ +from scipy.interpolate import RegularGridInterpolator +from scipy.sparse import csr_matrix +from scipy.sparse.csgraph import connected_components + +from datetime import datetime +import numpy as np +import netCDF4 as nc +import sys +import re + +#Convert Time to "seconds since 1970-01-01 00:00:00.0" +#eg 'seconds since 2024-04-04 12:00:00 ! NCDASE - BASE_DAT' +def ConvertTimeToUnixTime(flin,TimeVarName = None): + if TimeVarName == None: + TimeVarName="time" + data = nc.Dataset(flin,"r") + #print("TimeVarName = "+TimeVarName) + timevar=data[TimeVarName] + time=np.asarray(data[TimeVarName][:]) + epoch_1970 = datetime(1970, 1, 1, 0, 0, 0) + TimeUnitsString=timevar.units + TimeUnitsStrings=TimeUnitsString.split(" ") + tunits=TimeUnitsStrings[0] + dstr=TimeUnitsStrings[2] + dstr=dstr.split("-") + tstr=TimeUnitsStrings[3] + tstr=tstr.split(":") + print(dstr) + print(tstr) + secstr=tstr[2].rsplit(".", 1)[0] + base_date = datetime(int(dstr[0]),int(dstr[1]),int(dstr[2]),int(tstr[0]),int(tstr[1]),int(secstr)) + base_offset = int((base_date - epoch_1970).total_seconds()) + if tunits=="seconds": + unix_time = time + base_offset + if tunits=="days": + unix_time = time*24*60*60 + base_offset + if tunits=="hours": + unix_time = time*60*60 + base_offset + return unix_time + +def FileNameToUnixTime(flin,FcastPDY,FcastCYC): + year0=int(FcastPDY[0:4]) + month0=int(FcastPDY[4:6]) + day0=int(FcastPDY[6:8]) + hr0=int(FcastCYC) + #remove path and file + suffixp = flin.rfind(".") + dirp = flin.rfind("/") + flin=flin[dirp+1:suffixp] + + IsForecast=True + ntimep=flin.find(".f") + ntimeu=flin.find("_f") + ntime=max(ntimep,ntimeu) + print("ntime: "+str(ntime)) + if ntime<1: + IsForecast=False + ntimep=flin.find(".n") + ntimeu=flin.find("_n") + ntime=max(ntimep,ntimeu) + + print("ntime: "+str(ntime)) + ctime=flin[ntime+2:ntime+5] + print(ctime) + ctime=ctime.replace(".", "") #remove trailing "." in some file names + + print("FileNameToUnixTime A:") + print(flin) + print(ctime) + hrf=int(ctime) + epoch_1970 = datetime(1970, 1, 1, 0, 0, 0) + FileTime = datetime(year0,month0,day0, hr0, 0, 0) + base_offset = int((FileTime - epoch_1970).total_seconds()) + if IsForecast: + unix_time = base_offset + abs(hrf)*3600 + else: + unix_time = base_offset - abs(hrf)*3600 + + print("FileNameToUnixTime:") + print(str(year0)+" "+str(month0)+ " " +str(day0) + " " +str(hr0)) + print(flin) + print(ctime) + print(base_offset) + print(hrf) + return unix_time + +def loadWW3Mesh(fl): + print("mesh file="+fl) + f=open(fl, 'r') + header = f.readline() + header = f.readline() + header = f.readline() + header = f.readline() + header = f.readline() # number of nodes + nn=int(header) + print("nn = "+str(nn)) + xi=np.zeros(nn) + yi=np.zeros(nn) + zi=np.zeros(nn) + k=0 + for i in range(nn): + A = f.readline() + B=A.lstrip() + values = B.split(" ") +# print(values) + if len(values)>5: + xi[k]=values[2] + yi[k]=values[4] + zi[k]=values[6] + else: + xi[k]=values[1] + yi[k]=values[2] + zi[k]=values[3] + k=k+1 + print("number of nodes read: "+str(k)) + header = f.readline() + header = f.readline() + header = f.readline() # number of elements + ne=int(header)#includes boundary nodes and actual elements + print("ne="+str(ne)+" -includes boundary nodes") + nbnd=0 + bnd=[] + eix=np.zeros((ne,3), dtype=int) + k=0 + for i in range(ne): + A = f.readline() + #print(A) + values = A.split(" ") + #print(values) + if len(values) == 6: + if int(values[2])==2: + bnd.append(int(values[5])) + nbnd=nbnd+1 + if len(values)>15: + eix[k,0]=int(values[12]) + eix[k,1]=int(values[14]) + eix[k,2]=int(values[16]) + k=k+1 + elif len(values)>7: + eix[k,0]=int(values[6]) + eix[k,1]=int(values[7]) + eix[k,2]=int(values[8]) + k=k+1 + ei=eix[range(k),:] + print("number of open boundary nodes read: "+str(nbnd)) + print("number of elements read: "+str(k)) + return xi, yi, ei, zi + + +import numpy as np +from scipy.interpolate import griddata +#import matplotlib.pyplot as plt + +def interpolate_curvilinear_to_points(lon_in, lat_in, data_in, lon_out, lat_out): + """ + Performs bilinear-equivalent interpolation from a curvilinear grid to new points. + + Args: + lon_in (np.ndarray): 2D array of input longitudes. + lat_in (np.ndarray): 2D array of input latitudes. + data_in (np.ndarray): 2D array of data values corresponding to (lon_in, lat_in). + lon_out (np.ndarray or list): Longitudes of the target points. + lat_out (np.ndarray or list): Latitudes of the target points. + + Returns: + np.ndarray: Interpolated data values at the target points. + """ + # Flatten the input coordinates and data into 1D arrays + # griddata expects points as a list of (x, y) tuples or a 2D array + points_in = np.vstack((lon_in.flatten(), lat_in.flatten())).T + values_in = data_in.flatten() + + # Define the target points + points_out = np.vstack((lon_out, lat_out)).T + + # Perform the interpolation using scipy.interpolate.griddata with 'linear' method + # The 'linear' method in griddata is the appropriate choice for curvilinear data + # as it uses triangulation. + interpolated_data = griddata(points_in, values_in, points_out, method='linear') + + return interpolated_data + + +def interpolate_curvilinear_to_pointsMD(lon_in, lat_in, data_in, lon_out, lat_out): + """ + Performs bilinear-equivalent interpolation from a curvilinear grid to new points. + + Args: + lon_in (np.ndarray): 2D array of input longitudes. + lat_in (np.ndarray): 2D array of input latitudes. + data_in (np.ndarray): 3D array of data values corresponding to (lon_in, lat_in, ntimes). + lon_out (np.ndarray or list): Longitudes of the target points. + lat_out (np.ndarray or list): Latitudes of the target points. + + Returns: + np.ndarray: Interpolated data values at the target points (length(lon_out/lat_out x ntimes)). + """ + # Flatten the input coordinates and data into 1D arrays + # griddata expects points as a list of (x, y) tuples or a 2D array + points_in = np.vstack((lon_in.flatten(), lat_in.flatten())).T + + shp=data_in.shape + print(shp) + nx=shp[0] + ny=shp[1] + nt=shp[2] + + ns=nx*ny + S=np.zeros((ns,nt)) + s0=np.zeros((nx,ny)) + for k in range(nt): + s0[:,:]=np.transpose(data_in[:,:,k]) + S[:,k] = s0.flatten() + + # Define the target points + points_out = np.vstack((lon_out, lat_out)).T + + # Perform the interpolation using scipy.interpolate.griddata with 'linear' method + # The 'linear' method in griddata is the appropriate choice for curvilinear data + # as it uses triangulation. + interpolated_data = griddata(points_in, S, points_out, method='linear') + + return interpolated_data + + +def interpolate_curvilinear_to_pointsRRFS(lon_in, lat_in, data_in, lon_out, lat_out): + """ + Performs bilinear-equivalent interpolation from a curvilinear grid to new points. + + Args: + lon_in (np.ndarray): 2D array of input longitudes. + lat_in (np.ndarray): 2D array of input latitudes. + data_in (np.ndarray): 3D array of data values corresponding to (lon_in, lat_in, ntimes). + lon_out (np.ndarray or list): Longitudes of the target points. + lat_out (np.ndarray or list): Latitudes of the target points. + + Returns: + np.ndarray: Interpolated data values at the target points (length(lon_out/lat_out x ntimes)). + """ + # Flatten the input coordinates and data into 1D arrays + # griddata expects points as a list of (x, y) tuples or a 2D array + + points_in = np.vstack((lon_in.flatten(), lat_in.flatten())).T + # points_in = np.hstack((lon_in.flatten(), lat_in.flatten())) + + shp=data_in.shape + print(shp) + nx=shp[0] + ny=shp[1] + nt=shp[2] + + ns=nx*ny + S=np.zeros((ns,nt)) + s0=np.zeros((ny,nx)) + for k in range(nt): + s0[:,:]=np.transpose(data_in[:,:,k]) + S[:,k] = s0.flatten() + + # Define the target points + points_out = np.vstack((lon_out, lat_out)).T + + # Perform the interpolation using scipy.interpolate.griddata with 'linear' method + # The 'linear' method in griddata is the appropriate choice for curvilinear data + # as it uses triangulation. + print(S.shape) + print(points_in) + print(points_out) + + interpolated_data = griddata(points_in, S, points_out, method='linear') + + return interpolated_data + +def CopyAttributes(VarOld, VarNew): + #Copy attributes from old NetCDF file variable to new NetCDF file variable + att_names = VarOld.ncattrs() + for jatt in range(len(att_names)): + att_name=att_names[jatt] + if (not (att_name=="_FillValue")): + att_value = VarOld.getncattr(att_name) + VarNew.setncattr(att_name, att_value) + return + + +############################################################################################ +# BEGIN WIND TO RWPS INTERP ROUTINES +import esmpy +import scipy.sparse as sp + +def CurvilinearGridCreateInterpWeights(xi,yi,x1,y1, weights_file): +# Compute interpolation weights to interpolate from curvilinear grid (x1,y1) to points (xi,yi) +# and store in netcdf file using ESMPY + nx,ny=x1.shape + nn=len(xi) + n1=nx*ny + src_lon = x1 + src_lat = y1 + dst_lon=np.zeros((1,nn)) + dst_lat=np.zeros((1,nn)) + print("nn="+str(nn)) + dst_lon[0,:] = xi[:] + dst_lat[0,:] = yi[:] + + esmpy.Manager() + + src_grid = esmpy.Grid( + max_index=np.array([src_lon.shape[0], src_lon.shape[1]]), + staggerloc=esmpy.StaggerLoc.CENTER, + coord_sys=esmpy.CoordSys.SPH_DEG + ) + + dst_grid = esmpy.Grid( + max_index=np.array([nn, 1]), + staggerloc=esmpy.StaggerLoc.CENTER, + coord_sys=esmpy.CoordSys.SPH_DEG + ) + + # 4. Populate the coordinate data + src_lon_ptr = src_grid.get_coords(0) + src_lat_ptr = src_grid.get_coords(1) + src_lon_ptr[...] = src_lon + src_lat_ptr[...] = src_lat + + dst_lon_ptr = dst_grid.get_coords(0) + dst_lat_ptr = dst_grid.get_coords(1) + dst_lon_ptr[...] = dst_lon.T + dst_lat_ptr[...] = dst_lat.T + + # 5. Create esmpy Fields + src_field = esmpy.Field(src_grid, name="src_field") + dst_field = esmpy.Field(dst_grid, name="dst_field") + src_field.data[...]=np.sqrt(np.abs(x1/180))/(90+y1) # arbitrary function of x,y + + np.savetxt('F.txt', src_field.data[...]) + np.savetxt('X.txt', x1) + np.savetxt('Y.txt', y1) + print(f"Creating weights: {weights_file}") + regrid = esmpy.Regrid( + src_field, + dst_field, + filename=weights_file, + regrid_method=esmpy.RegridMethod.BILINEAR, + ignore_degenerate=True, # <--- Add this parameter + unmapped_action=esmpy.UnmappedAction.IGNORE # Optional: Ignores missing/masked points + ) +#Add number of rows and columns to weights file for clarity when constructing sparse matrix for interpolation + with nc.Dataset(weights_file, mode="a") as ds: + ds.Nrows = nn + ds.Ncols = n1 + + np.savetxt('Fi.txt', dst_field.data[...]) + np.savetxt('xi.txt', xi) + np.savetxt('yi.txt', yi) + + return + +def CalculateDistanceToBoundary(xi,yi,x1,y1): +# Distance to boundary calculation for use when interpolation envelope corresponds with +# interior of curvilinear grid boundary. +# +# Inputs: +# xi (nn): longitude of unstructured mesh nodes +# yi (nn): latitude of unstructured mesh nodes +# x1 (nx x ny): longitude for interpoltated field +# y1 (nx x ny): latitude for interpoltated field +# +# Outputs: +# dist2bnd (nn) : distance to edge of interpolation envelope. + nx=x1.shape[0] + ny=x1.shape[1] + xb=np.hstack((x1[1,:],x1[:,ny-1].T,x1[nx-1,:],x1[:,1].T)) + yb=np.hstack((y1[1,:],y1[:,ny-1].T,y1[nx-1,:],y1[:,1].T)) + np.savetxt('xbyb.txt', np.vstack((xb,yb))) + dist2bnd=np.zeros(nn) + for k in range(nn): + dist2bnd[k]=QuickDistance(yi[k],xi[k],yb,xb) + if k%10000==0: + print("calculating distance to boundary, "+str(k)+":"+ str(nn)+":"+str(k/nn) ) + return dist2bnd + +def CalculateDistanceToInterpEnvelope(xi,yi,fi,SearchWidth): +# Alternative distance to boundary calculation for use when interpolation envelope is +# distinctly interior to curvilinear grid boundary as happens for RRFS NA grid +# +# Inputs: +# xi (nn): longitude of unstructured mesh nodes +# yi (nn): latitude of unstructured mesh nodes +# fi (nn): interpolated field on mesh nodes with 'nan' values outside of interpolation envelope +# SearchWidth: Computational speed up to remove extra search points in distance to boundary +# +# Outputs: +# dist2bnd (nn) : distance to edge of interpolation envelope. dist2bnd[k]=0 if (xi[k],yi[k]) is outside +# of the interpolation envelope + nn=len(xi) + #dist2bnd=np.zeros(nn) + dist2bnd=np.full(nn,np.nan) + jin = np.where(~np.isnan(fi))[0].tolist()#points inside interpolation envelope + jout = np.where(np.isnan(fi))[0].tolist() #points outside interpolation envelope + xin=xi[jin] + yin=yi[jin] + xout=xi[jout] + yout=yi[jout] + jxU=np.where( xout < np.max(xin)+SearchWidth )[0].tolist() + jxD=np.where( xout > np.min(xin)-SearchWidth )[0].tolist() + jyU=np.where( yout < np.max(yin)+SearchWidth )[0].tolist() + jyD=np.where( yout > np.min(yin)-SearchWidth )[0].tolist() + j=list( set(jxU) & set(jxD) & set(jyU) & set(jyD) ) + xout=xout[j] + yout=yout[j] + din=np.zeros(len(jin)) + print(len(xin)) + print(len(xout)) + for k in range(len(xin)): + din[k]=QuickDistance(yin[k],xin[k],yout,xout) # distance from node to closest point not interpolated to + if k%10000==0: + print("calculating distance to boundary, "+str(k)+":"+ str(nn)+":"+str(k/nn) ) + dist2bnd[jin]=din + return dist2bnd + +def QuickDistance(lat1, lon1, lats2, lons2): + deg2kmY=111. + deg2kmX=np.cos( np.pi * lat1 / 180.)*deg2kmY + d= np.min( np.sqrt( ( (lat1-lats2)*deg2kmY)**2 + ((lon1-lons2)*deg2kmX)**2 ) ) + return d + +def VarianceLinearDistanceToBndy(DistanceToBoundary, InteriorVariance, VarianceOnBoundary, LengthScale): + InteriorNodeList=np.where(DistanceToBoundary**2 >= 0 ) + Variance=np.zeros(len(DistanceToBoundary))+np.inf + SpatialFunction=DistanceToBoundary/LengthScale + j=np.where(SpatialFunction>1.) + SpatialFunction[j]=1. + Variance[InteriorNodeList] = VarianceOnBoundary + ( InteriorVariance - VarianceOnBoundary ) * SpatialFunction[InteriorNodeList] + return Variance + +def VarianceInverseDistanceToBndy( DistanceToBoundary, InteriorVariance, LengthScale): + InteriorNodeList=np.where(DistanceToBoundary**2 >= 0 ) + Variance=np.zeros(len(DistanceToBoundary))+np.inf + SpatialFunction=LengthScale / DistanceToBoundary + j=np.where(SpatialFunction>1.) + SpatialFunction[j]=1. + Variance[InteriorNodeList] = InteriorVariance * SpatialFunction[InteriorNodeList] + return Variance + +def VarianceLinearDepth(zi,VarianceShallow,VarianceDeep,Zshallow,Zdeep): + Variance = VarianceShallow + (VarianceDeep-VarianceShallow)*(zi-Zshallow)/(Zdeep-Zshallow) + js=np.where(ziZdeep) + Variance[js]=VarianceShallow + Variance[jd]=VarianceDeep + return Variance + +import numpy as np +import netCDF4 as nc +def WriteInterpolationWeightsToNetCDF(weights_file,row,col,weights,Nrows,Ncols): + #create a esmpy style sparse matrix netcdf file + print(Nrows) + print(Ncols) + print(row) + print(col) + print(weights) + n_s=len(weights) + + if isinstance(row , list): + row = np.array( row ) + if isinstance( col, list): + col = np.array( col ) + if isinstance( weights, list): + weights = np.array( weights ) + + with nc.Dataset(weights_file, 'w', format='NETCDF4') as ncout: + + ncout.createDimension('n_s' , n_s) + ncout.setncattr("Nrows", Nrows) + ncout.setncattr("Ncols", Ncols) + + r_var=ncout.createVariable('row', 'i4', ('n_s',)) + r_var.long_name = 'row index' + r_var[:]=row[:] + + c_var=ncout.createVariable('col', 'i4', ('n_s',)) + c_var.long_name = 'column index' + c_var[:]=col[:] + + s_var=ncout.createVariable('S', 'f4', ('n_s',)) + s_var.long_name = 'matrix value' + s_var[:]=weights[:] + + + +def IsInElement(x,y,xp,yp): + IsIn=False + c1 = (x[1] - x[0]) * (yp - y[0]) - (y[1] - y[0]) * (xp - x[0]) + c2 = (x[2] - x[1]) * (yp - y[1]) - (y[2] - y[1]) * (xp - x[1]) + c3 = (x[0] - x[2]) * (yp - y[2]) - (y[0] - y[2]) * (xp - x[2]) + if ( ( c1 > 0 and c2 > 0 and c3 > 0) or ( c1 < 0 and c2 < 0 and c3 < 0) ): + IsIn=True + if (c1*c2*c3 == 0): # include points on triangle + IsIn=True + return IsIn + +def FindElement(x,y,e,xi,yi): + nd=len(e.shape) + e=np.squeeze(e) + j=-9999 + ne=e.shape[0] + if nd > 1: + xc = np.squeeze(np.mean(x[e], axis=1)) + yc = np.squeeze(np.mean(y[e], axis=1)) + DistanceToElements = np.abs((xi + 1j*yi) - (xc + 1j*yc)) + n=0 + while (j < 0 and n < ne) : + n=n+1 + j = np.argmin( DistanceToElements ) + xl = np.squeeze(x[e[j,:]]) + yl = np.squeeze(y[e[j,:]]) + IsIn=IsInElement(xl,yl,xi,yi) + if IsIn : + return j + else: + DistanceToElements[j]=float('inf') + j=-9999 + else: + xl = np.squeeze(x[e]) + yl = np.squeeze(y[e]) + IsIn=IsInElement(xl,yl,xi,yi) + if IsIn: + j=0 + return j + else: + j=-9999 + return j + + return j + + +def compute_mesh_to_mesh_interp_weights(x, y, e, xi, yi): + """ + Compute mesh-to-mesh interpolation weights using barycentric coordinates. + + Parameters: + ----------- + x : array-like + X coordinates of source mesh nodes + y : array-like + Y coordinates of source mesh nodes + e : array-like + Element connectivity matrix (n_elements x 3 for triangular elements) + (nodes indexed from 1 to nn) + xi : array-like + X coordinates of target points + yi : array-like + Y coordinates of target points + + Returns: + -------- + weights : ndarray + Interpolation weights (n_points x 3) + nodes : ndarray + Node indices for each point (n_points x 3) + elenum : ndarray + Element number for each point (n_points,) + Dist2EleCenter : ndarray + Distance to element center for each element + """ + + # Convert to numpy arrays + + deg2kmY=111. + UseNearestEle = False #if True just use closest element center + N = 12 # search N nearest elements (nearest by element center to target node distance) + + print(xi) + x = np.asarray(x) + y = np.asarray(y) + e = np.asarray(e) + xi = np.asarray(xi) + yi = np.asarray(yi) + print("e.shape") + print(e.shape) + # convert node indexs to 0 .. nn-1 + e0=e-1 + + xc = np.mean(x[e0], axis=1) + yc = np.mean(y[e0], axis=1) + + # Initialize output arrays + n_points = len(xi) + n_elements = len(e) + weights = np.zeros((n_points, 3)) + nodes = np.zeros((n_points, 3), dtype=int) + elenum = np.zeros(n_points, dtype=int) + Dist2EleCenter = np.zeros(n_points) + +# UseNearestEle=True + for k in range(n_points): + x0=xi[k] + y0=yi[k] + deg2kmX=np.cos( np.pi * y0 / 180.)*deg2kmY + distances = np.abs( deg2kmX*(x0 - xc) + 1j*deg2kmY*(y0 - yc)) + if UseNearestEle: + jg=np.argmin(distances) + else: + raw_indices = np.argpartition(distances, N)[:N] + sorted_sub_indices = np.argsort(distances[raw_indices]) + jg = raw_indices[sorted_sub_indices] + if UseNearestEle : + j=jg + else: + j=FindElement(x,y,e0[jg,:],x0,y0) + if j < 0: + print("couldn't find element for point : "+ str(k)) + print(str(x0)+" : "+str(y0)) + if len(jg)>0: + j = jg[0] + else: + j=0 + print("using closest element at distance : "+ str(distances[j])) + else: + j=jg[j] + + elenum[k] = j + Dist2EleCenter[k] = distances[j] + xl = x[e0[j, :]] + yl = y[e0[j, :]] + + # Compute barycentric coordinates using shoelace formula for areas + # Area of triangle formed by vertices 1, 2, and point (a3) + xt = np.array([xl[0], xl[1], x0, xl[0]]) + yt = np.array([yl[0], yl[1], y0, yl[0]]) + a3 = -np.dot(xt[1:4] - xt[0:3], yt[0:3] + yt[1:4]) / 2 + # Area of triangle formed by vertices 3, 1, and point (a2) + xt = np.array([xl[2], xl[0], x0, xl[2]]) + yt = np.array([yl[2], yl[0], y0, yl[2]]) + a2 = -np.dot(xt[1:4] - xt[0:3], yt[0:3] + yt[1:4]) / 2 + # Area of triangle formed by point, vertices 2, 3 (a1) + xt = np.array([x0, xl[1], xl[2], x0]) + yt = np.array([y0, yl[1], yl[2], y0]) + a1 = -np.dot(xt[1:4] - xt[0:3], yt[0:3] + yt[1:4]) / 2 + + # Normalize to get barycentric weights + total_area = a1 + a2 + a3 + weights[k, :] = [a1, a2, a3] / total_area + + nodes[k, :] = e[j, :] #<- indexed 1 .. nn + # Progress reporting every 100 iterations + if (k + 1) % 100 == 0: + print(f"Progress: {k+1}/{n_points}") + return weights, nodes, elenum, Dist2EleCenter + + +def InterpolateField2Nodes(nodes,weights, f): + fi=np.zeros(weights.shape[0]) + for k in range(weights.shape[0]): + fl=f[nodes[k,:]] + wl=weights[k,:] + fi[k]=np.dot(wl,fl) + if (not np.abs(fi[k]) > 0.): + fi[k]=0. + if ( fi[k]>np.max(fl) ): + fi[k]=np.max(fl) + if ( fi[k]cd $PBS_O_WORKDIR +#-->echo "Running on node: $(hostname)" +#--> python my_script.py --task $PBS_ARRAY_INDEX + +import os +def WriteInterpJobscriptPBS(fl,flin,mshfl,Njobs, ComputeNodes): + + meshslash=mshfl.rfind('/')+1 + TmpOutDir="STOFSInterpWeights."+mshfl[meshslash:len(mshfl)-4] + WghtFl="STOFS.wght."+mshfl[meshslash:len(mshfl)-4]+".txt" + WghtFl="InterpWeights."+mshfl[meshslash:len(mshfl)-4]+".stofs.txt" + WghtFlNetCDF="InterpWeights."+mshfl[meshslash:len(mshfl)-4]+".stofs.nc" + Njobs=32 #OVERWRITE FOR NOW + with open(fl, 'w') as f: + f.write("#PBS -N ESMPy\n") + f.write("#PBS -j oe\n") + f.write("#PBS -S /bin/bash\n") + f.write("#PBS -q dev\n") + f.write("#PBS -A NWPS-DEV\n") + f.write("#PBS -l walltime=01:00:00\n") + f.write("#PBS -J 1-"+str(Njobs)+"\n") + f.write("#PBS -l select=2:ncpus=32:mem=128gb\n") +# f.write("#PBS -l select=1:ncpus=1:mem=8G\n") + f.write("#PBS -l place=excl\n") + f.write("#PBS -l debug=true\n") + f.write("#PBS -r y\n") + + f.write("module reset\n") + f.write("module load PrgEnv-intel/8.5.0\n") + f.write("module load intel/19.1.3.304\n") + f.write("module load craype/2.7.17\n") + f.write("module load cray-mpich/8.1.19\n") + f.write("module load hdf5-C/1.14.0\n") + f.write("module load netcdf-C/4.9.2\n") + f.write("module load esmf-C/8.6.0\n") + f.write("module load ve/hafs/2.1\n") + + f.write("pip list -v\n") + + current_dir = os.getcwd() + f.write("cd "+current_dir+"\n") + + f.write("# calculate interpolation weights in parallel geographically \n") + f.write("python compute_unstr_to_rwps_interp_weights.py "+flin+" "+mshfl+" $PBS_ARRAY_INDEX " + str(Njobs)+" > InterpJob.$PBS_ARRAY_INDEX.out \n") + + +def WriteInterpJobscriptSLURM(fl,flin,mshfl,Njobs, ComputeNodes): + + meshslash=mshfl.rfind('/')+1 + TmpOutDir="STOFSInterpWeights."+mshfl[meshslash:len(mshfl)-4] + WghtFl="STOFS.wght."+mshfl[meshslash:len(mshfl)-4]+".txt" + WghtFlNetCDF="STOFS.wght."+mshfl[meshslash:len(mshfl)-4]+".nc" + with open(fl, 'w') as f: + f.write("#!/bin/bash \n") + f.write("#SBATCH --job-name=STOFS_interp_masterscript \n") + f.write("#SBATCH --ntasks=1 \n") # ntasks per interpolation + f.write("#SBATCH --time=08:00:00 \n") + f.write("#SBATCH --output=mpi_test_%j.log \n") + f.write("#SBATCH --error=%j.err \n") + f.write("#SBATCH --account=marine-cpu \n") + f.write("#SBATCH --nodes="+str(ComputeNodes)+" \n") + f.write("#SBATCH --ntasks-per-core=1"+" \n") + f.write("#SBATCH --array=0-"+str(Njobs-1)+" \n") + + f.write(" \n") + + f.write("module purge \n") + f.write("module use /scratch4/NCEPDEV/marine/Ali.Salimi/Hera_Data/HR4-OPT/FromJessica/Keston/ICunstructuredRuns15km-implicit-450s/global-workflow/sorc/ufs_model.fd/modulefiles \n") + f.write("module load ufs_ursa.intel \n") + f.write("module load py-scipy/1.14.1 \n") + f.write("module load py-netcdf4/1.7.1.post2 \n") + f.write("pip list \n") + f.write("# calculate interpolation weights in parallel geographically \n") + f.write("srun python GeoSubsetInterpolateSTOFS.py "+flin+" "+mshfl+" $SLURM_ARRAY_TASK_ID " + str(Njobs)+" > InterpJob.$SLURM_ARRAY_TASK_ID.out \n") + f.write("wait\n") + f.write("# concatonate different parts of the mesh to common text file \n") + f.write("cat "+TmpOutDir+"/Part.IntrpWghts.*.txt > "+WghtFl+" \n") + f.write("# convert output weights to netcdf file \n") + f.write("python convert_weights_to_netcdf.py "+flin+" "+mshfl+" \n") + + + flintrp="STOFS.to."+mshfl[meshslash:len(mshfl)-4]+".sh" + flout=flin[0:-2]+mshfl[meshslash:len(mshfl)-4]+".nc" + flinuv=flin[0:-3]+".vel.nc" + floutuv=flinuv[0:-2]+mshfl[meshslash:len(mshfl)-4]+".nc" + with open(flintrp, 'w') as f: + f.write("#!/bin/bash \n") + f.write("#SBATCH --job-name=STOFS_interp_masterscript \n") + f.write(" \n") + f.write("module purge \n") + f.write("module use /scratch4/NCEPDEV/marine/Ali.Salimi/Hera_Data/HR4-OPT/FromJessica/Keston/ICunstructuredRuns15km-implicit-450s/global-workflow/sorc/ufs_model.fd/modulefiles \n") + f.write("module load ufs_ursa.intel \n") + f.write("module load py-scipy/1.14.1 \n") + f.write("module load py-netcdf4/1.7.1.post2 \n") + f.write("pip list \n") + f.write("# calculate interpolation weights in parallel geographically \n") + f.write("python InterpolateSTOFS.py "+flin+" "+mshfl+" "+flout+" zeta 2\n") + f.write("python InterpolateSTOFS.py "+flinuv+" "+mshfl+" "+floutuv+" u-vel:v-vel 2\n") + diff --git a/ush/preprocess/interpolate_with_weights.py b/ush/preprocess/interpolate_with_weights.py new file mode 100644 index 0000000..0c49572 --- /dev/null +++ b/ush/preprocess/interpolate_with_weights.py @@ -0,0 +1,238 @@ +import numpy as np +import netCDF4 as nc +import sys +import interp_utilities as iutil +import xarray as xr +import scipy.sparse as sp +from scipy.interpolate import NearestNDInterpolator +import datetime + +# Engine for interpolating to WW3 unstructured mesh using precomputed interpolation weights from netcdf files with forecasts +# +# to call: +# python InterpolateSTOFS.py input_file meshpath outputfile variable1:variable2:variable3 ExtrapMethod +# +# example: +# python InterpolateSTOFS.py stofs.20260608.00/stofs.cwl.vel.nc meshes/RWPS.V0a.small.msh tesdtoZ.vel.nc u-vel:v-vel 2 +# or: +# python InterpolateSTOFS.py stofs.20260608.00/stofs.cwl.nc meshes/RWPS.V0a.small.msh tesdtoZ.vel.nc zeta 1 +# +# ExtrapMethod =-1 no extrapolation, NaN's potentially in output where source field is dry +# ExtrapMethod = 0 NaN values in interpolated field replaced with 0.0 +# ExtrapMethod = 1 Nearest Neighbor extrapolation from valid source values +# ExtrapMethod = 2 Nearest Neighbor extrapolation from valid interpolated values +# ExtrapMethod = 3 Nearest Neighbor extrapolation from interpolated nodes which allways have valid values (faster than 2 for larger source mesh) + +UseUnixTime=True +nargin = len(sys.argv) - 1 + +flin=sys.argv[1] + +weights_file=sys.argv[2] + +flout=sys.argv[3] +varname0=sys.argv[4] +varname=varname0.split(":") + +ExtrapMethod=-1 # no extrapolation +if nargin>4: + ExtrapMethod=int(sys.argv[5]) + +if ExtrapMethod==-1: + print("no extrapolation, nan left in place in output") +if ExtrapMethod==0: + print("Fill missing values in interpolated field with value 0") +if ExtrapMethod==1: + print("extrapolation from nearest valid point in source- can be slow if source mesh is much larger than destination mesh") +if ExtrapMethod==2: + print("extrapolation from nearest valid point in destination (interpolated field)") + +with xr.open_dataset(weights_file) as ds_s: + # Standard sparse storage uses 'row', 'col', and 'data' variables + row = ds_s['row'].values + col = ds_s['col'].values + weights = ds_s['S'].values + Nrows=ds_s.attrs.get('Nrows') + Ncols=ds_s.attrs.get('Ncols') + SrcFieldType=ds_s.attrs.get('SrcFieldType') + if ExtrapMethod>0: #these extrapolation methods need the source and destination nodes + x=ds_s['x_src'].values + y=ds_s['y_src'].values + xi=ds_s['x_dst'].values + yi=ds_s['y_dst'].values + +nni=Nrows +n1=Ncols +print("row") +print(row) +print("col") +print(col) + +matrix = sp.coo_matrix((weights, (row-1, col-1)), shape=(Nrows,Ncols)).tocsr() +print("sparse interpolation matrix") +print(matrix) +row_sum = matrix.sum(axis=1) +j0=np.where( row_sum==0 ) # destination nodes with no coverage from interpolation matrix + +data = nc.Dataset(flin,"r") +if "time" in data.variables: + time=iutil.ConvertTimeToUnixTime(flin,"time") +elif "MT" in data.variables: + time=iutil.ConvertTimeToUnixTime(flin,"MT") +else: + print("No time variable found in "+flin+" EXITING") + sys.exit(1) + +nt=len(time) + +print(time) + +nvar=len(varname) +vari=np.zeros((nvar,nt,nni)) + +if ExtrapMethod>=0: + IsExtrap=np.zeros((nvar,nt,nni),dtype=int) + +if ExtrapMethod==3: + AnyExtrap=np.zeros((nvar,nni),dtype=int) + +nan=float("nan") +for jv in range(nvar): + +# fill_value0=data[varname[jv]]._FillValue + try: + fill_value0=data[varname[jv]]._FillValue + except: + fill_value0=-99999 + print("fill value="+str(fill_value0)) + for k in range(nt): + print("interpolating for time step = "+str(k)+" of "+str(nt)) + vshp = data.variables[varname[jv]].shape +# if SrcFieldType=="unstructured": + if len(vshp)==1: + var=np.asarray(data[varname[jv]][:]) # No time dimension?, just spatial data to interpolate + elif len(vshp)==2: + var=np.asarray(data[varname[jv]][k,:]) + elif len(vshp)==3: # Wind field with dimensions time, x, y + if "wind" in flin: + var0=np.asarray(data[varname[jv]][k,:,:]) + print("var0.shape") + print(var0.shape) + print(n1) + var=np.transpose(var0).reshape(n1) + elif "ice" in flin: + var0=np.asarray(data[varname[jv]][k,:,:]) + if "rtofs" in flin: #remove bad geometry edges + var0=var0[1:-1,1:-1] + var=np.transpose(var0).reshape(n1) + + elif len(vshp)==4: # RTOFS field with 2nd dimensional "Level" and garbage boundries + var0=np.asarray(data[varname[jv]][k,0,:,:]) + if "rtofs" in flin: #remove bad geometry edges + var0=var0[1:-1,1:-1] + var=np.transpose(var0).reshape(n1) + + else: + print(vshp) + print(len(vshp)) + print("unkown data shape for "+varname[jv]+" terminating") + sys.exit() + + #replace fill with nan to avoid interpolating fill + j=np.where(var==fill_value0) + var[j]=nan + if ExtrapMethod==0: + j=np.where(np.isnan(var)) + var[j]=0. + + vari[jv,k,:] = matrix @ var # actual spatial interpolation step + vari[jv,k,j0]=nan # empty rows + if ExtrapMethod==3: # Fast posthoc nearest neighbor extrapolator + jd=np.where(np.isnan(vari[jv,k,:])) + AnyExtrap[jv,jd]=1. + elif ExtrapMethod>0:# and ExtrapMethod<3: + jd=np.where(np.isnan(vari[jv,k,:])) + dstp=np.array((xi[jd],yi[jd])) + if ExtrapMethod==1: + #extrapolate using nearest neighbor of source with valid value + js=np.where(~np.isnan(var)) + srcp=np.array((x[js],y[js])) + srcv=var[js] + if ExtrapMethod==2: + #extrapolate using nearest neighbor of interpolated field with valid value + js=np.where(~np.isnan(vari[jv,k,:])) + srcp=np.array((xi[js],yi[js])) + tmp=vari[jv,k,js] + srcv=tmp.flatten() + interp = NearestNDInterpolator(srcp.T,srcv) + ExtrapVals = interp( dstp.T ) + vari[jv,k,jd]=ExtrapVals + IsExtrap[jv,k,jd]=1 + +if ExtrapMethod==0: + jd=np.where(np.isnan(vari)) + vari[jd]==0. + IsExtrap[jd]=1 + +if ExtrapMethod==3: #posthoc extrapolation from points which are valid at all times + for jv in range(nvar): + jd=np.where(AnyExtrap[jv,:]==1) # nodes that have some "nan" intrepolated values + js=np.where(AnyExtrap[jv,:]==0) # nodes that have no "nan" intrepolated values + srcp = np.array((xi[js],yi[js])).T + srcv = vari[jv,0,js] #dummy input field + dstp = np.array((xi[jd],yi[jd])).T + srcv=srcv[0,:] + interpolator = NearestNDInterpolator(srcp, srcv) + distances, jsrc = interpolator.tree.query(dstp) + jd=jd[0] + for k in range(nt): + jdk=np.where(np.isnan(vari[jv,k,jd])) + jdk=jdk[0] + vari[jv,k,jd[jdk]]=vari[jv,k,jsrc[jdk]] # value of nearest "always valid" destination point + IsExtrap[jv,k,jd[jdk]]=1 + +print("nn(target mesh) = "+str(nni)+": Nrows = "+str(Nrows)) +print("nn(source mesh) = "+str(n1)+": Ncols = "+str(Ncols)) +if not ((nni==Nrows) and (n1==Ncols)): + print("WARNING: Wrong matrix weights: number of rows from "+ mshfl +" = "+str(nni)+ + " but number of rows in "+ weights_file +" = "+str(Nrows)+ + ", number of spatial points in "+ flin +" = "+str(n1)+ + " but number of columns in "+ weights_file +" = "+str(Ncols) ) + print(" You may need to regnerate file "+ weights_file +" with appropriate weights") + +with nc.Dataset(flout, 'w', format='NETCDF4') as ncout: + + ncout.createDimension('level' , 1) + ncout.createDimension('node' , nni) + ncout.createDimension('time', nt) + + time_var=ncout.createVariable('time', 'f8', ('time',)) + varin = data["time"] + iutil.CopyAttributes(varin, time_var) + if UseUnixTime: + time_var.units = 'seconds since 1970-01-01 00:00:00.0 0:00' + time_var.standard_name = 'time' + time_var[:]=time[:] + + for jv in range(nvar): + print("writing output for :"+varname[jv]) + varin = data[varname[jv]] + F_var=ncout.createVariable(varname[jv], 'f4', ('time','node'),fill_value = fill_value0) + iutil.CopyAttributes(varin, F_var) + F_var.location = 'node' + F_var[:,:] = vari[jv,:,:] + + if ExtrapMethod >= 0 : + xtrp_var=ncout.createVariable(varname[jv]+'IsExtrap', 'i1', ('time','node')) + xtrp_var.long_name = '==1 if the interpolated value extrapolated. 0 if interpolated' + xtrp_var.standard_name = 'is extrapolated' + xtrp_var.location = 'node' + if ExtrapMethod == 0: + xtrp_var.method = 'Interpolated nan values replaced with 0' + if ExtrapMethod == 1: + xtrp_var.method = 'nearest valid neighbor in source field' + if ExtrapMethod == 2: + xtrp_var.method = 'nearest valid neighbor in interpolated field' + xtrp_var[:,:] = IsExtrap[jv,:,:] + + ncout.close diff --git a/ush/preprocess/interpolate_with_weights_and_error_variance.py b/ush/preprocess/interpolate_with_weights_and_error_variance.py new file mode 100644 index 0000000..cb3c08e --- /dev/null +++ b/ush/preprocess/interpolate_with_weights_and_error_variance.py @@ -0,0 +1,329 @@ +import numpy as np +import netCDF4 as nc +import sys +import interp_utilities as iutil +import xarray as xr +import scipy.sparse as sp +from scipy.interpolate import NearestNDInterpolator +import datetime + +# Engine for interpolating to WW3 unstructured mesh using precomputed interpolation weights from netcdf files with forecasts +# +# to call: +# python InterpolateSTOFS.py input_file meshpath outputfile variable1:variable2:variable3 ExtrapMethod +# +# example: +# python InterpolateSTOFS.py stofs.20260608.00/stofs.cwl.vel.nc meshes/RWPS.V0a.small.msh tesdtoZ.vel.nc u-vel:v-vel 2 +# or: +# python InterpolateSTOFS.py stofs.20260608.00/stofs.cwl.nc meshes/RWPS.V0a.small.msh tesdtoZ.vel.nc zeta 1 +# +# ExtrapMethod =-1 no extrapolation, NaN's potentially in output where source field is dry +# ExtrapMethod = 0 NaN values in interpolated field replaced with 0.0 +# ExtrapMethod = 1 Nearest Neighbor extrapolation from valid source values +# ExtrapMethod = 2 Nearest Neighbor extrapolation from valid interpolated values +# ExtrapMethod = 3 Nearest Neighbor extrapolation from interpolated nodes which allways have valid values (faster than 2 for larger source mesh) + +UseUnixTime=True +nargin = len(sys.argv) - 1 + +flin=sys.argv[1] + +weights_file=sys.argv[2] + +flout=sys.argv[3] +varname0=sys.argv[4] +varname=varname0.split(":") +if nargin>4: + VarParam0=sys.argv[5] + IncludeErrorVariance=True +else: + print("no variance parameters provided- error variance will not be included in file "+flout) + IncludeErrorVariance=False + +VariableType='Wind' +if "vel" in varname[0]: + VariableType='Current' +if "zeta" in varname[0]: + VariableType='WaterLevel' +if "ice" in varname[0]: + VariableType='Ice' + +ExtrapMethod=-1 # no extrapolation +if nargin>4: + ExtrapMethod=int(sys.argv[5]) + +if ExtrapMethod==-1: + print("no extrapolation, nan left in place in output") +if ExtrapMethod==0: + print("Fill missing values in interpolated field with value 0") +if ExtrapMethod==1: + print("extrapolation from nearest valid point in source- can be slow if source mesh is much larger than destination mesh") +if ExtrapMethod==2: + print("extrapolation from nearest valid point in destination (interpolated field)") + +with xr.open_dataset(weights_file) as ds_s: + # Standard sparse storage uses 'row', 'col', and 'data' variables + row = ds_s['row'].values + col = ds_s['col'].values + weights = ds_s['S'].values + Nrows=ds_s.attrs.get('Nrows') + Ncols=ds_s.attrs.get('Ncols') + SrcFieldType=ds_s.attrs.get('SrcFieldType') + if ExtrapMethod>0: #these extrapolation methods need the source and destination nodes + x=ds_s['x_src'].values + y=ds_s['y_src'].values + xi=ds_s['x_dst'].values + yi=ds_s['y_dst'].values + #shift to RWPS convention + j=np.where(x>90) + x[j]=x[j]-360. + j=np.where(xi>90) + xi[j]=xi[j]-360. + +nni=Nrows +n1=Ncols + +matrix = sp.coo_matrix((weights, (row-1, col-1)), shape=(Nrows,Ncols)).tocsr() +print("sparse interpolation matrix") +print(matrix) + +data = nc.Dataset(flin,"r") +if "time" in data.variables: + time=iutil.ConvertTimeToUnixTime(flin,"time") +elif "MT" in data.variables: + time=iutil.ConvertTimeToUnixTime(flin,"MT") +else: + print("No time variable found in "+flin+" EXITING") + sys.exit(1) + +nt=len(time) + +print(time) + +nvar=len(varname) +vari=np.zeros((nvar,nt,nni)) + +if ExtrapMethod>=0: + IsExtrap=np.zeros((nvar,nt,nni),dtype=int) + +if ExtrapMethod==3: + AnyExtrap=np.zeros((nvar,nni),dtype=int) + +nan=float("nan") +for jv in range(nvar): + fill_value0=data[varname[jv]]._FillValue + print("fill value="+str(fill_value0)) + for k in range(nt): + print("interpolating for time step = "+str(k)+" of "+str(nt)) + vshp = data.variables[varname[jv]].shape +# if SrcFieldType=="unstructured": + if len(vshp)==1: + var=np.asarray(data[varname[jv]][:]) # No time dimension?, just spatial data to interpolate + if len(vshp)==2: + var=np.asarray(data[varname[jv]][k,:]) + if len(vshp)==3: # Wind field with dimensions time, x, y + var0=np.asarray(data[varname[jv]][k,:,:]) + var=np.transpose(var0).reshape(n1) + if len(vshp)==4: # RTOFS field with 2nd dimensional "Level" and garbage boundries + var0=np.asarray(data[varname[jv]][k,0,:,:]) + if "rtofs" in flin: #remove bad geometry edges + var0=var0[1:-1,1:-1] + var=np.transpose(var0).reshape(n1) + else: + print("unkown data shape for "+varname[jv]+" terminating") + sys.exit() + + #replace fill with nan to avoid interpolating fill + j=np.where(var==fill_value0) + var[j]=nan + + if ExtrapMethod==0: + j=np.where(np.isnan(var)) + var[j]=0. + + vari[jv,k,:] = matrix @ var # actual spatial interpolation step + + if ExtrapMethod==3: # Fast posthoc nearest neighbor extrapolator + jd=np.where(np.isnan(vari[jv,k,:])) + AnyExtrap[jv,jd]=1. + elif ExtrapMethod>0:# and ExtrapMethod<3: + jd=np.where(np.isnan(vari[jv,k,:])) + dstp=np.array((xi[jd],yi[jd])) + if ExtrapMethod==1: + #extrapolate using nearest neighbor of source with valid value + js=np.where(~np.isnan(var)) + srcp=np.array((x[js],y[js])) + srcv=var[js] + if ExtrapMethod==2: + #extrapolate using nearest neighbor of interpolated field with valid value + js=np.where(~np.isnan(vari[jv,k,:])) + srcp=np.array((xi[js],yi[js])) + tmp=vari[jv,k,js] + srcv=tmp.flatten() + interp = NearestNDInterpolator(srcp.T,srcv) + ExtrapVals = interp( dstp.T ) + vari[jv,k,jd]=ExtrapVals + IsExtrap[jv,k,jd]=1 + +if ExtrapMethod==0: + jd=np.where(np.isnan(vari)) + vari[jd]==0. + IsExtrap[jd]=1 + +if ExtrapMethod==3: #posthoc extrapolation from points which are valid at all times + for jv in range(nvar): + jd=np.where(AnyExtrap[jv,:]==1) # nodes that have some "nan" intrepolated values + js=np.where(AnyExtrap[jv,:]==0) # nodes that have no "nan" intrepolated values + srcp = np.array((xi[js],yi[js])).T + srcv = vari[jv,0,js] #dummy input field + dstp = np.array((xi[jd],yi[jd])).T + srcv=srcv[0,:] + interpolator = NearestNDInterpolator(srcp, srcv) + distances, jsrc = interpolator.tree.query(dstp) + jd=jd[0] + for k in range(nt): + jdk=np.where(np.isnan(vari[jv,k,jd])) + jdk=jdk[0] + vari[jv,k,jd[jdk]]=vari[jv,k,jsrc[jdk]] # value of nearest "always valid" destination point + IsExtrap[jv,k,jd[jdk]]=1 + +print("nn(target mesh) = "+str(nni)+": Nrows = "+str(Nrows)) +print("nn(source mesh) = "+str(n1)+": Ncols = "+str(Ncols)) +if not ((nni==Nrows) and (n1==Ncols)): + print("WARNING: Wrong matrix weights: number of rows from "+ mshfl +" = "+str(nni)+ + " but number of rows in "+ weights_file +" = "+str(Nrows)+ + ", number of spatial points in "+ flin +" = "+str(n1)+ + " but number of columns in "+ weights_file +" = "+str(Ncols) ) + print(" You may need to regnerate file "+ weights_file +" with appropriate weights") + +########################################################################################## +# START: CONSTRUCT PRESCRIBED ERROR COVARIANCE FOR UPDATING ESTIMATE +########################################################################################## +#grab variables for prescribing error variance +if IcludeErrorVariance: + dist2bnd_file=weights_file.replace("InterpolationWeights", "DistToBnd") + with xr.open_dataset(dist2bnd_file) as ds_s: + dist2bnd = ds_s['dist2bnd'].values + zi = ds_s['depth'].values + + if VariableType=="Current": + VarShallow=float(VarParam[0]) # variance (m/s)**2 for shallow regions + VarDeep=float(VarParam[1]) # variance (m/s)**2 for deep regions + BatShallow=float(VarParam[2]) # isobath (m) for shallow regions + BatDeep=float(VarParam[3]) # isobath (m) for deep regions + if "stofs" in flin: +# Variance = iutil.VarianceLinearDepth(zi,1.,100.,50.,250.) + Variance = iutil.VarianceLinearDepth(zi,VarShallow,VarDeep,BatShallow,BatDeep) + if "rtofs" in flin: #variance high in shallows and near boundary of coverage +# VarianceDepth = iutil.VarianceLinearDepth(zi,100.,1.,50.,250.) + VarianceDepth = iutil.VarianceLinearDepth(zi,VarShallow,VarDeep,BatShallow,BatDeep) + VarLambda= float(VarParam[4]) # lengthscale (km) for linear transition from bounadry variance(==VarShallow) to interior variance(==VarDeep) + VarianceBnd = iutil.VarianceLinearDistanceToBndy( dist2bnd, VarDeep,VarShallow, VarLambda) + Variance = np.maximum(VarianceDepth, VarianceBnd) + + if VariableType=="WaterLevel": + if "stofs" in flin: +# Variance = 1.+0*zi + VarInterior=float(VarParam[0]) # variance (m)**2 for stofs water level + Variance = VarInterior+0.*zi + if VariableType=="Wind": + ##LocalFS = [ rwps_pr, rwps_hi, rwps_ak, rwps_conus, rwps_na] # file names + ##VarFS = [ 4. , 4. , 9. , 16. , 25. ] # (m m /s /s) + ##LambdaFS = [ 150. , 200. , 500. , 1000. , 1500. ] # (km) + VarInterior=float(VarParam[0]) # variance (m/s)**2 for interior of forecast + if "nbm" in flin: + Variance = VarInterior + 0.*dist2bnd[:] + if "rrfs" in flin: + VarBoundary = float(VarParam[1]) # variance (m/s)**2 for boundary of forecast + VarLambda = float(VarParam[2]) # lengthscale (km) for linear transition from bounadry variance to interior variance + Variance = iutil.VarianceLinearDistanceToBndy( dist2bnd, VarInterior, VarBoundary,VarLambda ) + + ErrorVariance=np.zeros((nt,nni)) + for k in range(nt): + ErrorVariance[k,:]=Variance[:] + +########################################################################################## +# FINISHED: CONSTRUCT PRESCRIBED ERROR COVARIANCE FOR UPDATING ESTIMATE +########################################################################################## + +ne=ei.shape[0] +with nc.Dataset(flout, 'w', format='NETCDF4') as ncout: + + ncout.createDimension('level' , 1) + ncout.createDimension('node' , nni) + ncout.createDimension('element' , ne) + ncout.createDimension('time', nt) + ncout.createDimension('noel', 3) + + lon_var=ncout.createVariable('longitude', 'f8', ('node',)) + lon_var.units = 'degree_east' + lon_var.long_name = 'longitude' + lon_var.standard_name = 'longitude' + lon_var.axis = 'X' + lon_var[:]=xi[:] + + lat_var=ncout.createVariable('latitude', 'f8', ('node',)) + lat_var.units = 'degree_north' + lat_var.long_name = 'latitude' + lat_var.standard_name = 'latitude' + lat_var.axis = 'Y' + lat_var[:]=yi[:] + + if UseUnixTime: + units = 'seconds since 1970-01-01 00:00:00.0 0:00' + long_name = 'verification time generated by wgrib2 function verftime()' + standard_name = 'time' + else: + varin = data["time"] + units = varin.units + standard_name = varin.standard_name + long_name = varin.long_name + + time_var=ncout.createVariable('time', 'f8', ('time',)) + time_var.units = units + time_var.long_name = long_name + time_var.standard_name = standard_name + time_var[:]=time[:] + + tri_var=ncout.createVariable('tri', 'i4', ('noel','element')) + tri_var.long_name = 'element list' + tri_var.standard_name = 'element list' + tri_var[:]=np.transpose(ei) + + for jv in range(nvar): + if varname[jv]=='' + print("writing output for :"+varname[jv]) + varin = data[varname[jv]] + units = varin.units + standard_name = varin.standard_name + long_name = varin.long_name + location = 'node' + + F_var=ncout.createVariable(varname[jv], 'f4', ('time','node'),fill_value = fill_value0) + F_var.long_name = long_name + F_var.units = units + F_var.standard_name = standard_name + F_var.location = location + F_var[:,:] = vari[jv,:,:] + + if ExtrapMethod >= 0 : + xtrp_var=ncout.createVariable(varname[jv]+'IsExtrap', 'i1', ('time','node')) + xtrp_var.long_name = '==1 if the interpolated value extrapolated. 0 if interpolated' + xtrp_var.standard_name = 'is extrapolated' + xtrp_var.location = 'node' + if ExtrapMethod == 0: + xtrp_var.method = 'Interpolated nan values replaced with 0' + if ExtrapMethod == 1: + xtrp_var.method = 'nearest valid neighbor in source field' + if ExtrapMethod == 2: + xtrp_var.method = 'nearest valid neighbor in interpolated field' + xtrp_var[:,:] = IsExtrap[jv,:,:] + + if IcludeErrorVariance: + ErrorVariance_var=ncout.createVariable('ErrorVariance', 'f4', ('time','node'),fill_value = fill_value0) + ErrorVariance_var.long_name = 'forecast error variance' + ErrorVariance_var.units = "("+units+")**2" + ErrorVariance_var.standard_name = 'variance' + ErrorVariance_var[:,:]=ErrorVariance + + ncout.close diff --git a/ush/preprocess/interpolate_with_weights_with_fillvalue.py b/ush/preprocess/interpolate_with_weights_with_fillvalue.py new file mode 100644 index 0000000..fef0007 --- /dev/null +++ b/ush/preprocess/interpolate_with_weights_with_fillvalue.py @@ -0,0 +1,233 @@ +import numpy as np +import netCDF4 as nc +import sys +import interp_utilities as iutil +import xarray as xr +import scipy.sparse as sp +from scipy.interpolate import NearestNDInterpolator +import datetime + +# Engine for interpolating to WW3 unstructured mesh using precomputed interpolation weights from netcdf files with forecasts +# +# to call: +# python InterpolateSTOFS.py input_file meshpath outputfile variable1:variable2:variable3 ExtrapMethod +# +# example: +# python InterpolateSTOFS.py stofs.20260608.00/stofs.cwl.vel.nc meshes/RWPS.V0a.small.msh tesdtoZ.vel.nc u-vel:v-vel 2 +# or: +# python InterpolateSTOFS.py stofs.20260608.00/stofs.cwl.nc meshes/RWPS.V0a.small.msh tesdtoZ.vel.nc zeta 1 +# +# ExtrapMethod =-1 no extrapolation, NaN's potentially in output where source field is dry +# ExtrapMethod = 0 NaN values in interpolated field replaced with 0.0 +# ExtrapMethod = 1 Nearest Neighbor extrapolation from valid source values +# ExtrapMethod = 2 Nearest Neighbor extrapolation from valid interpolated values +# ExtrapMethod = 3 Nearest Neighbor extrapolation from interpolated nodes which allways have valid values (faster than 2 for larger source mesh) + +UseUnixTime=True +nargin = len(sys.argv) - 1 + +flin=sys.argv[1] + +weights_file=sys.argv[2] + +flout=sys.argv[3] +varname0=sys.argv[4] +varname=varname0.split(":") + +ExtrapMethod=-1 # no extrapolation +if nargin>4: + ExtrapMethod=int(sys.argv[5]) + +if ExtrapMethod==-1: + print("no extrapolation, nan left in place in output") +if ExtrapMethod==0: + print("Fill missing values in interpolated field with value 0") +if ExtrapMethod==1: + print("extrapolation from nearest valid point in source- can be slow if source mesh is much larger than destination mesh") +if ExtrapMethod==2: + print("extrapolation from nearest valid point in destination (interpolated field)") + +with xr.open_dataset(weights_file) as ds_s: + # Standard sparse storage uses 'row', 'col', and 'data' variables + row = ds_s['row'].values + col = ds_s['col'].values + weights = ds_s['S'].values + Nrows=ds_s.attrs.get('Nrows') + Ncols=ds_s.attrs.get('Ncols') + SrcFieldType=ds_s.attrs.get('SrcFieldType') + if ExtrapMethod>0: #these extrapolation methods need the source and destination nodes + x=ds_s['x_src'].values + y=ds_s['y_src'].values + xi=ds_s['x_dst'].values + yi=ds_s['y_dst'].values + +nni=Nrows +n1=Ncols +print("row") +print(row) +print("col") +print(col) + +matrix = sp.coo_matrix((weights, (row-1, col-1)), shape=(Nrows,Ncols)).tocsr() +print("sparse interpolation matrix") +print(matrix) +row_sum = matrix.sum(axis=1) +j0=np.where( row_sum==0 ) # destination nodes with no coverage from interpolation matrix + +data = nc.Dataset(flin,"r") +if "time" in data.variables: + time=iutil.ConvertTimeToUnixTime(flin,"time") +elif "MT" in data.variables: + time=iutil.ConvertTimeToUnixTime(flin,"MT") +else: + print("No time variable found in "+flin+" EXITING") + sys.exit(1) + +nt=len(time) + +print(time) + +nvar=len(varname) +vari=np.zeros((nvar,nt,nni)) + +if ExtrapMethod>=0: + IsExtrap=np.zeros((nvar,nt,nni),dtype=int) + +if ExtrapMethod==3: + AnyExtrap=np.zeros((nvar,nni),dtype=int) + +nan=float("nan") +for jv in range(nvar): + fill_value0=data[varname[jv]]._FillValue + print("fill value="+str(fill_value0)) + for k in range(nt): + print("interpolating for time step = "+str(k)+" of "+str(nt)) + vshp = data.variables[varname[jv]].shape +# if SrcFieldType=="unstructured": + if len(vshp)==1: + var=np.asarray(data[varname[jv]][:]) # No time dimension?, just spatial data to interpolate + elif len(vshp)==2: + var=np.asarray(data[varname[jv]][k,:]) + elif len(vshp)==3: # Wind field with dimensions time, x, y + if "wind" in flin: + var0=np.asarray(data[varname[jv]][k,:,:]) + print("var0.shape") + print(var0.shape) + print(n1) + var=np.transpose(var0).reshape(n1) + elif "ice" in flin: + var0=np.asarray(data[varname[jv]][k,:,:]) + if "rtofs" in flin: #remove bad geometry edges + var0=var0[1:-1,1:-1] + var=np.transpose(var0).reshape(n1) + + elif len(vshp)==4: # RTOFS field with 2nd dimensional "Level" and garbage boundries + var0=np.asarray(data[varname[jv]][k,0,:,:]) + if "rtofs" in flin: #remove bad geometry edges + var0=var0[1:-1,1:-1] + var=np.transpose(var0).reshape(n1) + + else: + print(vshp) + print(len(vshp)) + print("unkown data shape for "+varname[jv]+" terminating") + sys.exit() + + #replace fill with nan to avoid interpolating fill + j=np.where(var==fill_value0) + var[j]=nan + if ExtrapMethod==0: + j=np.where(np.isnan(var)) + var[j]=0. + + vari[jv,k,:] = matrix @ var # actual spatial interpolation step + vari[jv,k,j0]=nan # empty rows + if ExtrapMethod==3: # Fast posthoc nearest neighbor extrapolator + jd=np.where(np.isnan(vari[jv,k,:])) + AnyExtrap[jv,jd]=1. + elif ExtrapMethod>0:# and ExtrapMethod<3: + jd=np.where(np.isnan(vari[jv,k,:])) + dstp=np.array((xi[jd],yi[jd])) + if ExtrapMethod==1: + #extrapolate using nearest neighbor of source with valid value + js=np.where(~np.isnan(var)) + srcp=np.array((x[js],y[js])) + srcv=var[js] + if ExtrapMethod==2: + #extrapolate using nearest neighbor of interpolated field with valid value + js=np.where(~np.isnan(vari[jv,k,:])) + srcp=np.array((xi[js],yi[js])) + tmp=vari[jv,k,js] + srcv=tmp.flatten() + interp = NearestNDInterpolator(srcp.T,srcv) + ExtrapVals = interp( dstp.T ) + vari[jv,k,jd]=ExtrapVals + IsExtrap[jv,k,jd]=1 + +if ExtrapMethod==0: + jd=np.where(np.isnan(vari)) + vari[jd]==0. + IsExtrap[jd]=1 + +if ExtrapMethod==3: #posthoc extrapolation from points which are valid at all times + for jv in range(nvar): + jd=np.where(AnyExtrap[jv,:]==1) # nodes that have some "nan" intrepolated values + js=np.where(AnyExtrap[jv,:]==0) # nodes that have no "nan" intrepolated values + srcp = np.array((xi[js],yi[js])).T + srcv = vari[jv,0,js] #dummy input field + dstp = np.array((xi[jd],yi[jd])).T + srcv=srcv[0,:] + interpolator = NearestNDInterpolator(srcp, srcv) + distances, jsrc = interpolator.tree.query(dstp) + jd=jd[0] + for k in range(nt): + jdk=np.where(np.isnan(vari[jv,k,jd])) + jdk=jdk[0] + vari[jv,k,jd[jdk]]=vari[jv,k,jsrc[jdk]] # value of nearest "always valid" destination point + IsExtrap[jv,k,jd[jdk]]=1 + +print("nn(target mesh) = "+str(nni)+": Nrows = "+str(Nrows)) +print("nn(source mesh) = "+str(n1)+": Ncols = "+str(Ncols)) +if not ((nni==Nrows) and (n1==Ncols)): + print("WARNING: Wrong matrix weights: number of rows from "+ mshfl +" = "+str(nni)+ + " but number of rows in "+ weights_file +" = "+str(Nrows)+ + ", number of spatial points in "+ flin +" = "+str(n1)+ + " but number of columns in "+ weights_file +" = "+str(Ncols) ) + print(" You may need to regnerate file "+ weights_file +" with appropriate weights") + +with nc.Dataset(flout, 'w', format='NETCDF4') as ncout: + + ncout.createDimension('level' , 1) + ncout.createDimension('node' , nni) + ncout.createDimension('time', nt) + + time_var=ncout.createVariable('time', 'f8', ('time',)) + varin = data["time"] + iutil.CopyAttributes(varin, time_var) + if UseUnixTime: + time_var.units = 'seconds since 1970-01-01 00:00:00.0 0:00' + time_var.standard_name = 'time' + time_var[:]=time[:] + + for jv in range(nvar): + print("writing output for :"+varname[jv]) + varin = data[varname[jv]] + F_var=ncout.createVariable(varname[jv], 'f4', ('time','node'),fill_value = fill_value0) + iutil.CopyAttributes(varin, F_var) + F_var.location = 'node' + F_var[:,:] = vari[jv,:,:] + + if ExtrapMethod >= 0 : + xtrp_var=ncout.createVariable(varname[jv]+'IsExtrap', 'i1', ('time','node')) + xtrp_var.long_name = '==1 if the interpolated value extrapolated. 0 if interpolated' + xtrp_var.standard_name = 'is extrapolated' + xtrp_var.location = 'node' + if ExtrapMethod == 0: + xtrp_var.method = 'Interpolated nan values replaced with 0' + if ExtrapMethod == 1: + xtrp_var.method = 'nearest valid neighbor in source field' + if ExtrapMethod == 2: + xtrp_var.method = 'nearest valid neighbor in interpolated field' + xtrp_var[:,:] = IsExtrap[jv,:,:] + + ncout.close diff --git a/ush/preprocess/nbm/get_nbm_ice.sh b/ush/preprocess/nbm/get_nbm_ice.sh new file mode 100755 index 0000000..1bc14ad --- /dev/null +++ b/ush/preprocess/nbm/get_nbm_ice.sh @@ -0,0 +1,21 @@ +#!/bin/bash + +# This is just a place holder with notes for the time being. Notes on aws retrieval and processing are below +# This script takes rrfs grib2 forecast files, extracts 10m u and v wind +# components and outputs to netcdf. Comand line arguments are + +ICE_FILE="$COMINnbm/blend.t$2z.icec.ak.grib2" + +OUTPUT_DIR="$COMINlocal/ice.$1.$2" +OUTPUT_FILE="$OUTPUT_DIR/nbm.$1.$2.ice.ak.nc" + +mkdir -p "$OUTPUT_DIR" +# Remove existing output file to avoid mixing old data +rm -f "$OUTPUT_FILE" + +echo "writing ice from $INPUT_DIR to $OUTPUT_FILE" + +wgrib2 "$ICE_FILE" -match ":ICEC:" -netcdf "$OUTPUT_FILE" + +echo "nbm ice processing complete for forecast date $1, cycle $2, domain ak" +echo "output written to: $OUTPUT_FILE" diff --git a/ush/preprocess/nbm/make_nbm_wind.sh b/ush/preprocess/nbm/make_nbm_wind.sh new file mode 100755 index 0000000..ca031a0 --- /dev/null +++ b/ush/preprocess/nbm/make_nbm_wind.sh @@ -0,0 +1,25 @@ +#!/bin/bash + +# This script takes rrfs grib2 forecast files, extracts 10m u and v wind +# components and outputs to netcdf. Comand line arguments are + +WSPD_FILE="$COMINnbm/blend.t$2z.wspd.$3.grib2" +WDIR_FILE="$COMINnbm/blend.t$2z.wdir.$3.grib2" + +OUTPUT_DIR="wind.$1.$2" +OUTPUT_DIR="$COMINlocal/wind.$1.$2" +OUTPUT_FILE="$OUTPUT_DIR/nbm.$1.$2.wind10m.$3.nc" + +mkdir -p "$OUTPUT_DIR" + +# Remove existing output file to avoid mixing old data +rm -f "$OUTPUT_FILE" + +echo "writing 10m wind from $INPUT_DIR to $OUTPUT_FILE" + +wgrib2 "$WSPD_FILE" -match ":WIND:10 m" -netcdf "$OUTPUT_FILE" +wgrib2 "$WDIR_FILE" -match ":WDIR:10 m" -append -netcdf "$OUTPUT_FILE" + +echo "nbm processing complete for forecast date $1, cycle $2, domain $3" +echo "output written to: $OUTPUT_FILE" + diff --git a/ush/preprocess/process_current.sh b/ush/preprocess/process_current.sh new file mode 100644 index 0000000..242b8e7 --- /dev/null +++ b/ush/preprocess/process_current.sh @@ -0,0 +1,79 @@ +#!/bin/bash + +cd $HOMErwps/ush/preprocess + +meshname="${mesh##*/}" +meshname="${meshname: 0: -4}" + +stofscur="$tmp/stofs.$PDY.$cyc/stofs_2d_glo.t${cyc}z.fields.cwl.vel.nc" +rtofscur="$tmp/rtofs.$PDY.nc" + +varnames="u-vel:v-vel" + +#name of blended wind +rwps_current=$frc/$meshname.$PDY.$cyc.current.nc +#rwps_current=$frc/$meshname.$PDY.$cyc.vel.stofsxrtofs.nc +echo "outputting combined stofs and rtofs currents to $rwps_current" + +## STOFS interpolation +stofs_wghts="$fix/InterpolationWeights.$meshname.stofs.nc" +stofs_dists="$fix/DistToBndy.$meshname.stofs.nc" +stofs_rwps="$tmp/$meshname.$PDY.$cyc.vel.cwl.stofs.nc" +stofs_rwps_ti="$tmp/$meshname.$PDY.$cyc.vel.cwl.stofs.ti.nc" + + + +## RTOFS interpolation +rtofs_wghts="$fix/InterpolationWeights.$meshname.rtofs.current.nc" +rtofs_dists="$fix/DistToBndy.$meshname.rtofs.current.nc" +rtofs_rwps="$tmp/$meshname.$PDY.vel.rtofs.nc" +rtofs_rwps_ti="$tmp/$meshname.$PDY.$cyc.vel.cwl.rtofs.ti.nc" + + + +if [ ! -f "$stofs_wghts" ]; then + echo "missing stofs interpolation weights file: $stofs_wghts" + echo "compute with script compute_unstr_to_rwps_interp_weights.sh" + exit 1 +fi +if [ ! -f "$stofs_dists" ]; then + echo "missing stofs distance to boundary file: $stofs_dists" + echo "compute with script compute_unstr_to_rwps_interp_weights.sh" + exit 1 +fi + +# extrapolate with zero fill +python interpolate_with_weights.py $stofscur $stofs_wghts $stofs_rwps $varnames 0 & + +if [ ! -f "$rtofs_wghts" ]; then + echo "missing rtofs interpolation weights file: $stofs_wghts" + echo "compute with script ComputeGridToRWPSInterpWeights.py" + exit 2 +fi +if [ ! -f "$rtofs_dists" ]; then + echo "missing stofs distance to boundary file: $rtofs_dists" + echo "compute with script compute_unstr_to_rwps_interp_weights.py" + exit 2 +fi + +# no extrapolation +python interpolate_with_weights.py $rtofscur $rtofs_wghts $rtofs_rwps $varnames -1 & + +wait; + +python add_mesh_geom_to_file.py $rtofs_rwps $mesh +python add_mesh_geom_to_file.py $stofs_rwps $mesh + +python interp_time.py $stofs_rwps $rtofs_rwps $stofs_rwps_ti $varnames False & + +#interpolate from rtofs to common stofs and rtofs times within range of stofs time +#values out of range are extrapolated to assuming persistance +python interp_time.py $stofs_rwps $rtofs_rwps $rtofs_rwps_ti $varnames True & + +wait + +python add_err_var_to_file.py $rtofs_rwps_ti $rtofs_dists 100.:1.:50.:250.:50. +python add_err_var_to_file.py $stofs_rwps_ti $stofs_dists 1.:100.:50.:250. + +python bayes_forecast_update.py $stofs_rwps_ti $rtofs_rwps_ti $rwps_current $varnames + diff --git a/ush/preprocess/process_ice.sh b/ush/preprocess/process_ice.sh new file mode 100644 index 0000000..11394a1 --- /dev/null +++ b/ush/preprocess/process_ice.sh @@ -0,0 +1,93 @@ +#!/bin/bash + +cd $HOMErwps/ush/preprocess + +meshname="${mesh##*/}" +meshname="${meshname: 0: -4}" + + +inpdir=$tmp/ice.$PDY.$cyc + +nbmice=$inpdir/nbm.$PDY.$cyc.ice.ak.nc +rtofsice=$inpdir/rtofs.ice.$PDY.nc + + +rtofs_wghts="$fix/InterpolationWeights.$meshname.rtofs.ice.nc" +rtofs_dists="$fix/DistToBndy.$meshname.rtofs.ice.nc" +nbm_ak_wghts="$fix/InterpolationWeights.$meshname.nbm.ak.nc" +nbm_ak_dists="$fix/DistToBndy.$meshname.nbm.ak.nc" + +rtofs_rwps="$tmp/$meshname.$PDY.ice.rtofs.nc" +rtofs_rwps_ti="$tmp/$meshname.$PDY.$cyc.ice.rtofs.ti.nc" + +nbm_rwps="$tmp/$meshname.$PDY.$cyc.ice.nbm.ak.nc" +nbm_rwps_ti="$tmp/$meshname.$PDY.$cyc.ice.nbm.ak.ti.nc" + +varnames="ICEC_surface" + +#rwps_ice="$frc/$meshname.$PDY.$cyc.ice.rtofsxnbm.nc" +rwps_ice="$frc/$meshname.$PDY.$cyc.ice.nc" + +## NBM AK domain interpolation +## NOTE WE ARE USING THE RRFS ak WEIGHTS HERE NEED TO RELABEL AND RECOMPUTE FOR NBM AK DOMAIN +## DIFFERENCES ARE SMALL. posional distance ~200m on 3km grid roughly + +#../../fix/DistToBndy.rwps.oc_1500m_30km.nbm.ak.nc + + +if [ ! -f "$nbm_ak_wghts" ]; then + echo "missing nbm ak interpolation weights file: $nbm_ak_wghts" + echo "compute with script compute_unstr_to_rwps_interp_weights.sh" + exit 1 +fi + +if [ ! -f "$nbm_ak_dists" ]; then + echo "missing nbm ak distance to boundary file: $nbm_ak_dists" + echo "compute with script compute_unstr_to_rwps_interp_weights.sh" + exit 1 +fi + +# no extrapolation +python interpolate_with_weights.py $nbmice $nbm_ak_wghts $nbm_rwps $varnames -1 & + +## RTOFS interpolation + +if [ ! -f "$rtofs_wghts" ]; then + echo "missing rtofs interpolation weights file: $rtofs_wghts" + echo "compute with script ComputeGridToRWPSInterpWeights.py" + exit 2 +fi +if [ ! -f "$rtofs_dists" ]; then + echo "missing stofs distance to boundary file: $rtofs_dists" + echo "compute with script compute_unstr_to_rwps_interp_weights.py" + exit 2 +fi + +# extrapolate with 0 as fill +echo "$rtofsice" +python interpolate_with_weights.py $rtofsice $rtofs_wghts $rtofs_rwps $varnames 0 & + +wait; + +python add_mesh_geom_to_file.py $rtofs_rwps $mesh +python add_mesh_geom_to_file.py $nbm_rwps $mesh + +#python add_err_var_to_file.py $rtofs_rwps $rtofs_dists 100.:1.:50.:250.:50. +#python add_err_var_to_file.py $stofs_rwps $stofs_dists 1.:100.:50.:250. + +#interpolate from stofs to common stofs and rtofs times within range of stofs time +python interp_time.py $rtofs_rwps $nbm_rwps $rtofs_rwps_ti $varnames False & + +#interpolate from rtofs to common stofs and rtofs times within range of stofs time +#values out of range are extrapolated to assuming persistance +python interp_time.py $rtofs_rwps $nbm_rwps $nbm_rwps_ti $varnames True & + +wait + +#uniform variance of 100. +python add_err_var_to_file.py $rtofs_rwps_ti $rtofs_dists 100. + +#interior variance of 4., boundary variance fof 400., transition lengthscale 9. km +python add_err_var_to_file.py $nbm_rwps_ti $nbm_ak_dists 4.:400.:9. + +python bayes_forecast_update.py $rtofs_rwps_ti $nbm_rwps_ti $rwps_ice $varnames diff --git a/ush/preprocess/process_waterlevel.sh b/ush/preprocess/process_waterlevel.sh new file mode 100644 index 0000000..4179452 --- /dev/null +++ b/ush/preprocess/process_waterlevel.sh @@ -0,0 +1,32 @@ +#!/bin/bash + +cd $HOMErwps/ush/preprocess + +meshname="${mesh##*/}" +meshname="${meshname: 0: -4}" + +stofslev="$tmp/stofs.$PDY.$cyc/stofs_2d_glo.t${cyc}z.fields.cwl.nc" +stofs_wghts="$fix/InterpolationWeights.$meshname.stofs.nc" +stofs_dists="$fix/DistToBndy.$meshname.stofs.nc" +stofs_rwps="$frc/$meshname.$PDY.$cyc.cwl.waterlevel.nc" + +rwps_waterlevel="$frc/$meshname.$PDY.$cyc.waterlevel.nc" + +varnames="zeta" + +## STOFS interpolation + +if [ ! -f "$stofs_wghts" ]; then + echo "missing stofs interpolation weights file: $stofs_wghts" + echo "compute with script compute_unstr_to_rwps_interp_weights.sh" + exit 1 +fi +if [ ! -f "$stofs_dists" ]; then + echo "missing stofs distance to boundary file: $stofs_dists" + echo "compute with script compute_unstr_to_rwps_interp_weights.sh" + exit 1 +fi + +python interpolate_with_weights.py $stofslev $stofs_wghts $rwps_waterlevel $varnames 0 +python add_mesh_geom_to_file.py $rwps_waterlevel $mesh +python add_err_var_to_file.py $rwps_waterlevel $stofs_dists 1. diff --git a/ush/preprocess/process_wind.sh b/ush/preprocess/process_wind.sh new file mode 100644 index 0000000..525cd57 --- /dev/null +++ b/ush/preprocess/process_wind.sh @@ -0,0 +1,122 @@ +#!/bin/bash + +cd $HOMErwps/ush/preprocess + +winddir="$tmp/wind.$PDY.$cyc" +windvars="UGRD_10maboveground:VGRD_10maboveground" + +# extract mesh name from file path +meshname="${mesh##*/}" +# remove .msh suffix from mesh name +meshname="${meshname: 0: -4}" + +# incorporate meshname date and cycle into output directory name to avoid +# applying winds to wrong mesh + +outdir="$tmp/rwps_wind.$meshname.$PDY.$cyc" + +rwps_wind="$frc/$meshname.$PDY.$cyc.wind.nc" + +echo "outputing files to: $outdir" + +nbm_oc="$winddir/nbm.$PDY.$cyc.wind10m.oc.nc" +nbm_oc_uv="$winddir/nbm.$PDY.$cyc.wind10m.oc.uv.nc" + +rrfs_pr="$winddir/rrfs.$PDY.$cyc.wind10m.pr.nc" +rrfs_hi="$winddir/rrfs.$PDY.$cyc.wind10m.hi.nc" +rrfs_na="$winddir/rrfs.$PDY.$cyc.wind10m.na.nc" +rrfs_ak="$winddir/rrfs.$PDY.$cyc.wind10m.ak.nc" +rrfs_conus="$winddir/rrfs.$PDY.$cyc.wind10m.conus.nc" + +rwps_oc="$outdir/$meshname.$PDY.$cyc.wind10m.nbm.oc.nc" +rwps_oc_ti="$outdir/$meshname.$PDY.$cyc.wind10m.nbm.oc.ti.nc" +rwps_pr="$outdir/$meshname.$PDY.$cyc.wind10m.rrfs.pr.nc" +rwps_hi="$outdir/$meshname.$PDY.$cyc.wind10m.rrfs.hi.nc" +rwps_na="$outdir/$meshname.$PDY.$cyc.wind10m.rrfs.na.nc" +rwps_ak="$outdir/$meshname.$PDY.$cyc.wind10m.rrfs.ak.nc" +rwps_conus="$outdir/$meshname.$PDY.$cyc.wind10m.rrfs.conus.nc" + + +mkdir -p $outdir + +##LocalFS = [ rwps_pr, rwps_hi, rwps_ak, rwps_conus, rwps_na] # file names +##VarFS = [ 4. , 4. , 9. , 16. , 25. ] # (m m /s /s) +##LambdaFS = [ 150. , 200. , 500. , 1000. , 1500. ] # (km) +nbm_oc_wghts="$fix/InterpolationWeights.$meshname.nbm.oc.nc" +rrfs_hi_wghts="$fix/InterpolationWeights.$meshname.rrfs.hi.nc" +rrfs_pr_wghts="$fix/InterpolationWeights.$meshname.rrfs.pr.nc" +rrfs_ak_wghts="$fix/InterpolationWeights.$meshname.rrfs.ak.nc" +rrfs_na_wghts="$fix/InterpolationWeights.$meshname.rrfs.na.nc" +rrfs_conus_wghts="$fix/InterpolationWeights.$meshname.rrfs.conus.nc" + +nbm_oc_dist="$fix/DistToBndy.$meshname.nbm.oc.nc" +rrfs_hi_dist="$fix/DistToBndy.$meshname.rrfs.hi.nc" +rrfs_pr_dist="$fix/DistToBndy.$meshname.rrfs.pr.nc" +rrfs_ak_dist="$fix/DistToBndy.$meshname.rrfs.ak.nc" +rrfs_na_dist="$fix/DistToBndy.$meshname.rrfs.na.nc" +rrfs_conus_dist="$fix/DistToBndy.$meshname.rrfs.conus.nc" + + +( + #Convert NBM speed and direction to u,v + python spd_dir_to_uv_nbm.py $nbm_oc $nbm_oc_uv + # If the interpolation weights do not already exist for the domains create them + # also creates distance to boundary used in prescribed error covariance specification + [ ! -f "$nbm_oc_wghts" ] && python compute_gridded_to_rwps_interp_weights.py $nbm_oc_uv $mesh 1 + python interpolate_with_weights.py $nbm_oc_uv $nbm_oc_wghts $rwps_oc $windvars 3 > $tmp/intrp.rrfs.oc.out + python add_mesh_geom_to_file.py $rwps_oc $mesh +)& + +( + [ ! -f "$rrfs_hi_wghts" ] && python compute_gridded_to_rwps_interp_weights.py $rrfs_hi $mesh $rrfs_hi_wghts $rrfs_hi_dist + python interpolate_with_weights.py $rrfs_hi $rrfs_hi_wghts $rwps_hi $windvars -1 > $tmp/intrp.rrfs.hi.out + #add mesh geometry into interpolated file + python add_mesh_geom_to_file.py $rwps_hi $mesh + # Add error covariance field to files with interpolated fields for bayesian update + # Based on distance to boundary of input field and commant line parameters InternalVariance:BoundaryVariance:LengthScale(km) + python add_err_var_to_file.py $rwps_hi $rrfs_hi_dist 4.:40.:200. +)& + +( + [ ! -f "$rrfs_pr_wghts" ] && python compute_gridded_to_rwps_interp_weights.py $rrfs_pr $mesh $rrfs_pr_wghts $rrfs_pr_dist + python interpolate_with_weights.py $rrfs_pr $rrfs_pr_wghts $rwps_pr $windvars -1 > $tmp/intrp.rrfs.pr.out + python add_mesh_geom_to_file.py $rwps_pr $mesh + python add_err_var_to_file.py $rwps_pr $rrfs_pr_dist 4.:40.:150. +)& + +( + [ ! -f "$rrfs_ak_wghts" ] && python compute_gridded_to_rwps_interp_weights.py $rrfs_ak $mesh $rrfs_ak_wghts $rrfs_ak_dist + python interpolate_with_weights.py $rrfs_ak $rrfs_ak_wghts $rwps_ak $windvars -1 > $tmp/intrp.rrfs.ak.out + python add_mesh_geom_to_file.py $rwps_ak $mesh + python add_err_var_to_file.py $rwps_ak $rrfs_ak_dist 9.:90.:500. +)& + +( + [ ! -f "$rrfs_conus_wghts" ] && python compute_gridded_to_rwps_interp_weights.py $rrfs_conus $mesh $rrfs_conus_wghts $rrfs_conus_dist + python interpolate_with_weights.py $rrfs_conus $rrfs_conus_wghts $rwps_conus $windvars -1 > $tmp/intrp.rrfs.conus.out + python add_mesh_geom_to_file.py $rwps_conus $mesh + python add_err_var_to_file.py $rwps_conus $rrfs_conus_dist 16.:160.:1000. +)& + +( + [ ! -f "$rrfs_na_wghts" ] && python compute_gridded_to_rwps_interp_weights.py $rrfs_na $mesh $rrfs_na_wghts $rrfs_na_dist + python interpolate_with_weights.py $rrfs_na $rrfs_na_wghts $rwps_na $windvars -1 > $tmp/intrp.rrfs.na.out + python add_mesh_geom_to_file.py $rwps_na $mesh + python add_err_var_to_file.py $rwps_na $rrfs_na_dist 50.:500.:1500. +)& + +wait; + +#Interpolate NBM in time to times within the NBM forecast covered by the RRFS forecast +python interp_time.py $rwps_oc $rwps_pr $rwps_oc_ti $windvars +#add mesh geometry into file +python add_mesh_geom_to_file.py $rwps_oc_ti $mesh +#add prescribed error covariance for nbm oc domain (assumed constant 100. (m/s)^2 ) +python add_err_var_to_file.py $rwps_oc_ti $nbm_oc_dist 100. + +cp $rwps_oc_ti $rwps_wind +[ ! -f "$rwps_hi" ] && python bayes_forecast_update.py $rwps_wind $rwps_hi $rwps_wind $windvars +[ ! -f "$rwps_pr" ] && python bayes_forecast_update.py $rwps_wind $rwps_pr $rwps_wind $windvars +[ ! -f "$rwps_ak" ] && python bayes_forecast_update.py $rwps_wind $rwps_ak $rwps_wind $windvars +[ ! -f "$rwps_conus" ] && python bayes_forecast_update.py $rwps_wind $rwps_conus $rwps_wind $windvars +[ ! -f "$rwps_na" ] && python bayes_forecast_update.py $rwps_wind $rwps_na $rwps_wind $windvars diff --git a/ush/preprocess/rrfs/make_rrfs_wind.sh b/ush/preprocess/rrfs/make_rrfs_wind.sh new file mode 100755 index 0000000..337324c --- /dev/null +++ b/ush/preprocess/rrfs/make_rrfs_wind.sh @@ -0,0 +1,37 @@ +#!/bin/bash + +# This script takes rrfs grib2 forecast files, extracts 10m u and v wind +# components and outputs to netcdf. Comand line arguments are +# forecast time: YYYYMMDD +# forecast cycle: CC +# frecast region: hi,pr,na,ak or conus +# +# For example call as: +# sh make_rrfs_wind.sh 20260428 00 pr +# to produce output file: +# rrfs.20260428.00.wind10m.pr.nc +# +echo "make_rrfs_wind.sh fetching rrfs : time = $1, cycle = $2, domain = $3" + +INPUT_DIR="$COMINrrfs" +OUTPUT_DIR="$COMINlocal/wind.$1.$2" +OUTPUT_FILE="$OUTPUT_DIR/rrfs.$1.$2.wind10m.$3.nc" + +mkdir -p "$OUTPUT_DIR" +# Remove existing output file to avoid mixing old data +rm -f "$OUTPUT_FILE" + +echo "writing 10m wind from $INPUT_DIR to $OUTPUT_FILE" + +# Note: Ensure files are in the correct chronological order (e.g., sorted by name) +for file in $(ls "$INPUT_DIR"/rrfs.t$2z.2dfld.*km.f*.$3.grib2 | sort); do + echo "Processing: $file" + # Convert and append to the NetCDF file + # -netcdf: specifies the output format and filename + # -append: adds data to the existing netcdf file instead of overwriting + wgrib2 "$file" -match ":UGRD:10 m" -append -netcdf "$OUTPUT_FILE" + wgrib2 "$file" -match ":VGRD:10 m" -append -netcdf "$OUTPUT_FILE" +done + +echo "rrfs processing complete for forecast date $1, cycle $2, domain $3" +echo "output written to: $OUTPUT_FILE" diff --git a/ush/preprocess/rtofs/get_rtofs.sh b/ush/preprocess/rtofs/get_rtofs.sh new file mode 100755 index 0000000..b68930c --- /dev/null +++ b/ush/preprocess/rtofs/get_rtofs.sh @@ -0,0 +1,15 @@ +#!/bin/bash + +#Retrieve global RTOFS currents and consolidate into a single NetCDF file + +cd $HOMErwps/ush/preprocess + +tmpdir="$COMINlocal/tmp.rtofs.$PDY" +filesin="$COMINrtofs/*prog.nc" +flout="$COMINlocal/rtofs.$PDY.nc" + +mkdir -p $tmpdir +cp $filesin $tmpdir/ +python rtofs/get_rtofs_fcst.py $tmpdir $flout + +#rm -rf $tmpdir diff --git a/ush/preprocess/rtofs/get_rtofs_fcst.py b/ush/preprocess/rtofs/get_rtofs_fcst.py new file mode 100644 index 0000000..480fc68 --- /dev/null +++ b/ush/preprocess/rtofs/get_rtofs_fcst.py @@ -0,0 +1,95 @@ +import numpy as np +import os +import netCDF4 as nc +import sys + +sys.path.append(os.path.abspath(os.path.join(os.path.dirname(__file__), '..'))) +import interp_utilities as iutil + + + +# Consolidate u,v surface currents from global RTOFS into a single NetCDF file. + +dirin=sys.argv[1] +flout=sys.argv[2] + +ncfiles = os.listdir(dirin) +nt=len(ncfiles) +k=0 +flin=dirin+"/"+ncfiles[k] +data=nc.Dataset(flin,"r") +X=data["Longitude"][:,:] +Y=data["Latitude"][:,:] +fill_value0=data["u_velocity"]._FillValue +nx=X.shape[0] +ny=X.shape[1] +U=np.zeros((nt,1,nx,ny)) +V=np.zeros((nt,1,nx,ny)) +MT=np.zeros(nt) +time=np.zeros(nt) +print(U.shape) +for k in range(nt): + flin=dirin+"/"+ncfiles[k] + data=nc.Dataset(flin,"r") + U0=data["u_velocity"][:] + V0=data["v_velocity"][:] + T0=data["MT"][:] + print(k) + U[k,0,:,:]=U0[0,0,:,:] + V[k,0,:,:]=V0[0,0,:,:] + MT[k]=T0 + tmp=iutil.ConvertTimeToUnixTime(flin,"MT") + time[k]=tmp[0] + +tindx = np.argsort(time) +time=time[tindx] +MT=MT[tindx] +U=U[tindx,:,:,:] +V=V[tindx,:,:,:] + +Layers=data["Layer"][:] + +with nc.Dataset(flout, 'w', format='NETCDF4') as ncout: + + ncout.createDimension('Layer' , 1) + ncout.createDimension('X' , nx) + ncout.createDimension('Y' , ny) + ncout.createDimension('time', nt) + + lon_var=ncout.createVariable('longitude', 'f4', ('X','Y')) + var=data["Longitude"] + iutil.CopyAttributes(var, lon_var) + lon_var[:,:]=X[:,:] + + lat_var=ncout.createVariable('latitude', 'f4', ('X','Y')) + var=data["Latitude"] + iutil.CopyAttributes(var, lat_var) + lat_var[:,:]=Y[:,:] + + # Keep redundent time variables (time and MT) for conviniences + time_var=ncout.createVariable('time', 'f8', ('time',)) + time_var.units = 'seconds since 1970-01-01 00:00:00.0 0:00' + time_var.long_name = 'verification time generated by wgrib2 function verftime()' + time_var.standard_name = 'time' + time_var.axis = 'T' + time_var[:] = time[:] + + MT_var=ncout.createVariable("MT", 'f8', ('time',)) + var=data["MT"] + iutil.CopyAttributes(var, MT_var) + MT_var[:] = MT[:] + + u_var = ncout.createVariable('u-vel', 'f4', ('time','Layer','X','Y'),fill_value = fill_value0) + var=data["u_velocity"] + iutil.CopyAttributes(var, u_var) + u_var[:,:,:,:] = U + + v_var = ncout.createVariable('v-vel', 'f4', ('time','Layer','X','Y'),fill_value = fill_value0) + var=data["v_velocity"] + iutil.CopyAttributes(var, v_var) + v_var[:,:,:,:] = V + + l_var = ncout.createVariable('Layer', 'i4', ('Layer',)) + var=data["Layer"] + iutil.CopyAttributes(var, l_var) + l_var[:] = Layers[0] diff --git a/ush/preprocess/rtofs/get_rtofs_ice.sh b/ush/preprocess/rtofs/get_rtofs_ice.sh new file mode 100755 index 0000000..4621f30 --- /dev/null +++ b/ush/preprocess/rtofs/get_rtofs_ice.sh @@ -0,0 +1,18 @@ +#!/bin/bash + +#Retrieve global RTOFS ice forecast and consolidate in single NetCDF filesin + +tmpdir="$COMINlocal/tmp.rtofsIce.$PDY" +filesin="$COMINrtofs/*ice.nc" +dirout="$COMINlocal/ice.$PDY.$cyc" +flout="$dirout/rtofs.ice.$PDY.nc" + +mkdir $tmpdir +cp $filesin $tmpdir/ +echo $tmpdir +echo $flout + +PDYCC="${PDY}${cyc}" +echo $PDYCC +python rtofs/get_rtofs_ice_fcst.py $tmpdir $PDYCC $flout + diff --git a/ush/preprocess/rtofs/get_rtofs_ice_fcst.py b/ush/preprocess/rtofs/get_rtofs_ice_fcst.py new file mode 100644 index 0000000..8b3e769 --- /dev/null +++ b/ush/preprocess/rtofs/get_rtofs_ice_fcst.py @@ -0,0 +1,127 @@ +import numpy as np +import os +import netCDF4 as nc +import sys + +import datetime +import re + +def get_rtofs_unix_time(file_path: str) -> int: + """ + Extracts the valid UTC Unix timestamp from a NOAA RTOFS file path. + Example path: '/pub/data/nccf/com/rtofs/prod/rtofs.20260813/rtofs_glo.t12z.f024.archv.nc' + """ + # Regex to find 8-digit date, cycle hour (tXXz), and forecast hour (fXXX) + match = re.search(r'rtofs\.(\d{8}).*?\.t(\d{2})z\.f(\d{3})', file_path) + + if not match: + raise ValueError("Could not parse RTOFS naming convention from path.") + + date_str, cycle_str, forecast_str = match.groups() + + # 1. Parse base model initialization time to a UTC datetime object + base_time_str = f"{date_str} {cycle_str}" + base_dt = datetime.datetime.strptime(base_time_str, "%Y%m%d %H").replace(tzinfo=datetime.timezone.utc) + + # 2. Add the forecast hour offset + forecast_hours = int(forecast_str) + valid_dt = base_dt + datetime.timedelta(hours=forecast_hours) + + # 3. Convert valid UTC datetime to epoch integer + return int(valid_dt.timestamp()) + +# --- Example Usage --- +sample_path = "/pub/data/nccf/com/rtofs/prod/rtofs.20260813/rtofs_glo.t12z.f024.archv.nc" +unix_time = get_rtofs_unix_time(sample_path) + +print(f"Unix Timestamp: {unix_time}") +# Output: 1786708800 (Reflects: 2026-08-14 12:00:00 UTC) + +#Consolidate RTOFS ice forecast files to a single NetCDF file. + +# Get the path relative to this file and add to sys.path +sys.path.append(os.path.abspath(os.path.join(os.path.dirname(__file__), '..'))) +import interp_utilities as iutil + +dirin=sys.argv[1] +PDYCC=sys.argv[2] +flout=sys.argv[3] +print(PDYCC) +PDY=PDYCC[0:8] +CYC=PDYCC[8:10] +print(PDY) +print(CYC) + + +ncfiles = os.listdir(dirin) +print(ncfiles) +nt=len(ncfiles) + +nt=17 + + +k=0 +flin=dirin+"/"+ncfiles[k] +data=nc.Dataset(flin,"r") +X=data["Longitude"][:,:] +Y=data["Latitude"][:,:] +fill_value0=data["ice_coverage"]._FillValue +nx=X.shape[0] +ny=X.shape[1] +F=np.zeros((nt,nx,ny),dtype=np.float32) +MT=np.zeros(nt) +time=np.zeros(nt) +print("nx = "+str(nx)) +print("ny = "+str(ny)) +print(F.shape) +for k in range(nt): + flin=dirin+"/"+ncfiles[k] + data=nc.Dataset(flin,"r") + F0=data["ice_coverage"][:] + print(F0.shape) + print(F.shape) + print(str(k)+" of "+ str(nt)) + F[k,:,:]=F0[0,:,:] + tmp=iutil.FileNameToUnixTime(flin,PDY,CYC) + time[k]=tmp + print(flin) + print(time[k]) + +tindx = np.argsort(time) +time=time[tindx] +MT=MT[tindx] +F=F[tindx,:,:] + +with nc.Dataset(flout, 'w', format='NETCDF4') as ncout: + + ncout.createDimension('X' , nx) + ncout.createDimension('Y' , ny) + ncout.createDimension('time', nt) + + lon_var=ncout.createVariable('longitude', 'f4', ('X','Y')) + var=data["Longitude"] + iutil.CopyAttributes(var, lon_var) + lon_var[:,:]=X[:,:] + + lat_var=ncout.createVariable('latitude', 'f4', ('X','Y')) + var=data["Latitude"] + iutil.CopyAttributes(var, lat_var) + lat_var[:,:]=Y[:,:] + + # Keep redundent time variables (time and MT) for conviniences + time_var=ncout.createVariable('time', 'f8', ('time',)) + time_var.units = 'seconds since 1970-01-01 00:00:00.0 0:00' + time_var.long_name = 'verification time generated by wgrib2 function verftime()' + time_var.standard_name = 'time' + time_var.axis = 'T' + time_var[:] = time[:] + + MT_var=ncout.createVariable("MT", 'f8', ('time',)) + var=data["MT"] + iutil.CopyAttributes(var, MT_var) + MT_var[:] = MT[:] + + f_var = ncout.createVariable('ICEC_surface', 'f4', ('time','X','Y'),fill_value = fill_value0) + var=data["ice_coverage"] + iutil.CopyAttributes(var, f_var) + f_var[:,:,:] = F[:,:,:] diff --git a/ush/preprocess/spd_dir_to_uv_nbm.py b/ush/preprocess/spd_dir_to_uv_nbm.py new file mode 100644 index 0000000..91d46c5 --- /dev/null +++ b/ush/preprocess/spd_dir_to_uv_nbm.py @@ -0,0 +1,79 @@ +import os +import argparse +import numpy as np +import netCDF4 as nc +import sys +import math +import interp_utilities as iutil +# Converts NBM direction and magnitude wind variables to U,V +# Command line inputs: +# argumnent 1 is the netcdf input file name with variables: WIND_10maboveground and WDIR_10maboveground +# argumnent 2 is the netcdf output file name which will have vector wind variables: UGRD_10maboveground and VGRD_10maboveground +# call as: +# python spd_dir_to_uv_nbm.py nbm.20260502.00.wind10m.oc.nc nbm.20260502.00.wind10m.oc.uv.nc + +flin=sys.argv[1] +flout=sys.argv[2] + +fcst0 = nc.Dataset(flin,"r") + +x=np.asarray(fcst0["longitude"][:]) +y=np.asarray(fcst0["latitude"][:]) +t=np.asarray(fcst0["time"][:]) +spd=np.asarray(fcst0["WIND_10maboveground"][:,:,:]) +theta=np.asarray(fcst0["WDIR_10maboveground"][:,:,:]) + +spdV=fcst0["WIND_10maboveground"] +fill_value0 = spdV._FillValue + +u=-1*spd*np.sin(theta*np.pi/180.) +v=-1*spd*np.cos(theta*np.pi/180.) +u[np.where(spd==fill_value0)]=fill_value0 +v[np.where(spd==fill_value0)]=fill_value0 + +nx=len(x) +ny=len(y) +nt=len(t) + +with nc.Dataset(flout, 'w', format='NETCDF4') as ncout: + # Create dimensions + ncout.createDimension('lon' , nx) # Unlimited dimension + ncout.createDimension('lat' , ny) + ncout.createDimension('time', nt) + + lon_var=ncout.createVariable('longitude', 'f8', ('lon',)) + lon_var.units = 'degree_east' + lon_var.long_name = 'longitude' + lon_var.standard_name = 'longitude' + lon_var.axis = 'lon' + lon_var[:]=x[:] + + lat_var=ncout.createVariable('latitude', 'f8', ('lat',)) + lat_var.units = 'degree_north' + lat_var.long_name = 'latitude' + lat_var.standard_name = 'latitude' + lat_var.axis = 'lat' + lat_var[:]=y[:] + + time_var=ncout.createVariable('time', 'f8', ('time',)) + time_varin=fcst0["time"] + iutil.CopyAttributes(time_varin, time_var) + time_var[:]=t[:] + + u_var=ncout.createVariable('UGRD_10maboveground', 'f4', ('time','lat','lon'),fill_value = fill_value0) + u_var.long_name = 'U-Component of Wind' + u_var.units = 'm/s' + u_var.standard_name = 'UGRD_10maboveground' + u_var.level = '10 m above ground' + u_var[:,:,:]=u[:,:,:] + + v_var=ncout.createVariable('VGRD_10maboveground', 'f4', ('time','lat','lon'),fill_value = fill_value0) + v_var.long_name = 'V-Component of Wind' + v_var.units = 'm/s' + v_var.standard_name = 'VGRD_10maboveground' + v_var.level = '10 m above ground' + v_var[:,:,:]=v[:,:,:] + + ncout.close + + diff --git a/ush/preprocess/stofs/get_stofs.sh b/ush/preprocess/stofs/get_stofs.sh new file mode 100755 index 0000000..7985c6e --- /dev/null +++ b/ush/preprocess/stofs/get_stofs.sh @@ -0,0 +1,26 @@ +#!/bin/bash + +#script for retrieving stofs current and/or water level + +if [ "$#" -lt 3 ]; then + echo "No arguments field argument." + echo "Retrieving both stofs currents and stofs water level." + fields="current, waterlevel" +else + fields=$3 +fi + +outdir=$tmp/stofs.$PDY.$cyc + +mkdir -p $outdir + +if [[ "$fields" == *"current"* ]]; then + echo retrieving stofs current for $PDY cycle $cyc + cp $COMINstofs/stofs_2d_glo.t"$cyc"z.fields.cwl.vel.nc $outdir/ +fi + +if [[ "$fields" == *"level"* ]]; then + echo retrieving stofs water level for $PDY cycle $cyc + cp $COMINstofs/stofs_2d_glo.t"$cyc"z.fields.cwl.nc $outdir/ +fi + diff --git a/versions/preproc.wcoss2.ver b/versions/preproc.wcoss2.ver new file mode 100755 index 0000000..a53605b --- /dev/null +++ b/versions/preproc.wcoss2.ver @@ -0,0 +1,12 @@ +#! /usr/bin/env bash +export PrgEnv_intel_ver=8.5.0 +export interl_ver=19.1.3.304 +export craype_ver=2.7.17 +export cray_mpich_ver=8.1.19 +export hdf5_C_ver=1.14.0 +export netcdf_C_ver=4.9.2 +export esmf_C_ver=8.6.0 +export ve_hafs_ver=2.1 + +export intel_oneapi_ver=2022.2.0.262 +export wgrib2_ver=2.0.8