0

To warp image to GeoTIFF using GCPs method, I used the pyresample package. As the result, I did not find any example to save this result as a GeoTIFF file.

from pyhdf.SD import SD, SDC
from pyresample import image, geometry,save_quicklook
file="C:/Users/TLQ/Desktop/11/MYD09.A2011093.0410.006.2015217030905.hdf"
f1 = SD(file, SDC.READ)
lat = f1.select('Latitude')
latitude = lat[:,:]
lon = f1.select('Longitude')
longitude = lon[:,:]
b1 = f1.select('BAND31')
AOP = b1[:,:]
swath_def = geometry.SwathDefinition(lons=longitude, lats=latitude)
swath_con = image.ImageContainerNearest(AOP, swath_def, radius_of_influence=10000)
2
  • What does swath_con.image_data return? A numpy array? Apr 2, 2018 at 13:31
  • No, it is <class 'pyresample.image.ImageContainerNearest'> Apr 2, 2018 at 23:39

1 Answer 1

0

We had a similar but not exactly the same problem, here is my solution.

import os
import h5py
import numpy
import rasterio
import pyresample
import subprocess
from IPython.display import Image
#%% read data
filename = "C:/Datasets/FY3D-MWRI/FY3D_MWRIA_GBAL_L1_20220301_0423_010KM_MS.HDF"
f_obj = h5py.File(filename)

dem_arr = f_obj["Calibration/DEM"][:]
lat_arr = f_obj["Geolocation/Latitude"][:]
lon_arr = f_obj["Geolocation/Longitude"][:]

#%% Resample to equal lat/lon grid
# resolution can be calculated by converting 1 arc second to approximately 30m
# area_extent: (lower_left_x, lower_left_y, upper_right_x, upper_right_y)
area_def = pyresample.create_area_def(
    units="degrees",
    resolution=0.09259,
    area_id="my_area",
    area_extent=[73, 10, 140, 55],
    projection={"proj": "longlat", "datum": "WGS84"},
    description="China 10kmx10km degree lat-lon grid",

)

swath_def = pyresample.geometry.SwathDefinition(lons=lon_arr, lats=lat_arr)

result = pyresample.kd_tree.resample_nearest(
    swath_def, dem_arr, area_def, fill_value=9999, radius_of_influence=10000)

#%% output to tiff
output_filename=  "./hello.tif"
outds = rasterio.open("hello.tif",
                      mode="w",
                      count=1,
                      nodata=9999,
                      driver="GTiff",
                      dtype=numpy.int16,
                      crs='+proj=latlong',
                      width=result.shape[1],
                      height=result.shape[0],
                      transform=[0.09259, 0, 73, 0, -0.09259, 55],
                      )
outds.write(result, 1)
outds.close()

The results can be plotted using GMT

plot the tif

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.