diff --git a/docs/images/partial-cells.png b/docs/images/partial-cells.png new file mode 100644 index 0000000..a83a6de Binary files /dev/null and b/docs/images/partial-cells.png differ diff --git a/examples/rivers/make_river_clim.py b/examples/rivers/make_river_clim.py index 598053d..ff12052 100644 --- a/examples/rivers/make_river_clim.py +++ b/examples/rivers/make_river_clim.py @@ -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') diff --git a/pycnal/pycnal/grid.py b/pycnal/pycnal/grid.py index d1dc500..315e5a9 100644 --- a/pycnal/pycnal/grid.py +++ b/pycnal/pycnal/grid.py @@ -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). @@ -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 @@ -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] @@ -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): @@ -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 @@ -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): """ @@ -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() @@ -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') @@ -359,7 +359,7 @@ 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, \ @@ -367,9 +367,9 @@ def get_ROMS_hgrid(gridid): 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 @@ -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. """ @@ -418,8 +418,11 @@ 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 @@ -427,7 +430,7 @@ def get_ROMS_vgrid(gridid, zeta=None): vgrid = z_coordinate(h, depth, N) else: - raise ValueError('Unknow grid type') + raise ValueError('Unknown grid type') return vgrid @@ -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. """ @@ -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' @@ -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) diff --git a/pycnal/pycnal/vgrid.py b/pycnal/pycnal/vgrid.py index 70fdbe7..dada0e1 100644 --- a/pycnal/pycnal/vgrid.py +++ b/pycnal/pycnal/vgrid.py @@ -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): @@ -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): @@ -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]) / \ diff --git a/pycnal_toolbox/pycnal_toolbox/remapping.py b/pycnal_toolbox/pycnal_toolbox/remapping.py index da25184..9971f35 100644 --- a/pycnal_toolbox/pycnal_toolbox/remapping.py +++ b/pycnal_toolbox/pycnal_toolbox/remapping.py @@ -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) @@ -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') @@ -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