Earthdata Forum

Hi. My code is
#!/usr/bin/env python
# coding: utf-8

import os
import re
import numpy as np
import pandas as pd
from pyhdf.SD import SD, SDC
import rasterio
from rasterio.transform import from_origin

# Define file path
file_path = r"C:\Users\Utel\Downloads\x\MOD09GA.A2024328.h22v06.061.2024329014416.NRT.hdf"

# Read HDF file
hdf = SD(file_path, SDC.READ)

# Extract datasets
datasets = hdf.datasets().keys()

# Get bounding box coordinates
bounding_objects = ["NORTHBOUNDINGCOORDINATE", "WESTBOUNDINGCOORDINATE",
"EASTBOUNDINGCOORDINATE", "SOUTHBOUNDINGCOORDINATE"]

# Initialize bounding coordinates
bounding_coordinates = {}
for bounding_object in bounding_objects:
match = re.search(rf"OBJECT={bounding_object}(.*?)END_OBJECT={bounding_object}",
getattr(hdf, "ArchiveMetadata.0").replace(" ", ""), re.DOTALL)

if match:
value = re.search(r'VALUE=([\d.]+)', match.group(1)).group(1)
bounding_coordinates[bounding_object] = float(value)

# Read the dataset
dataset = hdf.select("sur_refl_b05_1")

# Get dataset attributes
dataset_attributes = dataset.attributes()
scale_factor = dataset_attributes["scale_factor"]
add_offset = dataset_attributes["add_offset"]
fill_value = dataset_attributes["_FillValue"]

# Get dimensions
num_columns = len(dataset.get())
num_rows = len(dataset[0])

# Define bounding box
left_x = bounding_coordinates["WESTBOUNDINGCOORDINATE"]
right_x = bounding_coordinates["EASTBOUNDINGCOORDINATE"]
lower_y = bounding_coordinates["SOUTHBOUNDINGCOORDINATE"]
upper_y = bounding_coordinates["NORTHBOUNDINGCOORDINATE"]

# Define scale
scale_x = (right_x - left_x) / num_columns
scale_y = (upper_y - lower_y) / num_rows

# Generate coordinate arrays
lon = np.linspace(left_x, right_x, num_columns)
lat = np.linspace(lower_y, upper_y, num_rows)

# Create the coordinate grid and invert latitude
lon, lat = np.meshgrid(lon, lat[::-1]) # Invert lat during meshgrid creation

# Flatten the dataset and replace fill values with NaN
dataset = dataset[:, :].astype(np.double)
dataset[np.isin(dataset, [fill_value])] = np.nan

# Convert the dataset from digital to analytical
dataset = (dataset - add_offset) * scale_factor

# Export to GeoTIFF
output_tiff_path = r"C:\Users\Utel\Downloads\output_dataset.tif"
transform = from_origin(left_x, upper_y, scale_x, scale_y) # top-left corner and pixel sizes

with rasterio.open(
output_tiff_path,
'w',
driver='GTiff',
height=num_rows,
width=num_columns,
count=1,
dtype='float64', # or np.float32 based on your data type
crs='EPSG:4326', # or your specific CRS
transform=transform,
) as dst:
dst.write(dataset, 1) # write the dataset to the first band

print(f"GeoTIFF saved to {output_tiff_path}")
while the output should be oblique I wanna be like wgs84.

I believe your questions has been answered by the HDF Group here: viewtopic.php?t=6218 but let us know if you have any follow on questions. Thanks -- Danielle