pgrid initial development for LO

pgrid/README.md

@@ -0,0 +1,54 @@
+# README for pgrid
+
+### This collection of programs is designed to make gridfiles for ROMS.  It works for analytical and realistic cases, and handles complex tasks like mask editing and smoothing.  It also creates other files associated with rivers, nudging to climatology, and vertical grid parameters, all in the form expected by LO/forcing and ROMS.
+
+---
+
+First edit the gridname and (if needed) a few directory locations at the top of gfun.py.  This gridname will then be used by all subsequent code.
+
+In order to keep track of several choices typically made about a grid, we use "dch," a dict of "default choices":
+```
+dch =  gfun.default_choices(Gr)
+```
+These are initialized in gfun.default_choices, but you typically override some of them in start_grid.py when doing the initial grid specification.  The choices are saved in a pickle file:
+```
+pickle.dump(dch, open(Gr['gdir'] + 'choices.p', 'wb'))
+```
+You can also go back and change things in dch later using Z_edit_dch.py.
+
+Each time you run a piece of code it makes a new grid file with the name altered to indicate what happened.  The names start as: grid_m00_r00_s00_x00.nc with the letters and numbers indicating changes to: mask, river, smoothing, or extras.
+
+You can use plot_grid.py to look at any of the grids.  You can override the default grid by calling it in ipython like:
+- run plot_grid.py -g 'sal0'
+
+Suggested order to run the code, and which dch items are used at each step:
+
+`start_grid.py` (can take awhile for big grids; be patient; e.g. 15 min for cas6)
+- analytical
+- t_dir/t_list
+- use_z_offset/z_offset
+
+`make_mask.py`
+- z_land
+- unmask_coast
+- remove_islands
+
+`edit_mask.py` (to get rid of obvious issues like Lake Washington and river channels)
+
+`carve_rivers.py`
+
+`edit_mask.py` (for real this time, perhaps running many times) NOTE: you can use an optional command line argument -d ## to change "dval" the carving depth used for lines or points from its default value of 5 m.
+
+`smooth_grid.py`
+- use_min_depth/min_depth
+
+`carve_rivers.py` (not needed?)
+
+`smooth_grid.py` (not needed?)
+
+`make_extras.py`
+- min_depth (enforced for whole grid)
+
+`grid_to_LiveOcean.py` (should be run on remote machine for large grids) (also saves S-coordinate info, so be careful what file you choose in the code)
+- nudging_edges
+- nudging_days 

pgrid/gfun.py

@@ -1,22 +1,26 @@
-# -*- coding: utf-8 -*-
 """
 Organizational functions for pgrid.
 """
 
 # **** USER EDIT ********
-#gridname = 'aestus3'
-gridname = 'cas6'
+gridname = 'sal0'; base_gridname = 'cas6'; base_tag = 'v3'
 # **** END USER EDIT ****
 
-import os; import sys
-sys.path.append(os.path.abspath('../../LiveOcean/alpha'))
-import Lfun
-Ldir = Lfun.Lstart()
+from pathlib import Path
+from lo_tools import Lfun
 
-sys.path.append(os.path.abspath('../../LiveOcean/plotting'))
+Ldir = Lfun.Lstart(gridname=base_gridname, tag=base_tag)
 
-dir0 = Ldir['parent']
-pgdir = dir0 + 'ptools_output/pgrid/'
+def gstart():
+    """
+    This returns a dict of Path objects that tell where various things are,
+    or where they should go.
+    """
+    pgdir = Ldir['LOo'] / 'pgrid'
+    gdir = pgdir / gridname # where grid.nc will end up
+    ri_dir = Ldir['LOo'] / 'pre' / 'river' / Ldir['gtag'] / 'tracks'
+    Gr ={'gridname': gridname,'pgdir': pgdir, 'gdir': gdir,'ri_dir': ri_dir}
+    return Gr
 
 def default_choices(Gr, wet_dry=False):
     # Default choices (can override in each case)    
@@ -34,12 +38,13 @@ def default_choices(Gr, wet_dry=False):
     dch['use_z_offset'] = True
     dch['z_offset'] = -1.06
     # specify topography files to use
-    dch['t_dir'] = Gr['dir0'] + 'ptools_data/topo/'    
+    t_dir = Ldir['data'] / 'topo'
+    dch['t_dir'] = Ldir['data'] / 'topo'
     # list of topo files: coarsest to finest
-    dch['t_list'] = ['srtm15/topo15.nc',
-              'cascadia/cascadia_gridded.nc',
-             'psdem/PS_183m.nc',
-             'ttp_patch/TTP_Regional_27m_patch.nc']
+    dch['t_list'] = [t_dir / 'srtm15' / 'topo15.nc',
+              t_dir / 'cascadia' / 'cascadia_gridded.nc',
+             t_dir / 'psdem' / 'PS_183m.nc',
+             t_dir / 'ttp_patch' / 'TTP_Regional_27m_patch.nc']
 
     # MASKING
     # list of existing masks to work from
@@ -67,46 +72,16 @@ def default_choices(Gr, wet_dry=False):
 
     return dch
 
-def gstart(gridname=gridname):
-
-    if gridname in ['aestus1', 'aestus2']:
-        ri_dir = dir0 + 'ptools_output/river/analytical/'
-    else:
-        ri_dir = dir0 + 'ptools_output/river/pnw_all_2016_07/'
-
-    gdir = pgdir + gridname + '/'
-    Gr ={'gridname': gridname, 'dir0': dir0, 'pgdir': pgdir, 'gdir': gdir,
-         'ri_dir': ri_dir}
-    return Gr
-
-def select_file(Gr, using_old_grid=False):
+def select_file(Gr):
     # interactive selection
