MODISL2 hdf to tiff and reprojection

Use this Forum to find information on, or ask a question about, NASA Earth Science data.
Post Reply
zhanggaoming
Posts: 4
Joined: Sun Sep 04, 2022 4:15 am America/New_York
Answers: 0

MODISL2 hdf to tiff and reprojection

by zhanggaoming » Sun Sep 04, 2022 5:08 am America/New_York

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
Attachments
Failed Projection2.JPG
Failed Projection2.JPG (103.42 KiB) Not viewed yet
Failed Projection.JPG
Failed Projection.JPG (58.13 KiB) Not viewed yet

Filters:

OB General Science - guoqingw
Subject Matter Expert
Subject Matter Expert
Posts: 79
Joined: Fri Jun 03, 2022 10:54 am America/New_York
Answers: 0
Location: NASA GSFC
Been thanked: 1 time
Contact:

Re: MODISL2 hdf to tiff and reprojection

by OB General Science - guoqingw » Thu Sep 08, 2022 11:19 am America/New_York

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
Posts: 4
Joined: Sun Sep 04, 2022 4:15 am America/New_York
Answers: 0

Re: MODISL2 hdf to tiff and reprojection

by zhanggaoming » Thu Sep 08, 2022 10:02 pm America/New_York

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
Subject Matter Expert
Subject Matter Expert
Posts: 657
Joined: Tue Feb 09, 2021 5:42 pm America/New_York
Answers: 1
Been thanked: 2 times

Re: MODISL2 hdf to tiff and reprojection

by OB SeaDAS - xuanyang02 » Fri Sep 09, 2022 10:11 am America/New_York

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
Posts: 4
Joined: Sun Sep 04, 2022 4:15 am America/New_York
Answers: 0

Re: MODISL2 hdf to tiff and reprojection

by zhanggaoming » Fri Sep 09, 2022 9:17 pm America/New_York

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
Subject Matter Expert
Subject Matter Expert
Posts: 657
Joined: Tue Feb 09, 2021 5:42 pm America/New_York
Answers: 1
Been thanked: 2 times

Re: MODISL2 hdf to tiff and reprojection

by OB SeaDAS - xuanyang02 » Mon Sep 12, 2022 10:59 am America/New_York


zhanggaoming
Posts: 4
Joined: Sun Sep 04, 2022 4:15 am America/New_York
Answers: 0

Re: MODISL2 hdf to tiff and reprojection

by zhanggaoming » Thu Sep 15, 2022 10:08 pm America/New_York

OB SeaDAS - xuanyang02 wrote:

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

Post Reply