Skip to content

Commit

Permalink
Merge branch 'jpender'
Browse files Browse the repository at this point in the history
  • Loading branch information
kshedstrom committed Jul 18, 2017
2 parents f3be425 + 2a3ae2a commit 59cd388
Show file tree
Hide file tree
Showing 5 changed files with 48 additions and 44 deletions.
Binary file added docs/images/partial-cells.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
2 changes: 1 addition & 1 deletion examples/rivers/make_river_clim.py
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@
import pycnal_toolbox


# load 2-dimentional discharge data
# load 2-dimentional discharge data
print('Load discharge data')
nc_data = netCDF.Dataset('CI_runoff.nc', 'r')
nc_rivers = netCDF.Dataset('Cook_Inlet_rivers.nc', 'a')
Expand Down
64 changes: 33 additions & 31 deletions pycnal/pycnal/grid.py
Original file line number Diff line number Diff line change
Expand Up @@ -43,10 +43,10 @@ class ROMS_gridinfo(object):
There are two ways to define the grid information. If grid_file
and hist_file are not passed to the object when it is created, the
information is retrieved from gridid.txt.
To add new grid please edit your gridid.txt. You need to define
an environment variable PYCNAL_GRIDID_FILE pointing to your
gridid.txt file. Just copy an existing grid and modify the
information is retrieved from gridid.txt.
To add new grid please edit your gridid.txt. You need to define
an environment variable PYCNAL_GRIDID_FILE pointing to your
gridid.txt file. Just copy an existing grid and modify the
definition accordingly to your case (Be carefull with
space and blank line).
Expand All @@ -69,7 +69,7 @@ def __init__(self, gridid,grid_file=None,hist_file=None):
#the grid and history files from the model
self.id = gridid
self._get_grid_info(grid_file,hist_file)

#now save the data in the dictionary, so we don't need to get it again
gridid_dictionary[gridid]=self

Expand Down Expand Up @@ -107,7 +107,7 @@ def _get_grid_info(self,grid_file,hist_file):
line_nb = line_nb + 1

if info == []:
raise ValueError('Unknow gridid. Please check your gridid.txt file')
raise ValueError('Unknown gridid. Please check your gridid.txt file')

if info[4] == 'roms':
self.name = info[1]
Expand All @@ -119,7 +119,7 @@ def _get_grid_info(self,grid_file,hist_file):
self.theta_b = np.float(info[7])
self.Tcline = np.float(info[8])

elif info[4] == 'z':
elif info[4] == 'z':
nline = len(info)
dep = info[5]
for line in range(6,nline):
Expand All @@ -133,14 +133,14 @@ def _get_grid_info(self,grid_file,hist_file):
self.depth = dep

else:
raise ValueError('Unknow grid type. Please check your gridid.txt file')
raise ValueError('Unknown grid type. Please check your gridid.txt file')

else: #lets get the grid information from the history and grid files
#print 'CJMP> getting grid info from ROMS history and grid files'
assert type(grid_file)!=type(None), 'if specify history file you must specify grid file'
assert type(hist_file)!=type(None), 'if specify grid file you must specify history file'

#open history file and get necessary grid information from it.
#open history file and get necessary grid information from it.
hist=netCDF.Dataset(hist_file,'r')

#put data into ROMS_gridinfo object
Expand All @@ -164,7 +164,7 @@ def _get_grid_info(self,grid_file,hist_file):
self.theta_s=np.float(hist.variables['theta_s'][:])
self.theta_b=np.float(hist.variables['theta_b'][:])
self.Tcline=np.float(hist.variables['Tcline'][:])


