This Python script takes NASA LANCE lvl 2 AMSR2 swath data from the GHRC DAAC and grids the data to a GeoTIFF output

grid_amsrL2Swath.py

#*****************************************************************************
# FILE: grid_amsrL2Swath.py
# AUTHOR: Kel Markert
# EMAIL: [email protected]
# MODIFIED BY: N/A
# ORGANIZATION: UAH/ESSC
# CREATION DATE: 02/10/2017
# LAST MOD DATE: 02/11/2017
# PURPOSE: This scripts takes LANCE L2 swath AMSR2 data from GHRC, grids the
#           data, and outputs a geotiff
# DEPENDENCIES: numpy, scipy, osgeo, gdal, pyproj, pillow, pytables
#*****************************************************************************

# Import dependencies
import tables
import numpy as np
from scipy.interpolate import griddata
from osgeo import gdal
from osgeo.gdalnumeric import *  
from osgeo.gdalconst import *
from pyproj import Proj, transform
from PIL import Image, ImageDraw

def get_swath_poly(xx,yy,ns_tolerance,density=1):
    """
    FUNCTION:  poly = get_swath_vertices(xx,yy,ns_tolerance)
    ARGUMENTS: xx - longitude geolocation coordinates for swath
               yy - latitude geolocation coodinates for swath
               ns_tolerance - latitude to clip data to (typically high 60
                              degress +)
    KEYWORDS:  density - default = 1; the sampling density along latitudes for 
                         the vertices, input is in degree latitude
    RETURNS:   pts - the polygon vertices calculated from the swath
    NOTES:     Walks up/down latitudes in between the ns_tolerance at an 
               interval of density and calculates a vextex at each interval.
    """
    
    # Convert latitude density to meters
    proj_density = density * (110.6*1E3)
    
    # Set input/output projections
    inProj = Proj(init='epsg:4326')     #GCS/WGS84
    outProj = Proj(init='epsg:3410')    #NSIDC EASE-Grid
    
    # Convert ns_tolerance to projected meter coordinates
    lontmp,lattmp = transform(inProj,outProj,0,ns_tolerance)
    
    # Create y-coordinate array for vertex sampling
    y_den = np.arange(-lattmp,lattmp+proj_density,proj_density)
    
    pts = []
    
    for i in range(2):
        # Perform two passes, one south-to-north another north-to-south
        # If it is the second pass, flip array and walk down y-coordinate
        if i == 1:
            y_den = y_den[::-1]
            
        # Walk through y-coordinates and find the vertex per interval
        for j in range(y_den.size-1):
            #get min/max y-coordinate for the interval
            y1 = y_den[j]
            y2 = y_den[j+1]
            
            # First pass (south-to-north)
            if i == 0:
                idx = np.where((yy>=y1)&(yy<=y2))
                pts.append([xx[idx].max(),yy[idx].mean()])
            
            # Second pass (north-to-south)
            else:
                idx = np.where((yy>=y2)&(yy<=y1))
                pts.append([xx[idx].min(),yy[idx].mean()])
        
    return pts        
        

def geoverts_2_imgverts(xx,yy,poly):
    """
    FUNCTION:  verts = geoverts_2_imgverts(xx,yy,poly)
    ARGUMENTS: xx - x-coordinate grid for output raster
               yy - y-coordinate grid for output raster
               poly - swath polygon vertices
    KEYWORDS:  N/A
    RETURNS:   vertsOut - the polygon vertices calculated from the swath
    NOTES:     Takes polygon geographic coordinates and converts them into
               image pixel coordinates for creating a mask image.
    """ 
    
    vertsOut = []
    
    # Iterate through each polygon vertex
    for i in range(len(poly)):
        # Extract x-y coordinate from vertex
        xval = poly[i][0]
        yval = poly[i][1]
        
        # Find closest image index to the x-y coordinate
        xidx = (np.abs(xx-xval)).argmin()
        yidx = (np.abs(yy-yval)).argmin()
        
        # Convert the 1-d index to 2-d
        ridx = yidx / xx.shape[1]
        cidx = xidx % xx.shape[1]
        
        vertsOut.append((cidx,ridx))
        
    return vertsOut

# Define constants
res = 6000           # output map spatial resolution in meters
ns_tolerance = 90    # NS latitudes to clid data to
noData = -1          # Output NoData value

# Specify output file path
outfile = r'~/AMSR2_L2_RainOcean_gridded.tif'
    
