MODISL2 hdf to tiff and reprojection

[zhanggaoming]

Hello, I am a newbie in using python to process MODIS data，now I am trying to covert MODISL2 hdf data(atmosphere corrected products using MOD02 and Seadas) to Geotiff files for futher band calculation.
I could not obtain geotransform parameters from MODISL2 hdf matadata, so I tried using geotransform parameters extract from a projected tif file(coverted by Seadas using same MODISL2 hdf), it was fine until I got a projection problem(Incorrect coordinate values), . My codes are as follows:

import gdal
import numpy as np
from numpy import shape
from pyhdf.SD import *
from osgeo import osr


in_tiff_path = r'D:\PythonTest\TSS calculation'
hdf_obj = SD('MODIS_2008097_SWIR.hdf')

# get attributes, layers, and single layer of Rrs
attr = hdf_obj.attributes()
dateset = hdf_obj.datasets()

columns = int(attr['Number of Pixel Control Points'])
rows = int(attr['Number of Scan Control Points'])

Rrs = hdf_obj.select('Rrs_1240').get()
d_Rrs = Rrs.astype(np.float32)

# convert negative values of hdf to real pixel value array
real_Rrs = (d_Rrs + 25000) / 500000
real_Rrs[real_Rrs < 0] = 0

# L2_flag
l2 = hdf_obj.select('l2_flags').get()
d_l2 = l2.astype(np.int32)

# get geotransform from tiff file converted by seadas
inpath = r'D:\PythonTest\TSS calculation\MODIS_2008097_SWIR.tif'
dataset = gdal.Open(inpath)

# creat Geotif file, set geotransform and projection
driver = gdal.GetDriverByName("GTiff")
out_tiff = driver.Create("MODIS_2008097.tif", columns, rows, 2, gdal.GDT_Float32)

out_tiff.SetGeoTransform(dataset.GetGeoTransform())
out_tiff.SetProjection('GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG",''"7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,''AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG",''"4326"]]')

# get bands
out_tiff.GetRasterBand(1).WriteArray(real_Rrs)
out_tiff.GetRasterBand(2).WriteArray(d_l2)

out_tiff.FlushCache()

Unfortunately, I got a result with failed projection.（Failed Projection.jpg）
Then I tried to calculate geotransform parameters manually using the image four corner coordinates from the HDF metadata and least squares method, then used this parameter to projection. Codes for calculation of geotransform:

# Calculation of affine parameter(Geo transformation)
class Point(object):
def __init__(self, x, y):
self.x = x
self.y = y

def affine_abc(pixel_points, geo_points):
n = len(pixel_points)
pixelX_square = 0.0
pixelX_pixelY = 0.0
pixelY_square = 0.0
pixelX = 0.0
pixelY = 0.0
pixelX_geoX = 0.0
pixelY_geoX = 0.0
geoX = 0.0
for i in range(0, n):
pixelX_square += math.pow(pixel_points[i].x, 2)
pixelX_pixelY += pixel_points[i].x * pixel_points[i].y
pixelY_square += math.pow(pixel_points[i].y, 2)
pixelX += pixel_points[i].x
pixelY += pixel_points[i].y
pixelX_geoX += pixel_points[i].x * geo_points[i].x
pixelY_geoX += pixel_points[i].y * geo_points[i].x
geoX += geo_points[i].x
a = np.array([[pixelX_square, pixelX_pixelY, pixelX], [pixelX_pixelY, pixelY_square, pixelY], [pixelX, pixelY, n]])
b = np.array([[pixelX_geoX], [pixelY_geoX], [geoX]])
at = np.linalg.inv(a)
result = at.dot(b)
return result[0, 0], result[1, 0], result[2, 0]


