I had created a code in python that converts the hdf file into nc file and clips it for a defined region. When I plot in the python everything seems alright. But if I import the .nc file into QGIS, the extent of the image is very large.
Plot in QGIS:
Plot in Python:
The code is as follows:
import numpy as np
import xarray as xr
from pyhdf.SD import SD, SDC
import rioxarray
from netCDF4 import Dataset as NetCDFFile
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import datetime
from shapely.geometry import mapping
import geopandas
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
#shp = geopandas.read_file('E:/Masters/Thesis/Data/MODIS/shapefiles/Ganges.shp')
# Open the HDF4 file
file_path = 'MOD11C3.A2021001.061.2021043155222.hdf' # replace with your file path
hdf_file = SD(file_path, SDC.READ)
# Read a specific dataset
dataset_name = 'LST_Day_CMG' # example dataset name, replace as needed
sds_obj = hdf_file.select(dataset_name)
data = sds_obj[:]
data = data * 0.02
data = np.ma.masked_where(data == 0, data) # Mask invalid data (e.g., 0 or any specific fill value)
data = np.expand_dims(data,axis=0)
#Extract time
# Extract the date part
date_str = file_path.split('.')[1][1:] # '2021001'
# Extract year and Julian day
year = int(date_str[:4])
julian_day = int(date_str[4:])
# Convert Julian day to a regular date
date = datetime.datetime(year, 1, 1) + datetime.timedelta(days=julian_day - 1)
# Print the result
print("Extracted date:", date)
ref_date = datetime.datetime(1980,1,1,0,0,0)
print(ref_date)
base_time = date - ref_date
sec_base_time = base_time.total_seconds()
# Generate the latitude and longitude arrays
lats = np.linspace(90, -89.95, data.shape[1])
lons = np.linspace(-179.95, 180, data.shape[2])
# Explore the file structure
datasets_dict = hdf_file.datasets()
print("Datasets in the HDF4 file:")
for idx, sds in enumerate(datasets_dict.keys()):
print(idx, sds)
# Create an xarray DataArray
time = np.array([sec_base_time])
data_array = xr.DataArray(
data,
name = "LST",
dims=["time","y", "x"],
coords={
"time": time,
"y": lats,
"x": lons
},
attrs={
"units": "K", # replace with appropriate units if different
"long_name": "Land Surface Temperature"
}
)
# Add attributes to the time coordinate
data_array.coords['time'].attrs = {
"units": "seconds since 1980-01-01 00:00:00",
"calendar": "gregorian",
"long_name": "time"
}
# #export to nc the whole thing
# # Set the _FillValue attribute
# data_array.attrs["_FillValue"] = np.nan
# # Write the DataArray to a NetCDF file with compression
# file_path_nc = 'check.nc'
# encoding = {data_array.name: {'zlib': True}} # Specify compression for the DataArray
# data_array.to_netcdf(file_path_nc, encoding=encoding)
# Define the range of latitudes and longitudes you want to clip
lat_min = 32 # replace with your minimum latitude
lat_max = 21 # replace with your maximum latitude
lon_min = 76 # replace with your minimum longitude
lon_max = 92 # replace with your maximum longitude
# Clip the DataArray
clipped = data_array.sel(y=slice(lat_min, lat_max), x=slice(lon_min, lon_max))
#print(clipped)
# Set the _FillValue attribute
clipped.attrs["_FillValue"] = np.nan
# Write the DataArray to a NetCDF file with compression
file_path_nc = f'LST_day_{sec_base_time}_check.nc'
#file_path_nc = 'check1.nc'
encoding = {clipped.name: {'zlib': True}} # Specify compression for the DataArray
clipped.to_netcdf(file_path_nc, encoding=encoding)
#ploting nc file
file_path_nc = f'LST_day_{sec_base_time}_check.nc'
#file_path_nc = 'check1.nc' # replace with your desired output file path
#file_path_nc = f'LST_day_{sec_base_time}_check.nc'
data_array = xr.open_dataarray(file_path_nc)
#if data_array.lat[0] < data_array.lat[-1]:
# data_array = data_array.reindex(lat=list(reversed(data_array.lat)))
lons = data_array.x
lats = data_array.y
data = data_array.values.squeeze()
time = data_array.time
# Plot the data
fig, ax = plt.subplots(figsize=(10, 5), subplot_kw={'projection': ccrs.PlateCarree()})
# Add coastlines and borders for context
ax.coastlines()
ax.add_feature(cfeature.BORDERS, linestyle=':')
# Plot the data
c = ax.pcolormesh(lons, lats, data, cmap='viridis', transform=ccrs.PlateCarree())
# Add a colorbar
fig.colorbar(c, ax=ax, orientation='vertical', label=data_array.units)
# Add title and labels
ax.set_title(data_array.long_name)
ax.set_xlabel('Longitude')
ax.set_ylabel('Latitude')
plt.show()