MODISL2 hdf to tiff and reprojection

 Posts: 4
 Joined: Sun Sep 04, 2022 4:15 am America/New_York
MODISL2 hdf to tiff and reprojection
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
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 (103.42 KiB) Not viewed yet

 Failed Projection.JPG (58.13 KiB) Not viewed yet
Tags:

 Subject Matter Expert
 Posts: 65
 Joined: Fri Jun 03, 2022 10:54 am America/New_York
 Location: NASA GSFC
 Been thanked: 1 time
 Contact:
Re: MODISL2 hdf to tiff and reprojection
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
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

 Posts: 4
 Joined: Sun Sep 04, 2022 4:15 am America/New_York
Re: MODISL2 hdf to tiff and reprojection
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 earthdataNASA, then I used l2genSeadas (Aerosol Model9)for atmospheric correction, and output format was set as hdf file, the output file has some Rrs bands and a L2flag band, latitude and longitude, not projected.
> 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 earthdataNASA, then I used l2genSeadas (Aerosol Model9)for atmospheric correction, and output format was set as hdf file, the output file has some Rrs bands and a L2flag band, latitude and longitude, not projected.

 Subject Matter Expert
 Posts: 595
 Joined: Tue Feb 09, 2021 5:42 pm America/New_York
 Been thanked: 1 time
Re: MODISL2 hdf to tiff and reprojection
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/help8.2.0/desktop/Reprojection.html
Are you able to use the reproject tool in SeaDAS to create the GeoTiff file? https://seadas.gsfc.nasa.gov/help8.2.0/desktop/Reprojection.html

 Posts: 4
 Joined: Sun Sep 04, 2022 4:15 am America/New_York
Re: MODISL2 hdf to tiff and reprojection
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!
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!

 Subject Matter Expert
 Posts: 595
 Joined: Tue Feb 09, 2021 5:42 pm America/New_York
 Been thanked: 1 time
Re: MODISL2 hdf to tiff and reprojection
You can do reprojection on command line  https://seadas.gsfc.nasa.gov/help8.2.0/GptCookbook/gptCookbookReproject.html
Here is the GPT tutorial  https://seadas.gsfc.nasa.gov/help8.2.0/GptCookbook/gptTutorial.html
Here is the GPT tutorial  https://seadas.gsfc.nasa.gov/help8.2.0/GptCookbook/gptTutorial.html

 Posts: 4
 Joined: Sun Sep 04, 2022 4:15 am America/New_York
Re: MODISL2 hdf to tiff and reprojection
OB SeaDAS  xuanyang02 wrote:
That really helped me a lot, thank you very much!
That really helped me a lot, thank you very much!