# Path to input file
infile = r'~/AMSR_2_L2_RainOcean_R00_201702071903_D.he5'

# Open HDF5 file
h5file = tables.open_file(infile)

# Load data into memory
lats = h5file.get_node('/HDFEOS/SWATHS/GPROF2010V2/Geolocation Fields/Latitude').read()
lons = h5file.get_node('/HDFEOS/SWATHS/GPROF2010V2/Geolocation Fields/Longitude').read()
swath = h5file.get_node('/HDFEOS/SWATHS/GPROF2010V2/Data Fields/surfaceRain').read()

# Close HDF5 file
h5file.close()

# Clip the data to only within +-60 degrees latitude
idx = np.where((lats<=ns_tolerance)&(lats>=-ns_tolerance))
swath = swath[idx]
lats = lats[idx]
lons = lons[idx]

# Set input/output projections
inProj = Proj(init='epsg:4326')     #GCS/WGS84
outProj = Proj(init='epsg:3410')    #NSIDC EASE-Grid

# Reproject the geolocation coordinates
lonproj,latproj = transform(inProj,outProj,lons,lats)

# Create rater grid coordinates and mesh
latgrd = np.arange(latproj.min(),latproj.max(),res)
longrd = np.arange(lonproj.min(),lonproj.max(),res)
xx,yy = np.meshgrid(longrd,latgrd)

# Format geolocation coordinates for gridding
pts = np.zeros((lons.size,2))
pts[:,0] = lonproj.ravel()
pts[:,1] = latproj.ravel()

# Calculate the swath footprint polygon
poly = get_swath_poly(lonproj,latproj,ns_tolerance,density=1)

# Convert the convex hull vertices from geographic points to image pixel indices
verts = geoverts_2_imgverts(xx,yy,poly)

# Create a mask from the convex hull vertices
img = Image.new('L', (xx.shape[1], xx.shape[0]), 0)
ImageDraw.Draw(img).polygon(verts, fill=1)
mask = numpy.array(img)
mask = np.flipud(mask)

# Grid the swath data and orient it correctly
grdData = griddata(pts, swath.ravel(), (xx,yy), method='nearest')
grdData = np.flipud(grdData)

# Set bad data and data outside of the polygon to the NoData value
grdData[np.where(mask==0)] = noData
grdData[np.where(grdData<=0)] = noData

# Load GDAL GeoTIFF dataset driver
drv = gdal.GetDriverByName("GTiff")

# Create blank ouput dataset
dsOut = drv.Create(outfile, xx.shape[1],  xx.shape[0], 1, gdal.GDT_Float32)

# Specify and set output file's geotransform
gt = [xx.min(), res, 0, yy.max(), 0, -res]
dsOut.SetGeoTransform(gt)

# Specify the NSIDC's EASE-Grid (EPSG:3410) projection parameters
dest_wkt = '''PROJCS["NSIDC EASE-Grid Global",
                GEOGCS["Unspecified datum based upon the International 1924 Authalic Sphere",
                    DATUM["Not_specified_based_on_International_1924_Authalic_Sphere",
                        SPHEROID["International 1924 Authalic Sphere",6371228,0,
                            AUTHORITY["EPSG","7057"]],
                        AUTHORITY["EPSG","6053"]],
                    PRIMEM["Greenwich",0,
                        AUTHORITY["EPSG","8901"]],
                    UNIT["degree",0.01745329251994328,
                        AUTHORITY["EPSG","9122"]],
                    AUTHORITY["EPSG","4053"]],
                UNIT["metre",1,
                    AUTHORITY["EPSG","9001"]],
                PROJECTION["Cylindrical_Equal_Area"],
                PARAMETER["standard_parallel_1",30],
                PARAMETER["central_meridian",0],
                PARAMETER["false_easting",0],
                PARAMETER["false_northing",0],
                AUTHORITY["EPSG","3410"],
                AXIS["X",EAST],
                AXIS["Y",NORTH]]'''

# Set output dataset projection
dsOut.SetProjection(dest_wkt)

# Write gridded swath to the output file
bandOut=dsOut.GetRasterBand(1)
bandOut.SetNoDataValue(noData)
BandWriteArray(bandOut, grdData)

# Flush output dataset
dsOut = None
bandOut = None 

Results from test:

Output file needs to have spatial projection defined...issues with GDAL. Output projection is in NSIDC EASE-Grid Global (EPSG: 3410)