w10n Level 2 Subsetting Output Repacking from 1D to 2D Array

w10n Level 2 Subsetting Output Repacking from 1D to 2D Array

Postby yiboj » Wed Dec 07, 2016 11:17 am

Webification (W10n; https://podaac.jpl.nasa.gov/podaac_labs) for science is a very useful tool to subset level 2, level 3 and level 4 data granules, especially for level 2 dataset since there are very limited tools available. But unfortunately the w10n output is in 1d array format which was lumped together from the original 2d matrix. Here we provide the code snippet to reformat the 1d array back to 2d matrix.

Code: Select all
def rePackw10nOutput(lats_1d, lons_1d, data_1d):

  """
  The algorithm uses the longitude value as the indicator for the end of each scanning line
  The longitude value difference between two ajacent points will change the sign from postive to negative
  at the end of the scanning line

  Input Parameters:
     lats_1d:  latitude input in 1d array
     lons_1d:  longitude input in 1d array
     data_1d:  data input in 1d array
  Output Parameters:
     ni:   latitude size
     nj:   longitude size
     lat:  latitude in 2d array (nj, ni)
     lon:  longitude in 2d array (nj, ni)
     data:  data in 2d array (nj, ni)
  Usage:
     ni, nj, lat, lon, data = rePackw10nOutput(lats_1d, lons_1d, data_1d)
  """

  # calculate the longitude difference between ajacent points
  lon_diff = np.diff(lons_1d)

  # find end of scan line
  scan_index = np.where( lon_diff < 0 )
  scan_index = scan_index[0]
  print scan_index

  # find the length of each scan line
  scan_length = np.diff(scan_index)
  scan_length = np.insert(scan_length, 0, scan_index[0])
  scan_length = np.append(scan_length, lons_1d.size - scan_index[-1])
  scan_index = np.append(scan_index, lons_1d.size)

  # find the size of new 2d array
  nj = max(scan_length)
  ni = scan_length.size

  # declare new 2d outout parameters
  lat = np.empty((nj, ni))
  lon = np.empty((nj, ni))
  data = np.empty((nj, ni))
  lat[:,:] = np.nan
  lon[:,:] = np.nan
  data[:,:] = np.nan

  # fill in the values of 2d parametrs using 1d parameters
  for i in range(0, ni):
    if i == 0:
      lat[0:scan_length[0], i] = lats_1d[0:scan_index[0]]
      lon[0:scan_length[0], i] = lons_1d[0:scan_index[0]]
      data[0:scan_length[0], i] = data_1d[0:scan_index[0]]
    else:
      lat[0:scan_length[i], i] = lats_1d[scan_index[i-1]:scan_index[i]]
      lon[0:scan_length[i], i] = lons_1d[scan_index[i-1]:scan_index[i]]
      data[0:scan_length[i], i] = data_1d[scan_index[i-1]:scan_index[i]]

  # return all the results
  return ni, nj, lat, lon, data


In the following code, we show the example how to use the subroutine to access the subsetted output from the L2 dataset AMSRE-REMSS-L2P-v7a from REMSS. The subsetted output is from the post "Python Dataset Subsetting Script Using w10n":

Code: Select all
from netCDF4 import Dataset
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
import numpy as np
import sys, os
import math
import pyresample as pr
from pyresample import image, geometry

###############################################################################

def rePackw10nOutput(lats_1d, lons_1d, data_1d):

  """
  The algorithm uses the longitude value as the indicator for the end of each scanning line
  The longitude value difference between two ajacent points will change the sign from postive to negative
  at the end of the scanning line

  Input Parameters:
     lats_1d:  latitude input in 1d array
     lons_1d:  longitude input in 1d array
     data_1d:  data input in 1d array
  Output Parameters:
     ni:   latitude size
     nj:   longitude size
     lat:  latitude in 2d array (nj, ni)
     lon:  longitude in 2d array (nj, ni)
     data:  data in 2d array (nj, ni)
  How to use the function:
     ni, nj, lat, lon, data = rePackw10nOutput(lats_1d, lons_1d, data_1d)
  """

  # calculate the longitude difference between ajacent points
  lon_diff = np.diff(lons_1d)

  # find end of scan line
  scan_index = np.where( lon_diff < 0 )
  scan_index = scan_index[0]
  print scan_index

  # find the length of each scan line
  scan_length = np.diff(scan_index)
  scan_length = np.insert(scan_length, 0, scan_index[0])
  scan_length = np.append(scan_length, lons_1d.size - scan_index[-1])
  scan_index = np.append(scan_index, lons_1d.size)

  # find the size of new 2d array
  nj = max(scan_length)
  ni = scan_length.size

  # declare new 2d outout parameters
  lat = np.empty((nj, ni))
  lon = np.empty((nj, ni))
  data = np.empty((nj, ni))
  lat[:,:] = np.nan
  lon[:,:] = np.nan
  data[:,:] = np.nan

  # fill in the values of 2d parametrs using 1d parameters
  for i in range(0, ni):
    if i == 0:
      lat[0:scan_length[0], i] = lats_1d[0:scan_index[0]]
      lon[0:scan_length[0], i] = lons_1d[0:scan_index[0]]
      data[0:scan_length[0], i] = data_1d[0:scan_index[0]]
    else:
      lat[0:scan_length[i], i] = lats_1d[scan_index[i-1]:scan_index[i]]
      lon[0:scan_length[i], i] = lons_1d[scan_index[i-1]:scan_index[i]]
      data[0:scan_length[i], i] = data_1d[scan_index[i-1]:scan_index[i]]

  # return all the results
  return ni, nj, lat, lon, data

def standalone_main():

  filename = '20100101020800-REMSS-L2P_GHRSST-SSTsubskin-AMSRE-l2b_v07a_r40759.dat-v02.0-fv01.0_subset.nc'
  ncin = Dataset(filename, 'r')
  lons = ncin.variables['lon'][:]
  lats = ncin.variables['lat'][:]
  sst = ncin.variables['sea_surface_temperature'][:]
  ncin.close()

  ni, nj, lat, lon, data = rePackw10nOutput(lats, lons, sst)

  fig = plt.figure()

  area_def = geometry.AreaDefinition('area1', 'Test1', 'area1',
        {'a': '6370998.0', 'units': 'm',
        'lat_0': '0.0',
        'lon_0': '0.0', 'proj': 'ortho'},
        640, 320,
        [-6370998.0, -6370998.0, 6370998.0, 6370998.0])
        #[-5326849.0625,-5326849.0625,5326849.0625,5326849.0625])
        #[-1370912.72, -909968.64, 1029087.28, 1490031.36])
  swath_def = pr.geometry.SwathDefinition(lon, lat)

  bmap = pr.plot.area_def2basemap(area_def)
  bmng = bmap.bluemarble()
  result = pr.kd_tree.resample_nearest(swath_def, data, area_def, radius_of_influence=20000, fill_value=None)
  col = bmap.imshow(result, origin='upper', vmin=250, vmax=310)

  plt.show()

###############################################################################

if __name__ == "__main__":
        standalone_main()




References:
1. Webification (W10n)
2. Python Dataset Subsetting Script Using w10n
yiboj
 
Posts: 124
Joined: Mon Mar 30, 2015 11:22 am

Return to Data Format