def print_ROMS_gridinfo(gridid):
"""
Expand Down Expand Up @@ -202,7 +202,7 @@ def list_ROMS_gridid():
data = open(gridid_file,'r')
lines = data.readlines()
data.close()

gridid_list = []
for line in lines:
s = line.split()
Expand All @@ -227,7 +227,7 @@ def get_ROMS_hgrid(gridid):
#Check for cartesian or geographical grid
spherical = nc.variables['spherical'][:]

#Get horizontal grid
#Get horizontal grid
if ((spherical == 0) or (spherical == 'F')):
#cartesian grid
print('Load cartesian grid from file')
Expand Down Expand Up @@ -359,17 +359,17 @@ def get_ROMS_hgrid(gridid):
else:
angle = None

#Get geographical grid
#Get geographical grid
hgrd = CGrid_geo(lon_vert, lat_vert, proj, \
lon_rho=lon_rho, lat_rho=lat_rho, \
lon_u=lon_u, lat_u=lat_u, lon_v=lon_v, lat_v=lat_v, \
lon_psi=lon_psi, lat_psi=lat_psi, dx=dx, dy=dy, \
dndx=dndx, dmde=dmde, angle_rho=angle)

#load the mask
try:
try:
hgrd.mask_rho = np.array(nc.variables['mask_rho'][:])
except:
except:
hgrd.mask_rho = np.ones(hgrd.lat_rho.shape)

return hgrd
Expand All @@ -382,11 +382,11 @@ def get_ROMS_vgrid(gridid, zeta=None):
Load ROMS vertical grid object. vgrid is a s_coordinate or
a z_coordinate object, depending on gridid.grdtype.
vgrid.z_r and vgrid.z_w (vgrid.z for a z_coordinate object)
can be indexed in order to retreive the actual depths. The
free surface time serie zeta can be provided as an optional
argument. Note that the values of zeta are not calculated
until z is indexed, so a netCDF variable for zeta may be passed,
even if the file is large, as only the values that are required
can be indexed in order to retreive the actual depths. The
free surface time serie zeta can be provided as an optional
argument. Note that the values of zeta are not calculated
until z is indexed, so a netCDF variable for zeta may be passed,
even if the file is large, as only the values that are required
will be retrieved from the file.
"""

Expand Down Expand Up @@ -418,16 +418,19 @@ def get_ROMS_vgrid(gridid, zeta=None):
vgrid = s_coordinate_2(h, theta_b, theta_s, Tcline, N, hraw=hraw, zeta=zeta)
elif Vtrans == 4:
vgrid = s_coordinate_4(h, theta_b, theta_s, Tcline, N, hraw=hraw, zeta=zeta)
elif Vtrans == 5:
vgrid = s_coordinate_5(h, theta_b, theta_s, Tcline, N, hraw=hraw, zeta=zeta)

else:
raise Warning('Unknow vertical transformation Vtrans')
raise Warning('Unknown vertical transformation Vtrans')

elif gridinfo.grdtype == 'z':
N = gridinfo.N
depth = gridinfo.depth
vgrid = z_coordinate(h, depth, N)

else:
raise ValueError('Unknow grid type')
raise ValueError('Unknown grid type')

return vgrid

