Commit 1e4888cf authored by Carl Schreck's avatar Carl Schreck
Browse files

Overdue catchup

parent 1b16241f
import xarray as xr
import numpy as np
level = 850
var_name = 'uwnd'
months = np.arange(6,11)
years = np.arange(1905, 2015)
lat_slice = slice(0, 60)
lon_slice = slice(180, 360)
if level is not None:
in_path = (
f'http://apdrc.soest.hawaii.edu:80/dods/public_data/Reanalysis_Data/'
f'NOAA_20th_Century/V3/daily/pressure/{var_name}'
)
else:
in_path = (
f'http://apdrc.soest.hawaii.edu:80/dods/public_data/Reanalysis_Data/'
f'NOAA_20th_Century/V3/daily/monolevel/{var_name}'
)
out_path = (
f'/store/carl/people/klotzbach/mjo_landfall/'
f'{var_name}{level}.daily.total.nc'
)
ds = xr.open_dataset(in_path)
print(ds)
is_in_months = ds['time.month'].isin(months)
is_in_years = ds['time.year'].isin(years)
ds = ds.sel(time=(is_in_years & is_in_months), lat=lat_slice, lon=lon_slice,)
print(ds)
ds = ds.sel(lev=level, drop=True)
print(ds)
ds = ds.rename({var_name: f'{var_name}{level}'})
print(ds)
ds.to_netcdf(out_path)
#!/bin/bash --login
#if [ -z "$1" ]; then
# QUEUE=`pick_queue`
#else
# QUEUE=$1
#fi
QUEUE=allnodes
echo $QUEUE
NCL_SCRIPT=calculate_dailies
NCL_DIR=`pwd`
LOG_DIR=$NCL_DIR/log
mkdir -p $LOG_DIR
# base_dir="/store/carl/era5/untar/1979/"
base_dir="/store/carl/era5/untar/"
for file in $(find $base_dir -type f -print)
do
filename=`basename $file`
JOB_NAME=$filename
echo $JOB_NAME `date`
LOG_FILE=$LOG_DIR/$JOB_NAME.log
ERR_FILE=$LOG_DIR/$JOB_NAME.err
NCL_OPTION="in_path=\"$file\""
rm $LOG_FILE $ERR_FILE
bsub \
-J $JOB_NAME \
-o $LOG_FILE \
-e $ERR_FILE \
-q $QUEUE \
-n 1 -R "span[hosts=1]" -W 48:00 \
-sp 75 \
~/template/run_ncl.sh $NCL_DIR $NCL_SCRIPT "$NCL_OPTION"
done
#!/bin/bash --login
#if [ -z "$1" ]; then
# QUEUE=`pick_queue`
#else
# QUEUE=$1
#fi
QUEUE=allnodes
echo $QUEUE
RUN_DIR=`pwd`
LOG_DIR=$RUN_DIR/log
mkdir -p $LOG_DIR
script_name=untar_years
for year in {1979..2021}
do
JOB_NAME=$script_name"_"$year
echo $JOB_NAME `date`
LOG_FILE=$LOG_DIR/$JOB_NAME.log
ERR_FILE=$LOG_DIR/$JOB_NAME.err
rm $LOG_FILE $ERR_FILE
bsub \
-J $JOB_NAME \
-o $LOG_FILE \
-e $ERR_FILE \
-q $QUEUE \
-n 1 -R "span[hosts=1]" -W 48:00 \
-sp 75 \
$RUN_DIR/$script_name".sh" $year
done
load "$CJS_NCL_LIB/lib.time.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/cd_inv_string.ncl"
begin
if .not.isvar("in_path") then
in_path = "/store/carl/era5/untar/1979/e5.oper.an.sfc.128_165_10u.ll025sc.1979010100_1979013123.grb.spasub.schreck555709"
end if
print(tostring(in_path))
path_fields = str_split(in_path, "._")
var_name = path_fields(6)
print(tostring(var_name))
if var_name.eq."10v" then
grib_var = "10V_GDS0_SFC"
elseif var_name.eq."10u" then
grib_var = "10U_GDS0_SFC"
end if
start_date_string = path_fields(8)
end_date_string = path_fields(9)
start_date = cd_inv_string(start_date_string, "%Y%N%D%H")
end_date = cd_inv_string(end_date_string, "%Y%N%D%H")
print(cd_string(start_date,"") + " -- " + cd_string(end_date,""))
time = new_time_coord(start_date, end_date, 1)
in_file = addfile(in_path + ".grb", "r")
data = in_file->$grib_var$
data!0 = "time"
data!1 = "lat"
data!2 = "lon"
data&time := time
printVarSummary(data)
printMinMax(data, True)
out_data = calculate_daily_values(data, "avg", 0, False)
printVarSummary(out_data)
printMinMax(out_data, True)
out_dir = getenv("DATA_DIR") + "/era5/daily/" + var_name + "/"
system("mkdir -p " + out_dir)
out_path = out_dir + "/" + var_name \
+ "." + cd_string(start_date, "%Y%N%D") \
+ "." + cd_string(end_date, "%Y%N%D") \
+ ".nc"
system("rm -f " + out_path)
out_file = addfile(out_path, "c")
out_file->$var_name$ = out_data
print("done")
end
;=========================================================================================
;
; Code to download surface data (ERA5) from NCAR RDA via opendap
; The data are also averaged over all hours of the day to get daily averaged values
;
; NCSU Large Scale and Tropical Dynamics
;
; Jan 26, 2022
;
;=========================================================================================
; For pressure levels, see
; https://rda.ucar.edu/thredds/catalog/files/g/ds633.0/e5.oper.an.pl/catalog.html
begin
;========================================================================================
year1=2021
year2=2021
month1=12
month2=12
outDir = "~/data/era5/incoming/"
minLat = -30
maxLat = 30
;========================================================================================
; ;
; choose variables below by uncommenting the relevant block
;=====U10
; define info for the variable that we need
; varName ="VAR_10U" ; this the variable name of the data in the netcdf file
; varid = "165" ; ERA5 id for the variable
; variable = "10u" ; short variable name that is part of the data file name
; ; we will rename varName with this in output data
; vartype = "sc" ; code for variable type
;=====V10
; define info for the variable that we need
varName ="VAR_10V" ; this the variable name of the data in the netcdf file
varid = "166" ; ERA5 id for the variable
variable = "10v" ; short variable name that is part of the data file name
; we will rename varName with this in output data
vartype = "sc" ; code for variable type
;======================================================================================
; No further user edits below
;======================================================================================
urlBase = "https://rda.ucar.edu/thredds/dodsC/files/g/ds633.0/e5.oper.an.sfc/"
varSpec = "/e5.oper.an.sfc.128_" + varid + "_" + variable + ".ll025" + vartype + "."
TimeUnits = "days since 1858-11-17 00:00:00"
do year = year1,year2
do month = month1,month2
day = 1
hour = 0
mins = 0
secs = 0
dateStart = dble2flt(cd_inv_calendar(year,month,day,hour,mins,secs,TimeUnits, 0))
dimo = days_in_month(year,month)
dateEnd = dble2flt(cd_inv_calendar(year,month,dimo,hour,mins,secs,TimeUnits, 0))
suffix = cd_calendar(dateStart,-2) + "00_" + cd_calendar(dateEnd,-2) + "23.nc"
url = urlBase + cd_calendar(dateStart,-1) + varSpec + suffix
print(" " + url)
exists = isfilepresent(url)
if(.not.exists) then
print("OPeNDAP test unsuccessful.")
print("Either the file doesn't exist, or NCL does")
print("not have OPeNDAP cabilities on this system.")
else
f = addfile(url,"r")
vnames = getfilevarnames(f)
print(vnames)
end if
ndim = 0
opt = False ; no need to set True as we should have no missing values
dat := f->$varName$(:,{minLat:maxLat},:)
; dat := calculate_daily_values(f->$varName$(:,{minLat:maxLat},:),"avg", ndim, opt)
printVarSummary(dat)
printMinMax(dat, True)
print(dat(0,:,:))
; Write out to file
fileName = outDir + "era5_sfc" + "_" + cd_calendar(dateStart,-2) + ".nc"
if (fileexists(fileName) ) then
print ("outfile exists")
outFile = addfile( fileName, "w" )
else
;setfileoption("nc", "format", "NETCDF4")
outFile = addfile( fileName, "c" )
end if
outFile->$variable$=dat
print ("wrote " + variable)
end do
end do
end
import xarray as xr
from glob import glob
import os
# file_string = '*2021*nc'
file_string = '*nc'
u_path = f'/store/carl/era5/daily/10u/{file_string}'
v_path = f'/store/carl/era5/daily/10v/{file_string}'
u_ds = xr.open_mfdataset(sorted(glob(u_path)))
print(u_ds)
# print(f'u min:{u_ds["10u"].min().values:.2f} u max: '
# f'{u_ds["10u"].max().values:.2f}')
v_ds = xr.open_mfdataset(sorted(glob(v_path)))
print(v_ds)
# print(f'v min:{v_ds["10v"].min().values:.2f} v max: '
# f'{v_ds["10v"].max().values:.2f}')
wind_speed = xr.ufuncs.sqrt(u_ds['10u']**2 + v_ds['10v']**2)
wind_speed.attrs = u_ds['10u'].attrs
wind_speed.name = 'wind_speed'
wind_speed.attrs['long_name'] = 'scalar wind speed'
print(wind_speed)
# print(f'speed min:{wind_speed.min().values:.2f} speed max: '
# f'{wind_speed.max().values:.2f}')
encoding={'time': {
'units': 'days since 1800-01-01 00:00:00',
'dtype': 'double'
}}
out_path = '/store/carl/era5/'
wind_speed.to_netcdf(f'{out_path}/wind_speed.total.nc', encoding=encoding)
u_ds['10u'].to_netcdf(f'{out_path}/10u.total.nc', encoding=encoding)
v_ds['10v'].to_netcdf(f'{out_path}/10v.total.nc', encoding=encoding)
print('done')
#!/bin/bash --login
if [ -z "$1" ]; then
year=1979
else
year=$1
fi
base_dir=$DATA_DIR/era5/
tar_dir=$base_dir/tar/$year
untar_dir=$base_dir/untar/$year
mkdir $tar_dir
mv $base_dir/incoming/*.$year* $tar_dir/
mkdir $untar_dir
files="$tar_dir/*"
for file in $files
do
echo $file
tar xfk $file -C $untar_dir
done
......@@ -26,28 +26,29 @@ begin
; These are some parameters that could be useful to have up top
if( .not.isvar("varName") ) then
varName = "u850"
varName = "10u"
end if
standardizing = False
basePath = "~/data/cfsr/"
; basePath = "~/data/cfsr/"
; basePath = "~/data/olr/current/"
basePath = "~/data/era5/"
; basePath = "/home/carl/data/nclimgrid/beta/merged/"
if( True ) then
if( False ) then
pathIn = basePath + "total/" + varName + ".total.nc"
pathClim = basePath + "clim/" + varName + ".clim.1991.2020.nc"
pathAnom = basePath + "anom/" + varName + ".anom.nc"
pathClim = basePath + "clim/" + varName + ".clim.1979.2021.nc"
pathAnom = basePath + "anom/" + varName + ".anom.1979.2021.nc"
pathStd = basePath + "std/" + varName + ".std.nc"
else
pathIn = basePath + varName + ".total.nc"
pathClim = basePath + varName + ".clim.1991.2020.nc"
pathAnom = basePath + varName + ".anom.nc"
pathClim = basePath + varName + ".clim.1979.2021.nc"
pathAnom = basePath + varName + ".anom.1979.2021.nc"
pathStd = basePath + varName + ".std.nc"
end if
climStart = 1991001
climEnd = 2020366
climStart = 1979001
climEnd = 2021366
; Open the input files
fin = addfile( pathIn, "r" )
......@@ -108,9 +109,9 @@ begin
print_clock( "Writing the standardized anomalies... " )
system( "rm " + pathStd )
fStd = addfile( pathStd, "c" )
fStd->time = time
fStd->lat = lat
fStd->lon = lon
; fStd->time = time
; fStd->lat = lat
; fStd->lon = lon
fStd->$varName$ = stdData
delete(fStd)
end if
......
......@@ -15,13 +15,13 @@ mkdir -p $LOG_DIR
EMPTY=\"\"
#for VAR_NAME in uShear vShear shear pwat u v chi psi
#for VAR_NAME in uwnd925 uShear vShear shear
for VAR_NAME in u v
for VAR_NAME in 10u 10v wind_speed
do
if [[ $VAR_NAME =~ ^(shear|uShear|vShear|pwat)$ ]]; then
# if [[ $VAR_NAME =~ ^(shear|uShear|vShear|pwat)$ ]]; then
ALL_LEVELS=( $EMPTY )
else
ALL_LEVELS=( 850 200 )
fi
# else
# ALL_LEVELS=( 850 200 )
# fi
for LEVEL in "${ALL_LEVELS[@]}"
do
if [ $LEVEL = $EMPTY ]; then
......@@ -42,7 +42,7 @@ do
-N \
-o $LOG_FILE \
-q $QUEUE \
-n 24 -R "span[hosts=1]" -W 48:00 \
-n 10 -R "span[hosts=1]" -W 48:00 \
-sp 10 \
~/template/run_ncl.sh $NCL_DIR $NCL_SCRIPT "$NCL_OPTION"
......
......@@ -7,26 +7,28 @@ fi
echo $QUEUE
NCL_SCRIPT=filter_waves
VARNAME="u850"
FILTNAME="k09"
NCL_DIR=`pwd`
LOG_DIR=$NCL_DIR/log
mkdir -p $LOG_DIR
JOB_NAME=$NCL_SCRIPT"."$VARNAME"."$FILTNAME
echo $JOB_NAME `date`
LOG_FILE=$LOG_DIR/$JOB_NAME.log
NCL_OPTION="varName=\"$VARNAME\" filtName=\"$FILTNAME\""
for VARNAME in olr 10u 10v wind_speed
do
JOB_NAME=$NCL_SCRIPT"."$VARNAME"."$FILTNAME
echo $JOB_NAME `date`
LOG_FILE=$LOG_DIR/$JOB_NAME.log
NCL_OPTION="varName=\"$VARNAME\" filtName=\"$FILTNAME\""
rm $LOG_FILE
bsub \
-J $JOB_NAME \
-N \
-o $LOG_FILE \
-q $QUEUE \
-n 1 -W 24:00 \
-sp 10 \
~/template/run_ncl.sh $NCL_DIR $NCL_SCRIPT "$NCL_OPTION"
rm $LOG_FILE
bsub \
-J $JOB_NAME \
-N \
-o $LOG_FILE \
-q $QUEUE \
-n 1 -W 24:00 \
-sp 10 \
~/template/run_ncl.sh $NCL_DIR $NCL_SCRIPT "$NCL_OPTION"
done
......@@ -42,25 +42,26 @@ begin
end if
print( varName )
basePath = "~/data/cfsr/"
; basePath = "~/data/nasa_ymc/fcst_verif/"
; basePath = "~/data/cfsr/"
basePath = "~/data/cygnss/kelvin/"
; basePath = "~/data/olr/current/"
; pathIn = basePath + ".std.nc"
; pathOut = basePath + ".waves.std.nc"
; pathIn = basePath + "trmm3b42.anom.nc"
; pathOut = basePath + "trmm3b42.wig.nc"
; pathOut = basePath + "trmm3b42." + filtName + ".nc"
pathIn = basePath + varName + ".anom.nc"
pathOut = basePath + varName + "." + filtName + ".nc"
pathIn = basePath + "anom/" + varName + ".anom.nc"
pathOut = basePath + "waves/" + varName + ".anom.waves.nc"
pathIn = basePath + varName + ".anom.1979.2021.nc"
pathOut = basePath + varName + ".kelvin.nc"
; pathOut = basePath + varName + "." + filtName + ".nc"
; pathIn = basePath + "anom/" + varName + ".anom.nc"
; pathOut = basePath + "waves/" + varName + ".anom.waves.nc"
calcHigh = False
calcLow = True
calcLow = False
calcMrg = False
calcMjo = True
calcMjo = False
calcKelvin = True
calcEr = True
calcEr = False
calcMtd = False
calcTd = False
calcWig = False
......
......@@ -26,10 +26,10 @@ begin
doCalcAnom = True
; basePath = "~/data/cfsr/monthly/" + varName
; basePath = "~/data/ersst/v5/sst.mnmean"
basePath = "~/data/oisst/oisst.month"
inPath = basePath + ".total.nc"
; inPath = basePath + ".nc"
basePath = "~/data/ersst/v5/sst.mnmean"
; basePath = "~/data/oisst/oisst.month"
; inPath = basePath + ".total.nc"
inPath = basePath + ".nc"
climPath = basePath + ".clim.1991.2020.nc"
anomPath = basePath + ".anom.nc"
stdPath = basePath + ".std.nc"
......
......@@ -24,6 +24,12 @@ begin ; main
print_clock( "Averaging" )
latWgtData = NewCosWeight( inData(:,{-30:30},:) )
globalMean = dim_avg_n_Wrap( latWgtData, (/ 1, 2 /) )
wgt = (/1.0, 3.0, 5.0, 3.0, 1.0/)
wgt = wgt/sum(wgt) ; normalize
; globalMean = wgt_runave_n_Wrap( globalMean, wgt, 1, 0 )
globalMean = runave_n_Wrap( globalMean, 11*12+1, 1, 0 )
print(globalMean)
printVarSummary(globalMean)
outData = inData - conform( inData, globalMean, 0 )
copy_VarMeta( inData, outData )
......
......@@ -23,22 +23,24 @@ begin ; main
; These are some parameters that could be useful to have up top
window = 96 * 2
overlap = 64 * 2
obsPerDay = 1
varName = "olr"
sensor = "hirs"
basePath = "/home/carl/data/cygnss/kelvin/"
timeUnits = "days since 1800-01-01 00:00:00"
minTime = cd_inv_calendar( 2001, 01, 01, 00, 0, 0, timeUnits, 0 )
maxTime = cd_inv_calendar( 2010, 12, 31, 18, 0, 0, timeUnits, 0 )
; timeUnits = "days since 1800-01-01 00:00:00"
; minTime = cd_inv_calendar( 1979, 01, 01, 00, 0, 0, timeUnits, 0 )
; maxTime = cd_inv_calendar( 2021, 12, 31, 18, 0, 0, timeUnits, 0 )
print( "Reading the data..." )
basePath = "/home/carl/data/olr/compare/" + sensor + "/"
inPath = basePath + "olr.anom.nc"
outPath = basePath + "olr.wk99.2000s.nc"
inPath = basePath + varName + ".anom.1979.2021.nc"
outPath = basePath + varName + ".wk99.nc"
inFile = addfile( inPath, "r" )
timeUnits = inFile->time@units
print( (/ timeUnits /) )
data = inFile->olr({minTime:maxTime },{-15:15},:)
obsPerDay = 1
; timeUnits = inFile->time@units
; print( (/ timeUnits /) )
; data = inFile->$varName$({minTime:maxTime },{-15:15},:)
data = inFile->$varName$(:,{-15:15},:)
; data@_FillValue = 0
printVarSummary( data )
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment