Add zone specific deltas computation

Plus merge land computation in compute_clip

Compute/compute_deltas_daily.py

@@ -15,6 +15,7 @@ from tomate import db_types as dt
 import tomate.scan_library as scanlib
 
 from lib import root_data
+from lib.zones import get_zones
 import lib.data.mask
 import lib.data.sst
 import lib.data.chl
@@ -22,7 +23,9 @@ import lib.data.chl
 
 step = 8
 scale = 10
-level = 4
+level = 3
+HI_number = 6
+zones = ['1', '2']
 box_lat = [20, 52]
 box_lon = [-82, -45]
 range_time = [[2007, 1, 1], [2007, 12, 31]]
@@ -40,21 +43,27 @@ cstr_lo.set_scan_filename(lambda *args, **kwargs:
                           'time', default_date={'hour': 9})
 db_lo = cstr_lo.make_data()
 
+dbz_lo = get_zones('lo')
+dbz_hi = get_zones('hi')
+
 for db in [db_lo, db_hi]:
     db.select_by_value('avail',
                        lat=slice(*box_lat), lon=slice(*box_lon),
                        time=slice(*range_time), by_day=True)
 
-
-## Compute
+for db, dbz in zip([db_lo, db_hi], [dbz_lo, dbz_hi]):
+    dbz.load_by_value(var=zones,
+                      lat=slice(*db.selected.lat.get_extent()),
+                      lon=slice(*db.selected.lon.get_extent()))
 
 
+## Compute
 def get_land_mask_sst():
     db_hi.load_selected(var='SST_mask', time=0)
     land = db_hi['SST_mask'][0] != 1
     mask = 1.*land.copy().filled(False)
 
-    remove_bahamas(db_hi, mask)
+    remove_bermudes(db_hi, mask)
     mask = extend_mask(mask, 10, 2)
     return mask > .5
 
@@ -64,15 +73,15 @@ def get_land_mask_chl():
     land = db_lo['CHL_flags'][0] == 9
     mask = 1.*land.copy().filled(False)
 
-    remove_bahamas(db_lo, mask)
+    remove_bermudes(db_lo, mask)
     mask = extend_mask(mask, 5, 1)
     return mask > .5
 
 
-def remove_bahamas(db, mask):
-    bahamas_lat = db.loaded.lat.subset(32.20, 32.45)
-    bahamas_lon = db.loaded.lon.subset(-64.95, -64.60)
-    mask[bahamas_lat, bahamas_lon] = False
+def remove_bermudes(db, mask):
+    bermudes_lat = db.loaded.lat.subset(32.20, 32.45)
+    bermudes_lon = db.loaded.lon.subset(-64.95, -64.60)
+    mask[bermudes_lat, bermudes_lon] = False
 
 
 def extend_mask(mask, neighbors, repeat):
@@ -102,8 +111,8 @@ def get_stats_init(db):
     coords_fg = [[var, 'in'], [time, 'in']]
     db_stats = dt.DataDisk([var, time],
                            root=path.join(root_data,
-                                          'HI/HI_{:.1f}/L{:d}_step{:d}'
-                                          .format(scale, level, step)),
+                                          'HI/HI_{:.1f}_{:d}/L{:d}_step{:d}'
+                                          .format(scale, HI_number, level, step)),
                            filegroups=[])
     db_stats.add_filegroup(FilegroupNetCDF, coords_fg, 'stats')
     db_stats.filegroups[-1].cs['var'].update_values(var[:], var[:])
@@ -121,60 +130,42 @@ def find_stats(values, stats, i):
     stats['q90'][i] = np.percentile(values, 90)
 
 
-def compute_land(db, var, land):
-    stats_ft = get_stats_init(db)
-    stats_bg = get_stats_init(db)
-
-    for i, time_slice in enumerate(db.selected.iter_slices('time', step)):
-        db.load_selected(var='front', time=time_slice)
-        ft = db['front'].astype('bool').filled(False) * ~land
-        db.load_selected(var='background', time=time_slice)
-        bg = db['background'].astype('bool').filled(False) * ~land
-
-        db.load_selected(var=var, time=time_slice)
-        values_ft = get_values(db[var], ft)
-        values_bg = get_values(db[var], bg)
-        db.unload_data()
-
-        find_stats(values_ft, stats_ft, i)
-        find_stats(values_bg, stats_bg, i)
-
-    return stats_ft, stats_bg
+def compute_clip(db, var, dbz, land):
 
