Add new file

Spot_polynyas.py

+# This script takes Sea Ice data, and spot polynyas with a fixed threshold
+# Authors: C. Pelletier edited by N. Pirlet
+# Cite this paper for credits: 
+
+import numpy as np
+import netCDF4 as nc
+import os
+import sys
+import time
+
+# Adapt------------
+year_start = 2001
+year_stop  = 2017
+
+siconc_thres = 1.00   
+sithic_thres = 0.25    
+
+lfi = False           # If you want polynyas close to lfi to be coastal polynyas (True) (NO_LFI => False)
+
+is_cyclic = True
+if(is_cyclic):
+    n_ovlp = 1        # West-East bondary overlap
+else:
+    n_ovlp = 0
+
+
+# Contient and area files
+input_geometry       = "/cofast/npirlet/eANT_025/eANT025_mesh_mask_v2.nc" 
+input_cell_area      = "/cofast/npirlet/eANT_025/eANT025_cell_area.nc"
+
+#--------------------
+
+fillv = -10
+open_oce_value = 0
+iced_value = 1
+antarctica_pol_value = 2
+island_pol_value = 3
+ocean_pol_value = 4
+
+
+
+fillv_float = -1.e20
+
+
+
+geom = nc.Dataset(input_geometry, mode='r')
+if(is_cyclic):
+    # Creation of "land" mask 
+    tmaskutil = geom.variables['tmaskutil'][0,0:370,:-n_ovlp] == 0
+    isf_draft = geom.variables['isf_draft'][0,0:370,:-n_ovlp] > 0
+    msk_cont = np.full(np.shape(tmaskutil),False)
+    msk_cont[tmaskutil == True] = True
+    msk_cont[isf_draft == True] = True
+    is_cont = msk_cont
+    is_ant = msk_cont
+    print(np.shape(msk_cont))
+
+    if lfi == True:
+        mask_LFI = np.load("/cofast/npirlet/eANT_025/mask_LFI_MAMJJASO_2001-2017.npy")
+        mask_LFI = mask_LFI[:,:-n_ovlp]
+        msk_cont[mask_LFI >= 50] = True
+else:
+    # Creation of "land" mask
+    tmaskutil = geom.variables['tmaskutil'][0,0:370,:] == 0
+    isf_draft = geom.variables['isf_draft'][0,0:370,:] > 0
+    msk_cont = np.full(np.shape(tmaskutil),False)
+    msk_cont[tmaskutil == True] = True
+    msk_cont[isf_draft == True] = True  
+    is_cont = msk_cont
+    is_ant = msk_cont
+    print(np.shape(msk_cont))
+
+    if lfi == True:
+        mask_LFI = np.load("/cofast/npirlet/eANT_025/mask_LFI_MAMJJASO_2001-2017.npy")
+        msk_cont[mask_LFI >= 50] = True
+
+lat = geom.variables['nav_lat'][0:370,:]
+lon = geom.variables['nav_lon'][0:370,:]
+geom.close()
+
+areac = nc.Dataset(input_cell_area, mode='r')
+if(is_cyclic):
+    cell_area = areac.variables['areacello'][0,0:370,:-n_ovlp]
+    print(np.shape(cell_area))
+else:
+    cell_area = areac.variables['areacello'][0,0:370,:]
+    print(np.shape(cell_area))
+areac.close()
+
+ny, nx = np.shape(cell_area)
+
+
+is_border = np.zeros(shape=[ny, nx], dtype=bool)
+
+is_border[:1,:] = True                
+is_border[-1:,:] = True
+if(not is_cyclic):
+    is_border[:,:1]  = True           
+    is_border[:,-1:] = True
+
+
+
+is_cont_neigh = np.ndarray(shape=[ny, nx, 4], dtype=bool)
+is_cont_neigh[:,:,0] = np.roll(is_cont, shift=+1, axis=-2)
+is_cont_neigh[:,:,1] = np.roll(is_cont, shift=-1, axis=-2)
+is_cont_neigh[:,:,2] = np.roll(is_cont, shift=+1, axis=-1)
+is_cont_neigh[:,:,3] = np.roll(is_cont, shift=-1, axis=-1)
+
+is_ant_neigh = np.ndarray(shape=[ny, nx, 4], dtype=bool)
+is_ant_neigh[:,:,0] = np.roll(is_ant, shift=+1, axis=-2)
+is_ant_neigh[:,:,1] = np.roll(is_ant, shift=-1, axis=-2)
+is_ant_neigh[:,:,2] = np.roll(is_ant, shift=+1, axis=-1)
+is_ant_neigh[:,:,3] = np.roll(is_ant, shift=-1, axis=-1)
+
+
+
+old_age_polynya = np.zeros(shape=[ny, nx], dtype=int)
+
+
+for year in range(year_start,year_stop+1):
+
+    # Data path to adapt
+    input_seaice_conc = "/cofast/npirlet/LUCIA/NO_LFI_40/NO_LFI_40_1d_siconc_y{}.nc".format(year)  
+    input_seaice_thick = "/cofast/npirlet/LUCIA/NO_LFI_40/NO_LFI_40_1d_sivolu_y{}.nc".format(year)  
+    outfile      = "/cofast/npirlet/LUCIA/NO_LFI_40/polynyas/polynyas_NO_LFI_40_y{}_30cm.nc".format(year)
+    #-------------------
+    
+    # Load sea ice concentration data
+    tmpdat = nc.Dataset(input_seaice_conc, mode='r')
+    
+    if(is_cyclic):
+        siconc = tmpdat.variables['siconc'][:,0:370,:-n_ovlp]
+    else:
+        siconc = tmpdat.variables['siconc'][:,0:370,:]
+  
+    tmpdat.close()
+    
+
+    # Load sea ice thickness data
+    tmpdat = nc.Dataset(input_seaice_thick, mode='r')
+    
+    if(is_cyclic):
+        sithic = tmpdat.variables['sivolu'][:,0:370,:-n_ovlp]
+    else:
+        sithic = tmpdat.variables['sivolu'][:,0:370,:]
+
+    time_var = tmpdat.variables['time_counter'][:]
+    long_name_time = tmpdat.variables['time_counter'].standard_name
+    units_time = tmpdat.variables['time_counter'].units
+    calendar_time = tmpdat.variables['time_counter'].calendar
+    tmpdat.close()
+
+
+    # Begin process
+    
+    nt = (np.shape(sithic))[0]
+
+    idx_open_patch = fillv * np.ones(shape=[nt, ny, nx], dtype=int)
+    
+    area_polynya = np.zeros(shape=[nt, ny, nx], dtype=float)
+    
+    age_polynya = np.zeros(shape=[nt, ny, nx], dtype=int)
+    
+    idx_polynya = np.zeros(shape=[nt,ny,nx], dtype=int)
+
+
+
+    
+    for t in range(0,nt):
+
+        cnt_pol = 0
+        cnt_oce_pol = 0
+        cnt_island_pol = 0
+        cnt_ant_pol = 0
+
+        unsorted_idx_polynya = np.zeros(shape=[ny, nx], dtype=int)
+        all_poly_size = []
+        
+        t_start = time.time()
+
+        is_open = np.logical_and(np.logical_not(is_cont), siconc[t,:,:] < siconc_thres, sithic[t,:,:] < sithic_thres)
+        idx_open_patch[t,:,:] = np.where(np.logical_and(np.logical_not(is_open), np.logical_not(is_cont)), iced_value, idx_open_patch[t,:,:])
+
+
+
+        is_open_neigh = np.ndarray(shape=[ny,nx,4], dtype=bool)
+
+        is_open_neigh[:,:,0] = np.roll(is_open, shift=+1, axis=-2)
+        is_open_neigh[:,:,1] = np.roll(is_open, shift=-1, axis=-2)
+        is_open_neigh[:,:,2] = np.roll(is_open, shift=+1, axis=-1)
+        is_open_neigh[:,:,3] = np.roll(is_open, shift=-1, axis=-1)
+
+
+        explored = np.zeros(shape=[ny, nx], dtype=bool)
+        candidates = np.transpose(np.nonzero(np.logical_and(is_open, np.logical_not(explored))))
+
+
+        curr_patch_idx = 0
+        
+        sizes_patch = []
+        touched_cont = []
+
+
+        while(np.size(candidates) > 0):
+
+            curr_patch_idx+=1
+
+            is_in_currpatch = np.zeros(shape=[ny, nx], dtype=bool)
+            is_in_currpatch[candidates[0,0], candidates[0,1]] = True
+
+            keep_going = True
+            cnt_explore = 0
+
+            while(keep_going):
+
+                n_cp = np.count_nonzero(is_in_currpatch)
+
+                # add x+ neighbors
+                is_in_currpatch = np.logical_or(is_in_currpatch, np.roll(np.logical_and(is_open_neigh[:,:,0], is_in_currpatch), shift=-1, axis=0))
+
+                # add x- neighbors
+                is_in_currpatch = np.logical_or(is_in_currpatch, np.roll(np.logical_and(is_open_neigh[:,:,1], is_in_currpatch), shift=1, axis=0))
+
+                # add y+ neighbors
+                is_in_currpatch = np.logical_or(is_in_currpatch, np.roll(np.logical_and(is_open_neigh[:,:,2], is_in_currpatch), shift=-1, axis=1))
+
+                # add y- neighbors
+                is_in_currpatch = np.logical_or(is_in_currpatch, np.roll(np.logical_and(is_open_neigh[:,:,3], is_in_currpatch), shift=1, axis=1))
+                
+
+                n_old = n_cp
+                n_cp = np.count_nonzero(is_in_currpatch)
+
+                n_new = n_cp - n_old
+                keep_going = (n_new > 0)
+                cnt_explore+=1
+
+
+            if( np.any( np.logical_and(is_in_currpatch, is_border))):
+                # patch touching border: open ocean
+
+                idx_open_patch[t,:,:] = np.where( is_in_currpatch, open_oce_value, idx_open_patch[t,:,:] )
+
+            else:
+                # patch doesn't touch border: polynya
+                
+                cnt_pol+=1
+                
+                curr_area_pol = np.sum(np.ma.array(data = cell_area, mask = np.logical_not(is_in_currpatch)))
+
+                all_poly_size.append(curr_area_pol)
+                
+
+                area_polynya[t,:,:] = np.where(is_in_currpatch, curr_area_pol, area_polynya[t,:,:])
+                
+                unsorted_idx_polynya = np.where(is_in_currpatch, cnt_pol, unsorted_idx_polynya)
+            
+
+                # polynya age
+                
+                if(np.any( np.logical_and( old_age_polynya > 0, is_in_currpatch ) ) ):
+                    # polynya previously existed: its age is the maximum age of any of the current polynya node taken from the old time step
+                    prev_age_tmp = np.ma.array(data = old_age_polynya, mask = np.logical_not(is_in_currpatch))
+
+                    age_polynya[t,:,:] = np.where(is_in_currpatch, np.ma.amax(prev_age_tmp) + 1, age_polynya[t,:,:])
+
+                else:
+                    # new polynya: aged 1
+                    age_polynya[t,:,:] = np.where(is_in_currpatch, 1, age_polynya[t,:,:])
+                
+
+
+                tmp_mask = is_in_currpatch[:,:,None] * (np.ones(shape=[4], dtype=bool))[None,None,:]
+                check_cont = np.logical_and(is_cont_neigh, tmp_mask)
+                    
+                if(np.any(check_cont)):
+                    # coastal polynya: Antarctica or island
+
+
+                    check_ant = np.logical_and(is_ant_neigh, tmp_mask)
+                    if(np.any(check_ant)):
+                        # polynya touches Antarctica
+                        idx_open_patch[t,:,:] = np.where(is_in_currpatch, antarctica_pol_value, idx_open_patch[t,:,:])
+
+
+                        cnt_ant_pol+=1
+                        
+                    else:
+                        # insular (non-Antarctica) polynya
+                        idx_open_patch[t,:,:] = np.where(is_in_currpatch, island_pol_value, idx_open_patch[t,:,:])    
+
+
+                        cnt_island_pol+=1
+
+
+                else: # np.any(check_cont))
+                    # the polynya doesn't touch masked cell: it's an ocean polynya
+
+                    idx_open_patch[t,:,:] = np.where(is_in_currpatch, ocean_pol_value, idx_open_patch[t,:,:])    
+
+                    cnt_oce_pol+=1
+
+
+            sizes_patch.append(n_cp)
+
+            # we have explored the current patch
+            explored = np.where(is_in_currpatch, True, explored)
+
+            # the remaining polynya candidates are the ice-free cells which have not been explored
+            candidates = np.transpose(np.nonzero(np.logical_and(is_open, np.logical_not(explored))))
+
+        # age: fillvalue over continent, 0 over free-ocean or sea-ice cover
+        age_polynya[t,:,:] = np.where(is_cont, fillv, np.where( np.logical_or(idx_open_patch[t,:,:] == open_oce_value, idx_open_patch[t,:,:] == iced_value),0, age_polynya[t,:,:]))
+        
+        # area polynya: fillvalue over continent
+        area_polynya[t,:,:] = np.where(is_cont, fillv_float, area_polynya[t,:,:])
+
+        # sort the polynya index by decreasing polynya extent
+        idx_polynya[t,:,:] = unsorted_idx_polynya
+        argsort = (np.argsort(all_poly_size))[::-1]
+        
+        for i in range(0,cnt_pol):
+            idx_polynya[t,:,:] = np.where(unsorted_idx_polynya == argsort[i]+1, i + 1, idx_polynya[t,:,:])
+
+        idx_polynya[t,:,:] = np.where(is_cont, fillv, idx_polynya[t,:,:])
+        
+        walltime_ts = time.time() - t_start
+        print('Time step %d'%t+': spotted %d'%cnt_pol+' polynya: %d'%cnt_ant_pol+' attached to Antarctica, %d'%cnt_island_pol+' attached to an island and %d'%cnt_oce_pol+' oceanic (walltime %.1f'%walltime_ts+' sec.)')
+
+
+
+    out_nc = nc.Dataset(outfile, mode='w')
+
+    out_nc.createDimension('time_counter', 0)
+    out_nc.createDimension('axis_nbounds', 2)
+    out_nc.createDimension('y', ny)
+    out_nc.createDimension('x', nx+n_ovlp)
+
+
+    tmp = out_nc.createVariable(varname="time_counter", datatype='f8', dimensions=['time_counter'])
+    tmp.standard_name = long_name_time
+    tmp.long_name     = long_name_time
+    tmp.units         = units_time
+    tmp.calendar      = calendar_time
+    tmp[:] = time_var[:]
+
+    tmp = out_nc.createVariable(varname="latitude", datatype='f8', dimensions=['y', 'x'])
+    tmp.standard_name = "latitude"
+    tmp.units = "degrees North"
+    tmp[:] = lat
+
+    tmp = out_nc.createVariable(varname="longitude", datatype='f8', dimensions=['y', 'x'])
+    tmp.standard_name = "longitude"
+    tmp.units = "degrees East"
+    tmp[:] = lon
+
+    
+    tmp = out_nc.createVariable(varname='polynya_status', datatype='f8', dimensions=['time_counter', 'y', 'x'], fill_value=fillv)
+    tmp.coordinates = 'time_counter latitude longitude'
+    tmp.description = "polnynya status. masked = continent; %d"%open_oce_value+' = open ocean; %d'%iced_value+" = sea-ice cover; %d"%antarctica_pol_value+" = polynya attached to Antarctica; %d"%island_pol_value+" = polynya attached to an island; %d"%ocean_pol_value+" = oceanic polynya"
+    if(is_cyclic):
+        tmp[:] = np.dstack((idx_open_patch, idx_open_patch[:,:,:n_ovlp]))
+    else:
+        tmp[:] = idx_open_patch
+
+    
+    tmp = out_nc.createVariable(varname='area_polynya', datatype='f4', dimensions=['time_counter', 'y', 'x'], fill_value =  fillv_float)
+    if(is_cyclic):
+        tmp[:] = np.dstack((area_polynya, area_polynya[:,:,:n_ovlp]))
+    else:
+        tmp[:] = area_polynya
+    tmp.coordinates = 'time_counter latitude longitude'
+    tmp.units = 'm2'
+    
+    
+    tmp = out_nc.createVariable(varname='age_polynya', datatype='i2', dimensions=['time_counter', 'y', 'x'], fill_value =  fillv)
+    tmp.coordinates = 'time_counter latitude longitude'
+    tmp.units = 'time step (days)'
+    if(is_cyclic):
+        tmp[:] = np.dstack((age_polynya, age_polynya[:,:,:n_ovlp]))
+    else:
+        tmp[:] = age_polynya
+    
+
+    tmp = out_nc.createVariable(varname='idx_polynya', datatype='i2', dimensions=['time_counter', 'y', 'x'], fill_value =  fillv)
+    tmp.coordinates = 'time_counter latitude longitude'
+    tmp.description = 'index counting polynyas (in decreasing order of extent)'
+    if(is_cyclic):
+        tmp[:] = np.dstack((idx_polynya, idx_polynya[:,:,:n_ovlp]))
+    else:
+        tmp[:] = idx_polynya
+     
+    out_nc.close()