## Marine and Coastal Phenology Index Calculation

### Marine and Coastal Phenology Index Calculation

Ocean phenology indices have been an emerging topic and closely linked as a environmental application of satellite derived essential climate variables (ECVs) to describe the ocean variability in a changing world. Typical metrics include trend analysis of SST, start and end days of summer, peak of spring transition etc., with predefined metric thresholds or criteria. We have developed an open source python software module (https://github.jpl.nasa.gov/ybj/PhenologyPy) to calculate satellite-based marine phenology products for the National Climate Assessment. Primarily prototyped on Sea Surface Temperature (SST) measurements from a number of Level 4 datasets over the North Atlantic Ocean it has been extended to handle any region of the world’s oceans. Some explorations have been made with long time series ocean color and other datasets in addition to SST.

In this recipe, we will show you the code snippet for the summer start index calculation from the project. The routine is called metric_summer_start_end in the phenologyalg.py module.

Code: Select all
`def metric_summer_start_end(lats, lons, startY, endY, summer_thresh_offset, variablename, outdir, outfilename, clim_filename, dirname):  ny = endY - startY + 1  # define max and min matrics  day_min = np.empty((lats.size, lons.size, ny))  day_max = np.empty((lats.size, lons.size, ny))  data_min = np.empty((lats.size, lons.size, ny))  data_max = np.empty((lats.size, lons.size, ny))  # define summer start/end matrics  day_summer_start = np.empty((lats.size, lons.size, ny))  day_summer_end = np.empty((lats.size, lons.size, ny))  data_summer_max_all = np.empty((lats.size, lons.size, ny))  data_summer_max = np.empty((lats.size, lons.size))  day_summer_max = np.empty((lats.size, lons.size))  sst_all        = np.empty((lats.size, lons.size, 366))  sst_all_number = np.empty((lats.size, lons.size, 366))  data_year_all  = np.empty((ny, 12, 31, lats.size, lons.size))  data_year      = np.empty((lats.size, lons.size, 366, ny))  # define summer start/end matrics trend  day_summer_start_trend = np.empty((lats.size, lons.size))  day_summer_end_trend = np.empty((lats.size, lons.size))  # initilization with NAN  data_year_all[:] = np.nan  data_year[:] = np.nan  day_summer_start[:] = np.nan  day_summer_end[:] = np.nan  data_summer_max[:] = np.nan  data_summer_max_all[:] = np.nan  day_summer_max[:] = np.nan  day_summer_start_trend[:] = np.nan  day_summer_end_trend[:] = np.nan  nlat = lats.size  nlon = lons.size  # Read in climatology  ncin = Dataset(clim_filename, 'r')  diff_clim = ncin.variables['diff_sst_climatology'][:]  ncin.close()  day_turn = np.empty([lats.size, lons.size])  day_turn[:,:] = np.nan  for i in range(0, lats.size):    for j in range(0, lons.size):      xrange = np.arange(366)      yrange = diff_clim[i,j,:].squeeze()      idx = np.argwhere(np.isfinite(yrange))      if idx.size > len(xrange)/2:        z = np.polyfit(xrange[idx[:,0].squeeze()], yrange[idx[:,0].squeeze()], 6)        polynomial = np.poly1d(z)        ys = polynomial(range(366))        for k in range(0, 365):          if(ys[k] > 0):            day_turn[i, j] = k + 1            break  for i in range(startY, endY+1):    print("Processing Year = " + str(i))    nd = 0    for j in range(1, 367):       for filename in glob.glob(dirname+'/'+str(i)+'/'+"{0:0>3}".format(j)+'/*.nc'):         ncfile = filename.rsplit( "/")[ -1 ]         ncin = Dataset(filename, 'r')         asst = ncin.variables[variablename][:]         sst = asst.squeeze()         #sst = asst[0,:,:]         lons = ncin.variables['lon'][:]         lats = ncin.variables['lat'][:]         amask = ncin.variables['mask'][:]         mask = amask[0,:,:]         #get global attributes         atts = ncin.__dict__         for att,val in atts.items():           if att.find('start_time') != -1:             start_time = val             continue         ncin.close()         aindex = np.where( (sst > -250.0) & (sst < 350.0) & (mask == 1) )         data_year[aindex[0][:],aindex[1][:],j-1,i-startY] = sst[aindex[0][:],aindex[1][:]]    # data max/date and min/date calculation    for i1 in range(0, lats.size):      for j1 in range(0, lons.size):           atemp = data_year[i1,j1,:,i-startY]           atemp1 = atemp.squeeze()           mx = np.ma.masked_array(atemp1,np.isnan(atemp1))           atemp = moving_average(mx)           if np.count_nonzero(np.isnan(atemp)) == atemp.size:             day_min[i1,j1,i-startY]  = np.nan             data_min[i1,j1,i-startY] = np.nan             day_max[i1,j1,i-startY]  = np.nan             data_max[i1,j1,i-startY] = np.nan           else:             day_min[i1,j1,i-startY]  = np.nanargmin(atemp)             if day_min[i1,j1,i-startY] > 300:               day_min[i1,j1,i-startY] = np.nan             data_min[i1,j1,i-startY] = np.nanmin(atemp)             day_max[i1,j1,i-startY]  = np.nanargmax(atemp)             data_max[i1,j1,i-startY] = np.nanmax(atemp)    # summer start/end matrics calculation    # Step 1: find summer maximum for each year     for i1 in range(0, lats.size):      for j1 in range(0, lons.size):           atemp = data_year[i1,j1,150:150+90,i-startY]           atemp = atemp.squeeze()           mx = np.ma.masked_array(atemp,np.isnan(atemp))           atemp = moving_average(mx)           if np.count_nonzero(np.isnan(atemp)) == atemp.size:             data_summer_max_all[i1,j1,i-startY]  = np.nan           else:             data_summer_max_all[i1,j1,i-startY] = np.nanmax(atemp)  # summer start/end matrics calculation  # Step 2: summer lowest maximum calculation  for i1 in range(0, lats.size):    for j1 in range(0, lons.size):       atemp = data_summer_max_all[i1,j1,:]       atemp = atemp.squeeze()       if np.count_nonzero(np.isnan(atemp)) == atemp.size:         data_summer_max[i1,j1]  = np.nan         day_summer_max[i1,j1]  = np.nan       else:         data_summer_max[i1,j1] = np.nanmin(atemp)         day_summer_max[i1,j1] = np.nanargmin(atemp)+startY  # summer start/end matrics calculation  # Step 3: summer start/end calculation  for i in range(startY, endY+1):    for i1 in range(0, lats.size):      for j1 in range(0, lons.size):           atemp = data_year[i1,j1,:,i-startY]           atemp = atemp.squeeze()           mx = np.ma.masked_array(atemp,np.isnan(atemp))           atemp = moving_average(mx)           """           if np.count_nonzero(np.isnan(atemp)) == atemp.size:             day_summer_start[i1,j1,i-startY]  = np.nan            elif np.isnan(day_turn[i1, j1]):             day_summer_start[i1,j1,i-startY]  = np.nan            else:            """           if np.count_nonzero(np.isnan(atemp)) <= atemp.size and ~np.isnan(day_min[i1,j1,i-startY]):             for k in range(int(day_turn[i1, j1]), atemp.size):               if atemp[k] >= data_summer_max[i1,j1]+summer_thresh_offset:                 day_summer_start[i1,j1,i-startY] = k                 break             for k in range(np.nanargmax(atemp), atemp.size):               if atemp[k] <= data_summer_max[i1,j1]+summer_thresh_offset:                 day_summer_end[i1,j1,i-startY] = k                 break  # Calculate the summer start/end day trend   x = np.arange(ny)  for i1 in range(0, lats.size):    for j1 in range(0, lons.size):      """      z = np.polyfit(range(ny), day_summer_start[i1,j1,:], 1)      day_summer_start_trend[i1,j1] = z[0]      z = np.polyfit(range(ny), day_summer_end[i1,j1,:], 1)      day_summer_end_trend[i1,j1] = z[0]      """      y = day_summer_start[i1,j1,:].squeeze()      idx = np.isfinite(y)      x1 = x[idx]      if x1.size > ny/2:        z = np.polyfit(x[idx], y[idx], 1)        day_summer_start_trend[i1,j1] = z[0]      else:        day_summer_start_trend[i1,j1] = np.nan      y = day_summer_end[i1,j1,:].squeeze()      idx = np.isfinite(y)      x1 = x[idx]      if x1.size > ny/2:        z = np.polyfit(x[idx], y[idx], 1)        day_summer_end_trend[i1,j1] = z[0]      else:        day_summer_end_trend[i1,j1] = np.nan  # Update the variables with _Fill_Value (-9999)  day_min[np.isnan(day_min)] = -9999  day_max[np.isnan(day_max)] = -9999  data_min[np.isnan(data_min)] = -9999  data_max[np.isnan(data_max)] = -9999  day_summer_start[np.isnan(day_summer_start)] = -9999  day_summer_end[np.isnan(day_summer_end)] = -9999  day_summer_start_trend[np.isnan(day_summer_start_trend)] = -9999  day_summer_end_trend[np.isnan(day_summer_end_trend)] = -9999  ### create output directory  data_info.createdir(outdir)  #output into netCDF file  fid = Dataset(outdir+'/'+outfilename,'w')  # Define the dimensions  nlat = fid.createDimension('lat', nlat) # Unlimited  nlon = fid.createDimension('lon', nlon) # Unlimited  nyear = fid.createDimension('year', ny) # Unlimited  nc_var = fid.createVariable('lat', 'f8',('lat'),zlib=True)  fid.variables['lat'][:] = lats  fid.variables['lat'].standard_name='latitude'  fid.variables['lat'].long_name='latitude'  fid.variables['lat'].axis='Y'  fid.variables['lat'].units='degrees_north'  nc_var = fid.createVariable('lon', 'f8',('lon'),zlib=True)  fid.variables['lon'][:] = lons  fid.variables['lon'].standard_name='longitude'  fid.variables['lon'].long_name='longitude'  fid.variables['lon'].units='degrees_east'  fid.variables['lon'].axis='X'  nc_var = fid.createVariable('year', 'i4',('year'),zlib=True)  fid.variables['year'][:] = range(startY,endY+1)  fid.variables['year'].standard_name='Year'  nc_var = fid.createVariable('day_min', 'i4',('lat','lon','year'),fill_value=-9999,zlib=True)  fid.variables['day_min'][:] = day_min  fid.variables['day_min'].standard_name='coldest day of the year'  fid.variables['day_min'].units='Day number in a year'  nc_var = fid.createVariable('day_max', 'i4',('lat','lon','year'),fill_value=-9999,zlib=True)  fid.variables['day_max'][:] = day_max  fid.variables['day_max'].standard_name='warmest day of the year'  fid.variables['day_max'].units='Day number in a year'  nc_var = fid.createVariable('data_min', 'f8',('lat','lon','year'),fill_value=-9999,zlib=True)  fid.variables['data_min'][:] = data_min  fid.variables['data_min'].standard_name='coldest data of the year'  fid.variables['data_min'].units=''  nc_var = fid.createVariable('data_max', 'f8',('lat','lon','year'),fill_value=-9999,zlib=True)  fid.variables['data_max'][:] = data_max  fid.variables['data_max'].standard_name='warmest data of the year'  fid.variables['data_max'].units=''  nc_var = fid.createVariable('day_summer_start', 'i4',('lat','lon','year'),fill_value=-9999,zlib=True)  fid.variables['day_summer_start'][:] = day_summer_start  fid.variables['day_summer_start'].standard_name='summer start day of the year'  fid.variables['day_summer_start'].units='Day number in a year'  nc_var = fid.createVariable('day_summer_end', 'i4',('lat','lon','year'),fill_value=-9999,zlib=True)  fid.variables['day_summer_end'][:] = day_summer_end  fid.variables['day_summer_end'].standard_name='summer start day of the year'  fid.variables['day_summer_end'].units='Day number in a year'  nc_var = fid.createVariable('day_summer_start_trend', 'f8',('lat','lon'),fill_value=-9999,zlib=True)  fid.variables['day_summer_start_trend'][:] = day_summer_start_trend  fid.variables['day_summer_start_trend'].standard_name='Summer End Day Trend per Year'  fid.variables['day_summer_start_trend'].units='Day number per year'  nc_var = fid.createVariable('day_summer_end_trend', 'f8',('lat','lon'),fill_value=-9999,zlib=True)  fid.variables['day_summer_end_trend'][:] = day_summer_end_trend  fid.variables['day_summer_end_trend'].standard_name='Summer Start Day Trend per Year'  fid.variables['day_summer_end_trend'].units='Day number per Year'  fid.cdm_data_type         = "Grid"  fid.close()  return`

Reference

Mills, Katherine E., Andrew J. Pershing, and Christina M. Hernández. "Forecasting the seasonal timing of Maine's lobster fishery." Frontiers in Marine Science 4 (2017): 337.
Wang, Lifei, Lisa A. Kerr, Nicholas R. Record, Eric Bridger, Benjamin Tupper, Katherine E. Mills, Edward M. Armstrong, and Andrew J. Pershing. "Modeling marine pelagic fish species spatiotemporal distributions utilizing a maximum entropy approach." Fisheries Oceanography 27, no. 6 (2018): 571-586.
Thomas, Andrew C., Andrew J. Pershing, Kevin D. Friedland, Janet A. Nye, Katherine E. Mills, Michael A. Alexander, Nicholas R. Record, Ryan Weatherbee, and M. Elisabeth Henderson. "Seasonal trends and phenology shifts in sea surface temperature on the North American northeastern continental shelf." Elem Sci Anth 5 (2017).
yiboj

Posts: 130
Joined: Mon Mar 30, 2015 11:22 am