Expand All @@ -448,15 +451,15 @@ def get_ROMS_grid(gridid, zeta=None, hist_file=None,grid_file=None):
grid information will be extracted from those files, and gridid
will be used to name that grid for the rest of the python
session.
grd.vgrid is a s_coordinate or
a z_coordinate object, depending on gridid.grdtype.
grd.vgrid.z_r and grd.vgrid.z_w (grd.vgrid.z for a
z_coordinate object) can be indexed in order to retreive the
actual depths. The free surface time serie zeta can be provided
as an optional argument. Note that the values of zeta are not
calculated until z is indexed, so a netCDF variable for zeta may
be passed, even if the file is large, as only the values that
grd.vgrid.z_r and grd.vgrid.z_w (grd.vgrid.z for a
z_coordinate object) can be indexed in order to retreive the
actual depths. The free surface time serie zeta can be provided
as an optional argument. Note that the values of zeta are not
calculated until z is indexed, so a netCDF variable for zeta may
be passed, even if the file is large, as only the values that
are required will be retrieved from the file.
"""

Expand Down Expand Up @@ -486,7 +489,6 @@ def write_ROMS_grid(grd, filename='roms_grd.nc'):

Mm, Lm = grd.hgrid.x_rho.shape


# Write ROMS grid to file
nc = netCDF.Dataset(filename, 'w', format='NETCDF3_64BIT')
nc.Description = 'ROMS grid'
Expand All @@ -498,7 +500,7 @@ def write_ROMS_grid(grd, filename='roms_grd.nc'):
nc.createDimension('xi_u', Lm-1)
nc.createDimension('xi_v', Lm)
nc.createDimension('xi_psi', Lm-1)

nc.createDimension('eta_rho', Mm)
nc.createDimension('eta_u', Mm)
nc.createDimension('eta_v', Mm-1)
Expand Down
17 changes: 8 additions & 9 deletions pycnal/pycnal/vgrid.py
Original file line number Diff line number Diff line change
Expand Up @@ -289,21 +289,18 @@ def __init__(self, h, theta_b, theta_s, Tcline, N, hraw=None, zeta=None):
self.z_r = z_r(self.h, self.hc, self.N, self.s_rho, self.Cs_r, self.zeta, self.Vtrans)
self.z_w = z_w(self.h, self.hc, self.Np, self.s_w, self.Cs_w, self.zeta, self.Vtrans)


def _get_s_rho(self):
lev = np.arange(1,self.N+1,1)
s = -(lev * lev - 2 * lev * self.N + lev + self.N * self.N - self.N) / \
(self.N * self.N - self.N) - \
0.01 * (lev * lev - lev * self.N) / (self.c1 - self.N)
# (self.c1 * self.N * self.N - self.N) - \
self.s_rho = s
lev = np.arange(1, self.N+1) - .5
self.s_rho = -(lev * lev - 2 * lev * self.N + lev + self.N * self.N - self.N) / \
(1.0 * self.N * self.N - self.N) - \
0.01 * (lev * lev - lev * self.N) / (1.0 - self.N)


def _get_s_w(self):
lev = np.arange(0,self.Np,1)
s = -(lev * lev - 2 * lev * self.N + lev + self.N * self.N - self.N) / \
(self.N * self.N - self.N) - \
0.01 * (lev * lev - lev * self.N) / (self.c1 - self.N)
# (self.c1 * self.N * self.N - self.N) - \
self.s_w = s

def _get_Cs_r(self):
Expand Down Expand Up @@ -347,6 +344,7 @@ def __init__(self, h, hc, N, s_rho, Cs_r, zeta, Vtrans):
self.zeta = zeta
self.Vtrans = Vtrans


def __getitem__(self, key):

if isinstance(key, tuple) and len(self.zeta.shape) > len(self.h.shape):
Expand Down Expand Up @@ -412,12 +410,13 @@ def __getitem__(self, key):

ti = zeta.shape[0]
z_w = np.empty((ti, self.Np) + self.h.shape, 'd')

if self.Vtrans == 1:
for n in range(ti):
for k in range(self.Np):
z0 = self.hc * self.s_w[k] + (self.h - self.hc) * self.Cs_w[k]
z_w[n,k,:] = z0 + zeta[n,:] * (1.0 + z0 / self.h)
elif self.Vtrans == 2 or self.Vtrans == 4:
elif self.Vtrans == 2 or self.Vtrans == 4 or self.Vtrans == 5:
for n in range(ti):
for k in range(self.Np):
z0 = (self.hc * self.s_w[k] + self.h * self.Cs_w[k]) / \
Expand Down
9 changes: 6 additions & 3 deletions pycnal_toolbox/pycnal_toolbox/remapping.py
Original file line number Diff line number Diff line change
Expand Up @@ -53,6 +53,7 @@ def remapping(varname, srcfile, wts_files, srcgrd, dstgrd, \
srcgrdz = pycnal.grid.ROMS_Grid(srcgrd.name+'_Z', srcgrd.hgrid, src_zcoord)
dstgrdz = pycnal.grid.ROMS_Grid(dstgrd.name+'_Z', dstgrd.hgrid, dst_zcoord)


# varname argument
if type(varname).__name__ == 'list':
nvar = len(varname)
Expand Down Expand Up @@ -214,14 +215,11 @@ def remapping(varname, srcfile, wts_files, srcgrd, dstgrd, \
else:
src_varz = src_var[nt,jjrange[0]:jjrange[1],iirange[0]:iirange[1]]

# print datetime.datetime.now()
# horizontal interpolation using scrip weights
print('horizontal interpolation using scrip weights')
dst_varz = pycnal.remapping.remap(src_varz, wts_file, \
spval=spval)

# print datetime.datetime.now()

if ndim == 3:
# vertical interpolation from standard z level to sigma
print('vertical interpolation from standard z level to sigma')
Expand All @@ -230,6 +228,11 @@ def remapping(varname, srcfile, wts_files, srcgrd, dstgrd, \
else:
dst_var = dst_varz

if varname[nv] == 'u':
dst_u = dst_var
if varname[nv] == 'v':
dst_v = dst_var

# write data in destination file
print('write data in destination file')
nc.variables[varname[nv]][nctidx] = dst_var
Expand Down

0 comments on commit 59cd388

Please sign in to comment.