Skip to content

Commit

Permalink
adding 2ds ice for rtofs
Browse files Browse the repository at this point in the history
  • Loading branch information
@rgrumbine
rgrumbine committed Oct 22, 2024
1 parent 4e62fcd commit f74e3ac
Show file tree
Hide file tree
Showing 4 changed files with 64 additions and 76 deletions.
0 nwp/new_cafs.py → nwp/cafs/new_cafs.py
View file
File renamed without changes.
24 changes: 11 additions & 13 deletions nwp/cafs/runup.sh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,25 +14,23 @@
# Path to cafs output:
export base=/home/Robert.Grumbine/clim_data/cafs/

# location of python and support
cd /home/Robert.Grumbine/rgdev/toolbox/nwp

# This must be more or less exactly this:
source /home/Robert.Grumbine/rgref/env3.12/bin/activate
export PYTHONPATH=$PYTHONPATH:/home/Robert.Grumbine/rgdev/toolbox/nwp
# location of python and support
cd /home/Robert.Grumbine/rgdev/toolbox/nwp/cafs

#These can be anything of convenience
export MPLCONFIGDIR=$HOME/rgexpt/
export OUTDIR=$HOME/rgdev/cafs_nwp

#------------------------------------------------------
#tag=20220401
#tag=20241015
tag=20240711
#debug: tag=20241001
tag=`date +"%Y%m%d"`
#debug:
tag=20241022

#reverse -- now to past
end=`date +"%Y%m%d"`
end=20220401
end=20230401
#debug: end=20241011

while [ $tag -ge $end ]
Expand All @@ -41,14 +39,14 @@ do
yy=`echo $tag | cut -c1-4`
mm=`echo $tag | cut -c5-6`
dd=`echo $tag | cut -c7-8`
if [ ! -f $HOME/rgdev/cafs_nwp/out.$tag ] ; then
time python3 new_cafs.py $yy $mm $dd > $OUTDIR/out.$tag
if [ ! -f $HOME/rgdev/cafs_nwp/${yy}/out.$tag ] ; then
time python3 new_cafs.py $yy $mm $dd > $OUTDIR/${yy}/out.$tag
fi
if [ -f nwp_${tag}_6.png ] ; then
mv nwp_${tag}_*.png $OUTDIR
mv nwp_${tag}_*.png $OUTDIR/${yy}
fi
if [ -f path_${tag}_6.kml ] ; then
mv path_${tag}_*.kml $OUTDIR
mv path_${tag}_*.kml $OUTDIR/${yy}
fi

tag=`expr $tag - 1`
Expand Down
4 changes: 0 additions & 4 deletions nwp/functions.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,10 +6,6 @@
import numpy as np
import numpy.ma as ma

#debug: import netCDF4
#debug: from netCDF4 import Dataset
#debug: import networkx as netx

#--------------------------------------------------------
def find(lons, lats, lonin, latin):
#debug print("lon, lat in:",lonin, latin, flush=True)
Expand Down
112 changes: 53 additions & 59 deletions nwp/rtofs/a.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,23 +10,25 @@
#--------------------------------------------------------------
#fin = Dataset("../dcom/rtofs_glo_2ds_f000_prog.nc","r")
#fin = Dataset("../dcom/rtofs_glo_2ds_f000_diag.nc","r")
#fin = Dataset("../dcom/rtofs_glo_2ds_f000_ice.nc","r")
fin = Dataset("rtofs_glo.t00z.n00.cice_inst.nc","r")

nx = len(fin.dimensions["ni"])
ny = len(fin.dimensions["nj"])

print(nx, ny)
#debug: exit(0)

#extract longitude, latitude, aice
# -- node properties
lons = fin.variables["TLON"][:,:]
lats = fin.variables["TLAT"][:,:]