def affine_def(pixel_points, geo_points):
n = len(pixel_points)
pixelX_square = 0.0
pixelX_pixelY = 0.0
pixelY_square = 0.0
pixelX = 0.0
pixelY = 0.0
pixelX_geoY = 0.0
pixelY_geoY = 0.0
geoY = 0.0
for i in range(0, n):
pixelX_square += math.pow(pixel_points[i].x, 2)
pixelX_pixelY += pixel_points[i].x * pixel_points[i].y
pixelY_square += math.pow(pixel_points[i].y, 2)
pixelX += pixel_points[i].x
pixelY += pixel_points[i].y
pixelX_geoY += pixel_points[i].x * geo_points[i].y
pixelY_geoY += pixel_points[i].y * geo_points[i].y
geoY += geo_points[i].y
a = np.array([[pixelX_square, pixelX_pixelY, pixelX], [pixelX_pixelY, pixelY_square, pixelY], [pixelX, pixelY, n]])
b = np.array([[pixelX_geoY], [pixelY_geoY], [geoY]])
at = np.linalg.inv(a)
result = at.dot(b)
return result[0, 0], result[1, 0], result[2, 0]


if __name__ == '__main__':
p0 = Point(0.0, 0.0)
p1 = Point(1353, 0.0)
p2 = Point(0.0, 2029)
p3 = Point(1353.0, 2029.0)
g0 = Point(attr['Upper Left Longitude'], attr['Upper Left Latitude'])
g1 = Point(attr['Upper Right Longitude'], attr['Upper Right Latitude'])
g2 = Point(attr['Lower Left Longitude'], attr['Lower Left Latitude'])
g3 = Point(attr['Lower Right Longitude'], attr['Lower Right Latitude'])
ps = [p0, p1, p2, p3]
gs = [g0, g1, g2, g3]
abc = affine_abc(ps, gs)
def0 = affine_def(ps, gs)
# Geo transformation
affine_MODIS = (abc[2], abc[0], abc[1], def0[2], def0[0], def0[1])


The resulting image looks better this time, but upon closer inspection I fuond that the coordinates are still incorrect.(Failed Projection2.jpg)
I would be grateful if someone could provide some advice or hints(for example on how to calculate the MODIS geotransform parameters correctly, etc.)
thanks

[OB General Science - guoqingw]

Hi,

Do you mind providing more information about the MODIS L2 data that you are using? Where did you download the file? What is the original MODIS L2 filename?

Thank you,
Guoqing

[zhanggaoming]

OB General Science - guoqingw wrote:
> Hi,
>
> Do you mind providing more information about the MODIS L2 data that you are
> using? Where did you download the file? What is the original MODIS L2
> filename?
>
> Thank you,
> Guoqing

Hi, thank you for your reply, firstly I download MOD021KM V6 from earthdata-NASA, then I used l2gen-Seadas (Aerosol Model-9)for atmospheric correction, and output format was set as hdf file, the output file has some Rrs bands and a L2-flag band, latitude and longitude, not projected.

[OB SeaDAS - xuanyang02]

So MODIS_2008097_SWIR.hdf is the L2 file created by l2gen? What SeaDAS tool did you use to create the tiff file MODIS_2008097.tif?

Are you able to use the reproject tool in SeaDAS to create the GeoTiff file? https://seadas.gsfc.nasa.gov/help-8.2.0/desktop/Reprojection.html

[zhanggaoming]

OB SeaDAS - xuanyang02 wrote:

Yes, exactly I used reprojection tool to create that tif file to retrive geotransform parameters, since I got thousands of image to process, I can not finish them one by one using Seadas GUI. That's why I try to find
geotransform parameters, I will be very grateful if you know some better way for batch reprojection, thank you!

[OB SeaDAS - xuanyang02]

You can do reprojection on command line -- https://seadas.gsfc.nasa.gov/help-8.2.0/GptCookbook/gptCookbookReproject.html

Here is the GPT tutorial -- https://seadas.gsfc.nasa.gov/help-8.2.0/GptCookbook/gptTutorial.html

[zhanggaoming]

OB SeaDAS - xuanyang02 wrote:

That really helped me a lot, thank you very much!