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

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

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():

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

{'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: 130
Joined: Mon Mar 30, 2015 11:22 am