#sst = fin.variables["sst"][0,:,:]
#sst = fin.variables["ice_thickness"][0,:,:]
aice = fin.variables["aice"][0,:,:]
# cice_inst
#fin = Dataset("rtofs_glo.t00z.n00.cice_inst.nc","r")
#nx = len(fin.dimensions["ni"])
#ny = len(fin.dimensions["nj"])
#lons = fin.variables["TLON"][:,:]
#lats = fin.variables["TLAT"][:,:]
#sst = fin.variables["sst"][0,:,:]
#aice = fin.variables["aice"][0,:,:]

# 2ds_ice
fin = Dataset("rtofs_glo_2ds_f000_ice.nc","r")
nx = len(fin.dimensions['X'])
ny = len(fin.dimensions['Y'])
lons = fin.variables["Longitude"][:,:]
lats = fin.variables["Latitude"][:,:]
#not available: sst = fin.variables["sst"][0,:,:]
aice = fin.variables["ice_coverage"][0,:,:]
#debug: check distance of ice-free arctic aice.fill(0.)

#debug:
print("lons: ",lons.max(), lons.min() )
Expand Down Expand Up @@ -78,36 +80,9 @@
#debug: exit(0)

#----------------------------------------------------------------
#print("nx,ny,nx*ny:",nx,ny,nx*ny)
def find(lons, lats, lonin, latin):
#debug print("lon, lat in:",lonin, latin, flush=True)
tmpx = lons - lonin
tmpy = lats - latin
#debug print("x ",tmpx.max(), tmpx.min(), lons.max(), lons.min(), flush=True )

xmask = ma.masked_outside(tmpx, -0.5, 0.5)
xin = xmask.nonzero()
wmask = ma.logical_and(xmask, ma.masked_outside(tmpy, +0.5, -0.5) )
win = wmask.nonzero()

imin = -1
jmin = -1
dxmin = 999.
dymin = 999.
dmin = 999.
for k in range(0, len(win[0]) ):
i = win[1][k]
j = win[0][k]
#debug print(k,i,j,abs(tmpx[j,i]), abs(tmpy[j,i]), dxmin, dymin, dmin, flush=True)
#if (abs(tmpx[j,i]) < dxmin and abs(tmpy[j,i]) < dymin):
if (sqrt(tmpx[j,i]**2 + tmpy[j,i]**2) < dmin):
imin = i
jmin = j
dxmin = abs(tmpx[j,i])
dymin = abs(tmpy[j,i])
dmin = sqrt(tmpx[j,i]**2 + tmpy[j,i]**2)
#print("dmin:",imin, jmin, dmin, dxmin, dymin)
return (imin,jmin)
from functions import *

#----------------------------------------------------------------

#tlat = 74.0
#for iii in range (0, 400):
Expand All @@ -131,18 +106,19 @@ def find(lons, lats, lonin, latin):
#--------------------------------------------------------------

# Construct nodes -- limit area to keep run time manageable:
latmin = 65.0
latmax = 88.0
latmin = 64.0
latmax = 82.0
#lonmin = 185.0-360.
#lonmax = 290.0-360.
lonmin = -175.0
lonmax = -70.0
xmask = ma.masked_outside(lons, lonmin, lonmax)
xin = xmask.nonzero()
#print(len(xin), len(xin[0]))
#debug: print(len(xin), len(xin[0]), flush=True)
xmask = ma.logical_and(xmask, ma.masked_outside(lats, latmin, latmax))
xin = xmask.nonzero()
#print(len(xin), len(xin[0]))
#debug:
print("number of points:", len(xin[0]), flush=True)

xmask = ma.logical_and(xmask, aice < 1000.)
xin = xmask.nonzero()
Expand Down Expand Up @@ -173,6 +149,14 @@ def find(lons, lats, lonin, latin):
# RG: tripolar grid means adjacent geographic points aren't always i,j adjacent
# fix!
# Construct edges between nodes:

#cost_type =
#1 -> steps in i,j space
#2 -> meters
#3 -> weighted by polar class
#4 -> weight by 1./(1.1-aice)
cost_type = 4