+    for zone in dbz.loaded.var[:].tolist() + ['total']:
+        stats_ft = get_stats_init(db)
+        stats_bg = get_stats_init(db)
 
-def compute_clip(db, var):
-    stats_ft = get_stats_init(db)
-    stats_bg = get_stats_init(db)
+        for i, time_slice in enumerate(db.selected.iter_slices('time', step)):
+            db.load_selected(var='front', time=time_slice)
+            ft = db['front'].astype('bool').filled(False)
+            db.load_selected(var='background', time=time_slice)
+            bg = db['background'].astype('bool').filled(False)
 
-    for i, time_slice in enumerate(db.selected.iter_slices('time', step)):
-        db.load_selected(var='front', time=time_slice)
-        ft = db['front'].astype('bool').filled(False)
-        db.load_selected(var='background', time=time_slice)
-        bg = db['background'].astype('bool').filled(False)
+            if zone != 'total':
+                ft *= dbz[zone]
+                bg *= dbz[zone]
 
-        db.load_selected(var=var, time=time_slice)
-        values_ft = get_values(db[var], ft)
-        values_bg = get_values(db[var], bg)
-        db.unload_data()
+            ft *= ~land
+            bg *= ~land
 
-        if var == 'CHL':
-            values_ft = np.delete(values_ft, np.where(values_ft > 3.))
-            values_bg = np.delete(values_bg, np.where(values_bg > 3.))
+            db.load_selected(var=var, time=time_slice)
+            values_ft = get_values(db[var], ft)
+            values_bg = get_values(db[var], bg)
+            db.unload_data()
 
-        find_stats(values_ft, stats_ft, i)
-        find_stats(values_bg, stats_bg, i)
+            if var == 'CHL':
+                values_ft = np.delete(values_ft, np.where(values_ft > 3.))
+                values_bg = np.delete(values_bg, np.where(values_bg > 3.))
 
-    return stats_ft, stats_bg
+            find_stats(values_ft, stats_ft, i)
+            find_stats(values_bg, stats_bg, i)
 
+        stats_ft.write(path.join(f'zone_{zone}', f'stats_{var}_ft.nc'))
+        stats_bg.write(path.join(f'zone_{zone}', f'stats_{var}_bg.nc'))
 
-# land_mask = get_land_mask_sst()
-# stats_hi_ft, stats_hi_bg = compute_land(db_hi, 'SST', land_mask)
-stats_hi_ft, stats_hi_bg = compute_clip(db_hi, 'SST')
-stats_hi_ft.write('stats_sst_ft.nc')
-stats_hi_bg.write('stats_sst_bg.nc')
 
-# land_mask = get_land_mask_chl()
-# stats_lo_ft, stats_lo_bg = compute_land(db_lo, 'CHL', land_mask)
-stats_lo_ft, stats_lo_bg = compute_clip(db_lo, 'CHL')
-stats_lo_ft.write('stats_chl_ft.nc')
-stats_lo_bg.write('stats_chl_bg.nc')
+land_mask = get_land_mask_sst()
+compute_clip(db_hi, 'SST', dbz_hi, land_mask)
+land_mask = get_land_mask_chl()
+compute_clip(db_lo, 'CHL', dbz_lo, land_mask)

README.md

@@ -18,7 +18,8 @@ Collocate submesoscale fronts to phytoplankton levels using satellite imagery
     matplotlib
     cftime>=1.1.3
     netcdf4
+    shapely
 
 All loading of data from disk (and most of writing results to disk) is done
