Reprojecting swath to lat lon grid?

Question

I have satellite swath data in the form of an HDF5 file containing 3 arrays: 'data_array', 'latitude' and 'longitude'. All three are 32-bit floats, are arrays of the same dimensions (X-by-Y) and the names are self-explanatory. All three datasets use a fill value of -32767

What I want to do it reproject data_array from native satellite swath onto a latitude/longitude grid with a pixel size of 0.009 degrees. I want to save the resulting data array as a geoTIFF containing project information (for use in QGis/ENVI/etc).

I have tried using gdal based upon this but gdalwarp segfaults (maybe because there are some fill value pixels in the lat/lon grid? Unsure how to deal with this).

I also tried using python and the pytroll libraries, but find the documentation somewhat impenetrable. All my attempts thus far have failed.

So, can anyone recommend a simple method to perform this reprojection?

I assumed it would be something quite common, but am having difficulty finding a solution.

I can plot project the data onto a grid using NASA Panoply, but this appears to be for visualisation only. I can't see how to save a GEOTiff output.

Answer 1

You might be interested by the Orfeo Toolbox which might have some useful feature for you (i only know about the tool, not the detailed procedure).

Orfeo ToolBox (OTB) is an open-source project for state-of-the-art remote sensing. Built on the shoulders of the open-source geospatial community, it can process high resolution optical, multispectral and radar images at the terabyte scale. A wide variety of applications are available: from ortho-rectification or pansharpening, all the way to classification, SAR processing, and much more!

All of OTB’s algorithms are accessible from Monteverdi, QGIS, Python, the command line or C++. Monteverdi is an easy to use visualization tool with an emphasis on hardware accelerated rendering for high resolution imagery (optical and SAR). With it, end-users can visualize huge raw imagery products and access all of the applications in the toolbox. From resource limited laptops to high performance MPI clusters, OTB is available on Windows, Linux and Mac. It is community driven, extensible and heavily documented.

https://www.orfeo-toolbox.org/features-2/

You say in one of your comment : "the data is in raw instrument swath projection". You might specifically be interested by the Orfeo Toolbox recipes structured in articles and with code examples :

From raw image to calibrated product

Pre-processing tasks

This section presents various pre-processing tasks that are presented in a classical order to obtain a calibrated, pan-sharpened image.

• Optical radiometric calibration

• Pan-sharpening

• Digital Elevation Model management

• Ortho-rectification and map projections

In the same area, one more page should help you out :

Residual registration

Image registration is a fundamental problem in image processing. The aim is to align two or more images of the same scene often taken at different times, from different viewpoints, or by different sensors. It is a basic step for orthorectification, image stitching, image fusion, change detection, and others. [...] Sensor model is generally not sufficient to provide image registrations. Indeed, several sources of geometric distortion can be contained in optical remote sensing images including earth rotation, platform movement, non linearity, etc.

They result in geometric errors on scene level, image level and pixel level. It is critical to rectify the errors before a thematic map is generated, especially when the remote sensing data need to be integrated together with other GIS data.

Answer 2

Try this with GDAL:

gdal_translate hdf5_file -sds bands.tif

Or this in R:

library(gdalUtils)
band=1
file_path<-"/myfolder/myfile.hdf"
system(paste0('gdal_translate ', get_subdatasets(file_path)[band],' myfile_', band','.tif'))

Answer 3

You could try georeferencing through thinplate spline (tps).

In Python should it be something like this code line:

import gdal
FOLD = 'C:/FOLD'
img1=gdal.Open(FOLD+"/Imagery.hdf5")
img2=gdal.Open(img1.GetSubDatasets()[3][0]) #where 3 is the band of interest (is an example)
gdal.Warp(FOLD+'/FILE.tif',img2,tps=True) #the -tps algorithm georreferences your data