for k in range(0, len(xin[0])):
i = xin[1][k]
j = xin[0][k]
Expand All @@ -186,20 +170,27 @@ def find(lons, lats, lonin, latin):

if (im >= 0):
if (nodemap[j,im] != -1):
G.add_edge(n, nodemap[j,im], weight = 1.)
weight = cost(cost_type, lat1 = lats[j,i], lon1 = lons[j,i],
lat2 = lats[j,im], lon2 = lons[j,im], aice = aice[j,i])
G.add_edge(n, nodemap[j,im], weight = weight)
if (ip < nx):
if (nodemap[j,ip] != -1):
G.add_edge(n, nodemap[j,ip], weight = 1.)
weight = cost(cost_type, lat1 = lats[j,i], lon1 = lons[j,i],
lat2 = lats[j,ip], lon2 = lons[j,ip], aice = aice[j,i])
G.add_edge(n, nodemap[j,ip], weight = weight)

if (jp < ny ):
if (nodemap[jp,i] != -1):
G.add_edge(n, nodemap[jp,i], weight = 1.)
weight = cost(cost_type, lat1 = lats[j,i], lon1 = lons[j,i], lat2 = lats[jp,i], lon2 = lons[jp,i], aice = aice[j,i])
G.add_edge(n, nodemap[jp,i], weight = weight)
if (im >= 0):
if (nodemap[jp,im] != -1):
G.add_edge(n, nodemap[jp,im], weight = 1.)
weight = cost(cost_type, lat1 = lats[j,i], lon1 = lons[j,i], lat2 = lats[jp,im], lon2 = lons[jp,im], aice = aice[j,i])
G.add_edge(n, nodemap[jp,im], weight = weight)
if (ip < nx):
if (nodemap[jp,ip] != -1):
G.add_edge(n, nodemap[jp,ip], weight = 1.)
weight = cost(cost_type, lat1 = lats[j,i], lon1 = lons[j,i], lat2 = lats[jp,ip], lon2 = lons[jp,ip], aice = aice[j,i])
G.add_edge(n, nodemap[jp,ip], weight = weight)
#RG: a guess about the archipelago seam
else:
tmpi = i
Expand All @@ -210,13 +201,16 @@ def find(lons, lats, lonin, latin):

if (jm >= 0 ):
if (nodemap[jm,i] != -1):
G.add_edge(n, nodemap[jm,i], weight = 1.)
weight = cost(cost_type, lat1 = lats[j,i], lon1 = lons[j,i], lat2 = lats[jm,i], lon2 = lons[jm,i], aice = aice[j,i])
G.add_edge(n, nodemap[jm,i], weight = weight)
if (im >= 0):
if (nodemap[jm,im] != -1):
G.add_edge(n, nodemap[jm,im], weight = 1.)
weight = cost(cost_type, lat1 = lats[j,i], lon1 = lons[j,i], lat2 = lats[jm,im], lon2 = lons[jm,im], aice = aice[j,i])
G.add_edge(n, nodemap[jm,im], weight = weight)
if (ip < nx):
if (nodemap[jm,ip] != -1):
G.add_edge(n, nodemap[jm,ip], weight = 1.)
weight = cost(cost_type, lat1 = lats[j,i], lon1 = lons[j,i], lat2 = lats[jm,ip], lon2 = lons[jm,ip], aice = aice[j,i])
G.add_edge(n, nodemap[jm,ip], weight = weight)

#debug:
print("Have constructed graph, number of nodes, edges =",k, len(G.edges), flush=True)
Expand Down Expand Up @@ -249,7 +243,7 @@ def find(lons, lats, lonin, latin):

pseudo_length = netx.dijkstra_path_length(G,start, finish)

print("dijkstra length ", len(path), flush=True)
print("dijkstra length ", len(path), pseudo_length, flush=True)

tlons = np.zeros((len(path)))
tlats = np.zeros((len(path)))
Expand Down

0 comments on commit f74e3ac

Please sign in to comment.