-using [tomate](https://github.com/Descanonge/tomate), a custom made python
-package to load data from multiple files.
\ No newline at end of file
+using [tomate (v1.\*)](https://github.com/Descanonge/tomate), a custom made python
+package to load data from multiple files.

lib/data/sst.py

@@ -41,13 +41,13 @@ def add_sst(cstr):
     cstr.add_post_loading_func(mask_sst_land, ['SST', 'SST_mask'],
                                all_variables=True, current_fg=True)
 
+
 def k2c(db):
     """Change SST from Kelvin to Celsius."""
     i = db.loaded.idx('SST')
     db.data[i] -= 273.15
 
 
-
 def mask_sst_error(db):
     mask = (db['SST'].mask + db['SST_error'].data >= .40)
     db.set_mask('SST', mask)

lib/zones.py

 """Define, create, or load zones."""
 
 from os import path
-from glob import glob
 
 import numpy as np
 
 from tomate import Lat, Lon, Constructor
 from tomate.filegroup import FilegroupNetCDF
 import tomate.scan_library as scanlib
+import tomate.db_types as dt
 
 from lib import root_data
 
@@ -19,30 +19,25 @@ import lib.data.chl
 import lib.data.hi
 
 
-
 def make_zones():
-    lon_min, lon_max, lat_min, lat_max = -84., -40., 18., 62.
-
-    gs = np.array([[-75., 34.], [lon_max, 40.],
-                   [lon_max, 60], [-75., 40.]])
-
-    gy = np.array([[lon_min, lat_min], [lon_min, 34.], [-75., 34.],
-                   [lon_max, 40.], [lon_max, lat_min]])
+    lon_min, lon_max, lat_min, lat_max = -84., -40., -18., 62.
 
-    gs2 = np.array([[-75., 32.], [lon_max, 38.],
-                    [lon_max, 60], [-75., 40.]])
+    zones = {}
+    # 1 - GS
+    zones['1'] = np.array([[-75., 32.], [lon_max, 38.],
+                           [lon_max, 60], [-75., 38.]])
+    # 2 - Gyre, Full
+    zones['2'] = np.array([[lon_min, lat_min], [lon_min, 32.], [-75., 32.],
+                           [lon_max, 38.], [lon_max, lat_min]])
 
-    gy2 = np.array([[lon_min, lat_min], [lon_min, 32.], [-75., 32.],
-                    [lon_max, 38.], [lon_max, lat_min]])
-
-    zones = {'gulfstream': gs, 'gyre': gy, 'gulfstream_2': gs2, 'gyre_2': gy2}
+    # 3 - Gyre, removing coast and Bahamas
     return zones
 
 
 def rasterize_zones(db, zones):
-    lat = db.selected.lat.copy()
-    lon = db.selected.lon.copy()
-    cstr = Constructor(root_data, [lat, lon])
+    lat = db.avail.lat.copy()
+    lon = db.avail.lon.copy()
+    cstr = Constructor(path.join(root_data, 'zones'), [lat, lon])
     cstr.add_filegroup(FilegroupNetCDF, [['lat', 'in'], ['lon', 'in']],
                        'ZONES')
     cstr.set_variables_infile(**{name: name for name in zones})
@@ -63,11 +58,13 @@ def get_zones(kind='hi'):
     lat = Lat()
     lon = Lon()
     cstr = Constructor(root_data, [lat, lon])
-    cstr.add_filegroup(FilegroupNetCDF, [['lat', 'in'], ['lon', 'in']], 'ZONES')
-    cstr.set_fg_regex(rf'zones_{kind}\.nc')
+    cstr.add_filegroup(FilegroupNetCDF, [['lat', 'in'], ['lon', 'in']], 'ZONES',
+                       variables_shared=True)
+    cstr.set_fg_regex(rf'%(var:char)_{kind}\.nc')
     cstr.set_scan_in_file(scanlib.nc.scan_in_file, 'lat', 'lon')
     cstr.set_scan_in_file(scanlib.nc.scan_variables, 'var')
     cstr.add_post_loading_func(to_bool, slice(None))
+    cstr.set_data_types([dt.DataPlot])
     db = cstr.make_data()
 
     return db
@@ -79,6 +76,7 @@ if __name__ == '__main__':
     for db, kind in zip([lib.data.chl.get_data(),
                          lib.data.sst.get_data()],
                         ['lo', 'hi']):
-        db.select_by_value(lat=slice(17, 63), lon=slice(-85, -39))
         db_ = rasterize_zones(db, zones)
-        db_.write(f'zones_{kind}.nc', var_kw={'_all': {'datatype': 'i1'}})
+        for name, mask in zones.items():
+            db_.write(f'{name}_{kind}.nc', var_kw={'_all': {'datatype': 'i1'}},
+                      var=name)