PODAAC

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).