Saturday 8 July 2017

Patch processing of multiple HDF5 files using Python and Gdal (gdalwarp)


I am trying to extend this process


http://geoinformaticstutorial.blogspot.com/2014/06/reading-and-map-projecting-amsr-2-data.html



wherein a set of raster data with separate geolocation arrays are map re-projected using virtual rasters and gdal.


Goal: Create a global image mosaic from ~15 HDF5 "swath" grided data files. The final image must have rows = 4872 & columns = 11568. The latter is based on NSIDC EASE Grid Version 2 (http://nsidc.org/data/ease/ease_grid2.html)


The command line tags and options for gdalwarp are:


-geoloc -overwrite -co "COMPRESS=LZW" -tr 3000 3000 -srcnodata -9999 -dstnodata -9999 -r average-s_srs EPSG:4326 -t_srs '+proj=cea +x_0=0 +y_0=0 +lon_0=0 +lat_ts=30 +ellps=WGS84 +units=m +no_defs'


Each file has 6 data bands, thus at the end should get 6 different image mosaics. Files are ~2GB large.


Progress so far:


Able to open and read each file using h5py


indata = h5py.File(filename)
data_band_1 = indata["cell_sigma0_1"]
data_band_1 = np.where(data_band_1 <=0,-9999,data_band_1 )

data_band_1 [data_flg_band_1 == TRUE] = -9999

Save numpy arrays as Tiff using the following, but with out geotransform information


srs = osr.SpatialReference()
srs.SetWellKnownGeogCS("WGS84")
outFileName = mydates[iter_day]+'.tif'
print(outFileName)
gdalDriver = getGDALFormatFromExt(outFileName)
format="GTiff"
driver = gdal.GetDriverByName(format)

metadata = driver.GetMetadata()
newDataset = driver.Create(outFileName, inXSize, inYSize, numBands,gdal.GDT_Float32)
newDataset.SetProjection(srs.ExportToWkt())
newDataset.GetRasterBand(i+1).SetDescription(temp_name)
newDataset.GetRasterBand(i+1).WriteArray(temp_data)
newDataset.GetRasterBand(i+1).SetNoDataValue(nodat)
stat = newDataset.GetRasterBand(i+1).ComputeStatistics(0) newDataset.GetRasterBand(i+1).SetStatistics(stat[0], stat[1], stat[2], stat[3])
newDataset.FlushCache()
newDataset = None


Created virtual rasters (.vrt) similar to the method in the link above. But replaced the source data with the new Tiff.


Question: How can I merge all the .vrt files to a raster with a specific size (e.g., known rows and columns)


Note: the lat/lon steps are not necessary uniform (can't setup a world file)




No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...