change the river discharge source point of Fraser and create new Fraser River flow file
from __future__ import division
from salishsea_tools import rivertools
from salishsea_tools import nc_tools
import numpy as np
import matplotlib.pyplot as plt
import netCDF4 as nc
import arrow
import numpy.ma as ma
import sys
sys.path.append('/ocean/klesouef/meopar/tools/I_ForcingFiles/Rivers')
%matplotlib inline
filename = '/ocean/jieliu/research/meopar/nemo-forcing/rivers/rivers_month.nc'
clim_rivers = nc.Dataset(filename, 'r')
nc_tools.show_dimensions(clim_rivers)
nc_tools.show_variables(clim_rivers)
criverflow = clim_rivers.variables['rorunoff']
# get other variables so we can put them in new files
lat = clim_rivers.variables['nav_lat']
lon = clim_rivers.variables['nav_lon']
riverdepth = clim_rivers.variables['rodepth']
<type 'netCDF4.Dimension'>: name = 'x', size = 398 <type 'netCDF4.Dimension'>: name = 'y', size = 898 <type 'netCDF4.Dimension'> (unlimited): name = 'time_counter', size = 12 [u'nav_lat', u'nav_lon', u'time_counter', u'rorunoff', u'rodepth', u'rotemper']
rivertype = 'constant' ## monthly or constant(yearly)
if rivertype == 'monthly':
fluxfile = nc.Dataset('/ocean/sallen/allen/research/MEOPAR/nemo-forcing/rivers/Salish_allrivers_monthly.nc','r')
#inialise the runoff and run_depth arrays
runoff, run_depth, run_temp = rivertools.init_runoff_array_monthly()
#get river fluxes from netcdf file
if rivertype == 'constant':
fluxfile = nc.Dataset('/ocean/sallen/allen/research/MEOPAR/nemo-forcing/rivers/Salish_allrivers_cnst.nc','r')
#inialise the runoff and run_depth arrays
runoff, run_depth, run_temp = rivertools.init_runoff_array()
#list of watersheds we are including
names = ['skagit','fraser','evi_n','howe','bute','puget','jdf','evi_s','jervis','toba']
for name in range(0,len(names)):
watershedname = names[name]
Flux = fluxfile.variables[watershedname][:]
if rivertype == 'constant':
Flux = float(Flux)
runoff_orig = np.copy(runoff)
runoff, run_depth, run_temp = rivertools.put_watershed_into_runoff(rivertype,
watershedname, Flux, runoff, run_depth, run_temp)
if rivertype == 'monthly':
rivertools.check_sum_monthly(runoff_orig, runoff, Flux)
if rivertype == 'constant':
rivertools.check_sum(runoff_orig, runoff, Flux)
skagit has 10 rivers (908.05509995059549, 937.3491821289062) fraser has 10 rivers (3578.4171797583772, 3549.736083984375) evi_n has 21 rivers (255.62751641603765, 638.462158203125) howe has 2 rivers (588.99650864276578, 572.8439331054688) bute has 3 rivers (609.26456547060798, 550.5059204101562) puget has 43 rivers (480.71352669585275, 502.76959228515625) jdf has 27 rivers (399.59934187151299, 409.8137512207031) evi_s has 17 rivers (330.71661134751884, 329.566162109375) jervis has 17 rivers (307.47534529284559, 296.90533447265625) toba has 1 rivers (285.36386496863395, 270.2106018066406)
print run_depth[414,334]
-1.0
if rivertype == 'monthly':
nemo = nc.Dataset('/ocean/jieliu/research/meopar/river-treatment/rivers_month_edit.nc', 'w')
nemo.description = 'Monthly Averages, All Rivers, modify on depth and runoff grid point'
# dimensions
nemo.createDimension('x', 398)
nemo.createDimension('y', 898)
nemo.createDimension('time_counter', None)
# variables
# latitude and longitude
nav_lat = nemo.createVariable('nav_lat','float32',('y','x'),zlib=True)
nav_lat = lat
x = nemo.createVariable('nav_lon','float32',('y','x'),zlib=True)
nav_lon = lon
# time
time_counter = nemo.createVariable('time_counter', 'float32', ('time_counter'),zlib=True)
time_counter.units = 'non-dim'
time_counter[0:12] = range(1,13)
# runoff
rorunoff = nemo.createVariable('rorunoff', 'float32', ('time_counter','y','x'), zlib=True)
rorunoff._Fillvalue = 0.
rorunoff._missing_value = 0.
rorunoff._units = 'kg m-2 s-1'
rorunoff[0:12,:] = runoff
# depth
rodepth = nemo.createVariable('rodepth','float32',('y','x'),zlib=True)
rodepth._Fillvalue = -1.
rodepth.missing_value = -1.
rodepth.units = 'm'
rodepth[:] = run_depth[0,:,:]
# temperature
rotemper = nemo.createVariable('rotemper','float32',('time_counter','y','x'),zlib=True)
rotemper._Fillvalue = -99.
rotemper.missing_value = -99.
rotemper.units = 'deg C'
rotemper[0:12,:] = run_temp
nemo.close()
if rivertype == 'constant':
nemo = nc.Dataset('/ocean/jieliu/research/meopar/river-treatment/rivers_cnst_edit.nc', 'w')
nemo.description = 'Constant Yearly Average, All Rivers, modify on depth and runoff grid point'
# dimensions
nemo.createDimension('x', 398)
nemo.createDimension('y', 898)
nemo.createDimension('time_counter', None)
# variables
# latitude and longitude
nav_lat = nemo.createVariable('nav_lat','float32',('y','x'),zlib=True)
nav_lat = lat
x = nemo.createVariable('nav_lon','float32',('y','x'),zlib=True)
nav_lon = lon
# time
time_counter = nemo.createVariable('time_counter', 'float32', ('time_counter'),zlib=True)
time_counter.units = 'non-dim'
time_counter[0] = 1
# runoff
rorunoff = nemo.createVariable('rorunoff', 'float32', ('time_counter','y','x'), zlib=True)
rorunoff._Fillvalue = 0.
rorunoff._missing_value = 0.
rorunoff._units = 'kg m-2 s-1'
rorunoff[0,:] = runoff
# depth
rodepth = nemo.createVariable('rodepth','float32',('y','x'),zlib=True)
rodepth._Fillvalue = -1.
rodepth.missing_value = -1.
rodepth.units = 'm'
rodepth[:] = run_depth
nemo.close()
# Constant and data ranges etc
year = 2015
smonth = 06
emonth = 06
startdate = arrow.get(year,smonth,14)
enddate = arrow.get(year,emonth,14)
print startdate, enddate
2015-06-14T00:00:00+00:00 2015-06-14T00:00:00+00:00
# get Fraser Flow data
filename = '/data/dlatorne/SOG-projects/SOG-forcing/ECget/Fraser_flow'
fraserflow = np.loadtxt(filename)
print fraserflow
[[ 2013. 10. 17. 1709.085] [ 2013. 10. 18. 1676.078] [ 2013. 10. 19. 1653.68 ] ..., [ 2015. 6. 28. 4622.91 ] [ 2015. 6. 29. 4586.816] [ 2015. 6. 30. 4717.231]]
#Fraser watershed
pd = rivertools.get_watershed_prop_dict('fraser')
totalfraser = (pd['Fraser1']['prop'] + pd['Fraser2']['prop'] +
pd['Fraser3']['prop'] + pd['Fraser4']['prop'])
fraser has 10 rivers
# Climatology, Fraser Watershed
fluxfile = nc.Dataset('/ocean/jieliu/research/meopar/nemo-forcing/rivers/Salish_allrivers_monthly.nc','r')
climFraserWaterShed = fluxfile.variables['fraser'][:]
# Fraser River at Hope Seasonal Climatology (found in matlab using Mark's mean daily data)
climFraseratHope = (931, 878, 866, 1814, 4097, 6970, 5538, 3539, 2372, 1937, 1595, 1119)
NonHope = climFraserWaterShed - climFraseratHope
otherratio = 0.016
fraserratio = 1-otherratio
nonFraser = (otherratio * climFraserWaterShed.sum()/NonHope.sum()) * NonHope
afterHope = NonHope - nonFraser
print pd['Fraser1']['i'],pd['Fraser1']['j']
418 397
def calculate_daily_flow(r,criverflow):
'''interpolate the daily values from the monthly values'''
print r.day, r.month
if r.day < 16:
prevmonth = r.month-1
if prevmonth == 0:
prevmonth = 12
nextmonth = r.month
else:
prevmonth = r.month
nextmonth = r.month + 1
if nextmonth == 13:
nextmonth = 1
fp = r - arrow.get(year,prevmonth,15)
fn = arrow.get(year,nextmonth,15) - r
ft = fp+fn
fp = fp.days/ft.days
fn = fn.days/ft.days
print ft, fp, fn
driverflow = fn*criverflow[prevmonth-1] + fp*criverflow[nextmonth-1]
return driverflow
def write_file(r,flow,lat,lon,riverdepth):
''' given the flow and the riverdepth and the date, write the nc file'''
directory = '.'
# set up filename to follow NEMO conventions
filename = 'RFraserCElse_y'+str(year)+'m'+'{:0=2}'.format(r.month)+'d'+'{:0=2}'.format(r.day)+'.nc'
# print directory+'/'+filename
nemo = nc.Dataset(directory+'/'+filename, 'w')
nemo.description = 'Real Fraser Values, Daily Climatology for Other Rivers'
# dimensions
ymax, xmax = lat.shape
nemo.createDimension('x', xmax)
nemo.createDimension('y', ymax)
nemo.createDimension('time_counter', None)
# variables
# latitude and longitude
nav_lat = nemo.createVariable('nav_lat','float32',('y','x'),zlib=True)
nav_lat = lat
x = nemo.createVariable('nav_lon','float32',('y','x'),zlib=True)
nav_lon = lon
# time
time_counter = nemo.createVariable('time_counter', 'float32', ('time_counter'),zlib=True)
time_counter.units = 'non-dim'
time_counter[0:1] = range(1,2)
# runoff
rorunoff = nemo.createVariable('rorunoff', 'float32', ('time_counter','y','x'), zlib=True)
rorunoff._Fillvalue = 0.
rorunoff._missing_value = 0.
rorunoff._units = 'kg m-2 s-1'
rorunoff[0,:] = flow
# depth
rodepth = nemo.createVariable('rodepth','float32',('y','x'),zlib=True)
rodepth._Fillvalue = -1.
rodepth.missing_value = -1.
rodepth.units = 'm'
rodepth = riverdepth
nemo.close()
return
def fraser_correction(pd, fraserflux, r, afterHope, NonFraser, fraserratio, otherratio,
runoff):
''' for the Fraser Basin only, replace basic values with the new climatology after Hope and the
observed values for Hope. Note, we are changing runoff only and not using/changing river
depth '''
for key in pd:
if "Fraser" in key:
flux = calculate_daily_flow(r,afterHope) + fraserflux
subarea = fraserratio
else:
flux = calculate_daily_flow(r,NonFraser)
subarea = otherratio
river = pd[key]
runoff = rivertools.fill_runoff_array(flux*river['prop']/subarea,river['i'],
river['di'],river['j'],river['dj'],river['depth'],
runoff,np.empty_like(runoff))[0]
return runoff
##open climatolgy file with modified fresh water point source
clim_rivers_edit = nc.Dataset('rivers_month_edit.nc','r' )
criverflow_edit = clim_rivers_edit.variables['rorunoff']
for r in arrow.Arrow.range('day', startdate, enddate):
print r
driverflow = calculate_daily_flow(r, criverflow_edit)
storeflow = calculate_daily_flow(r, criverflow_edit)
step1 = fraserflow[fraserflow[:,0] == r.year]
step2 = step1[step1[:,1] == r.month]
step3 = step2[step2[:,2] == r.day]
# print r.year, r.month, r.day, step3[0,3]
runoff = fraser_correction(pd, step3[0,3] , r, afterHope, nonFraser, fraserratio, otherratio,
driverflow)
write_file(r,runoff,lat,lon,riverdepth)
ig = 418
jg = 397
print criverflow_edit[7:10,418,397], driverflow[ig,jg]
print storeflow[ig,jg], driverflow[ig,jg]
ig = 351; jg = 345
print storeflow[ig,jg], driverflow[ig,jg]
ig = 749; jg=123
print storeflow[ig,jg], driverflow[ig,jg]
# jan 0, feb 1, mar 2, apr 3, may 4, jun 5
# jul 6, aug 7, sep 8
2015-06-14T00:00:00+00:00 14 6 31 days, 0:00:00 0.967741935484 0.0322580645161 14 6 31 days, 0:00:00 0.967741935484 0.0322580645161 14 6 31 days, 0:00:00 0.967741935484 0.0322580645161 14 6 31 days, 0:00:00 0.967741935484 0.0322580645161 14 6 31 days, 0:00:00 0.967741935484 0.0322580645161 14 6 31 days, 0:00:00 0.967741935484 0.0322580645161 14 6 31 days, 0:00:00 0.967741935484 0.0322580645161 14 6 31 days, 0:00:00 0.967741935484 0.0322580645161 14 6 31 days, 0:00:00 0.967741935484 0.0322580645161 14 6 31 days, 0:00:00 0.967741935484 0.0322580645161 14 6 31 days, 0:00:00 0.967741935484 0.0322580645161 14 6 31 days, 0:00:00 0.967741935484 0.0322580645161 [ 14.77230835 10.16569614 8.63957405] 26.5695 27.8402 26.5695 0.0237942 0.0143087 0.203718 0.203718