- 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