-    if using_old_grid==True:
-        fn_list = []
-        dir0 = Ldir['parent'] + 'LiveOcean_data/grids/'
-        gn_list = ['cascadia1', 'cascadia2']
-        for gn in gn_list:
-            fn_list.append(dir0 + gn + '/grid.nc')
-    elif using_old_grid==False:
-        print('\n** %s in <<%s>> **\n' % ('Choose file to edit', Gr['gridname']))
-        fn_list_raw = os.listdir(Gr['gdir'])
-        fn_list = []
-        for item in fn_list_raw:
-            if item[-3:] == '.nc':
-                fn_list.append(item)
-        fn_list.sort()
-    Nfn = len(fn_list)
-    fn_dict = dict(zip(range(Nfn), fn_list))
-    for nfn in range(Nfn):
-        print(str(nfn) + ': ' + fn_list[nfn])
-    my_nfn = int(input('-- Input number -- '))
-    fn = fn_dict[my_nfn]
+    fn = choose_item(Gr['gdir'], tag='.nc', itext='** Choose grid from list **')
     return fn
 
 def increment_filename(fn, tag='_m'):
+    fns = str(fn)
     # create the new file name
-    gni = fn.find(tag)
-    new_num = ('00' + str(int(fn[gni+2: gni+4]) + 1))[-2:]
-    fn_new = fn.replace(fn[gni:gni+4], tag + new_num)
-
+    gni = fns.find(tag)
+    new_num = ('00' + str(int(fns[gni+2: gni+4]) + 1))[-2:]
+    fns_new = fns.replace(fns[gni:gni+4], tag + new_num)
+    fn_new = Path(fns_new)
     return fn_new

pgrid/gfun_utility.py

@@ -1,26 +1,11 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
 """
-Created on Wed Mar 29 15:55:45 2017
-
-@author: PM5
-
 Utility functions for pgrid.
 """
 import numpy as np
-import zfun
-import h5py
-import matfun
-import netCDF4 as nc
+import xarray as xr
 import seawater as sw
 
-import os
-import sys
-pth = os.path.abspath('../../LiveOcean/alpha')
-if pth not in sys.path:
-    sys.path.append(pth)
-import zrfun
-import Lfun
+from lo_tools import Lfun, zfun, zrfun
 
 def simple_grid(aa, res):
     dlat = aa[3] - aa[2]
@@ -74,46 +59,18 @@ def stretched_grid(lon_list, x_res_list, lat_list, y_res_list):
     return np.array(plon_list), np.array(plat_list)
 
 def load_bathy_nc(t_fn):
-    ds = nc.Dataset(t_fn)
-    tlon_vec = ds['lon'][:]
-    tlat_vec = ds['lat'][:]
-    tz = ds['z'][:]
+    ds = xr.open_dataset(t_fn)
+    tlon_vec = ds['lon'].values
+    tlat_vec = ds['lat'].values
+    tz = ds['z'].values
     ds.close()
     return tlon_vec, tlat_vec, tz
-
-def load_bathy(t_fn):
-    try:
-        a = h5py.File(t_fn)
-        aa = dict()
-        print(' using h5py')
-        for item in a.keys():
-            aa[item] = a[item][:]
-            #print('    ' + item + ' ' + str(aa[item].shape))
-        # reshaping because of how things are packed
-        tlon_vec = a['lon'][:].flatten()
-        tlat_vec = a['lat'][:].flatten()
-        tz = a['z'][:].T
-        #tmask = a['mask'][:].T
-        # The field "mask" was created in the topo processing
-        # as a matrix of ints with 1=ocean, 0=land.
-        # We will supersede it here with our own masking.
-        a.close()
-    except:
-        # for some reason h5py does not work with, for example,
-        # psdem/PS_183m.mat
-        tmat = matfun.loadmat(t_fn)
-        print(' using matfun')
-        # reshaping because of how things are packed
-        tlon_vec = tmat['lon'][:].flatten()
-        tlat_vec = tmat['lat'][:].flatten()
-        tz = tmat['z'][:]
-    return tlon_vec, tlat_vec, tz
 
 def load_bathy2(t_fn, lon_vec, lat_vec):
     # load a section of the new NetCDF Smith-Sandwell bathy
-    ds = nc.Dataset(t_fn)
-    Lon = ds['lon'][:]
-    Lat = ds['lat'][:]
+    ds = xr.open_dataset(t_fn)
+    tlon_vec = ds['lon'].values
+    tlat_vec = ds['lat'].values
     i0 = zfun.find_nearest_ind(Lon, lon_vec[0])
     i1 = zfun.find_nearest_ind(Lon, lon_vec[-1])
     j0 = zfun.find_nearest_ind(Lat, lat_vec[0])
@@ -125,14 +82,10 @@ def load_bathy2(t_fn, lon_vec, lat_vec):
     return tlon_vec, tlat_vec, tz
 
 def make_nc(out_fn, plon, plat, lon, lat, z, m, dch):
-    # get rid of old version
-    try:
-        os.remove(out_fn)
-    except OSError:
-        pass # assume error was because the file did not exist
 
+    out_fn.unlink(missing_ok=True)
     # create new NetCDF file
-    foo = nc.Dataset(out_fn, 'w', format='NETCDF3_CLASSIC')
+    foo = xr.Dataset(out_fn)
 
     # create dimensions
     M, L = lon.shape # use ROMS teminology