01 Description
It contains three pro processes, of which modis_grid_warp is the main program, modis_grid_geo_transform and img_warp are regular pro processes, and read_h4 is a custom function.
modis_grid_geo_transform generates a projection coordinate data set based on corner information and converts it into a longitude and latitude data set, and img_warp performs reprojection (GLT correction) based on the longitude and latitude data set
02 complete code
2.1 modis_grid_geo_transform process
; +
; Function usage:
; This program is used to read the corner information of MODIS GRID products and return the latitude and longitude data set
; Function parameters:
; h4_path: path of MODIS GRID file
; extract_lon: the returned longitude data set
; extract_lat: returned latitude data set
; ds_size: (keyword parameter) binary array containing column number and row number, if not passed in, read from file
; helpme: View function usage
;-
pro modis_grid_geo_transform, h4_path, extract_lon, extract_lat, ds_size=ds_size
; Get metadata
h4_id = hdf_sd_start(h4_path, /read)
metadata_ix = hdf_sd_attrfind(h4_id, 'StructMetadata.0')
hdf_sd_attrinfo, h4_id, metadata_ix, data=metadata
; Get the row and column number
if ~keyword_set(ds_size) then begin
ds_size = make_array(2, /double)
ix_one = strpos(metadata, 'XDim')
ix_two = strpos(metadata, 'YDim')
ix_three = strpos(metadata, 'UpperLeftPointMtrs')
ds_size[0] = (strsplit(strmid(metadata, ix_one, ix_two - ix_one), '=', /extract))[1]
ds_size[1] = (strsplit(strmid(metadata, ix_two, ix_three - ix_two), '=', /extract))[1]
endif else begin
ds_size = double(ds_size); Convert to double precision to avoid subsequent out-of-range
endelse
; Get corner point information
ix_one = strpos(metadata, 'UpperLeftPointMtrs')
ix_two = strpos(metadata, 'LowerRightMtrs')
ix_three = strpos(metadata, 'Projection')
ul_coor = strsplit(strmid(metadata, ix_one, ix_two - ix_one), '=(,)', /extract)
dr_coor = strsplit(strmid(metadata, ix_two, ix_three - ix_two), '=(,)', /extract)
ul_coor = double(ul_coor[1:2])
dr_coor = double(dr_coor[1:2])
x_res = (dr_coor[0] - ul_coor[0]) / ds_size[0]
y_res = (ul_coor[1] - dr_coor[1]) / ds_size[1]
; Calculation-projected coordinate data set
x_prj_coor = make_array(ds_size, /double)
y_prj_coor = make_array(ds_size, /double)
for ix=0, ds_size[0] - 1 do x_prj_coor[ix, *] = ul_coor[0] + x_res * ix + x_res / 2.0; represents the center position of the pixel
for ix=0, ds_size[1] - 1 do y_prj_coor[*, ix] = ul_coor[1] - y_res * ix - y_res / 2.0
x_prj_coor = reform(x_prj_coor, ds_size[0] * ds_size[1], /overwrite); /overwrite avoids copying and speeds up copying
y_prj_coor = reform(y_prj_coor, ds_size[0] * ds_size[1], /overwrite)
; The reform is performed to adapt to the input of the map_proj_inverse function: An n-element vector containing the x values.
; Projected coordinate data set ==> Latitude and longitude data set
prj_info = map_proj_init('Sinusoidal', /gctp, sphere_radius=6371007.181, $
center_longitude=0.0, false_easting=0.0, false_northing=0.0)
; It should be noted that gctp represents U.S. Geological Survey's General Cartographic Transformation Package (GCTP) for projection rather than IDL's own projection library
geo_coor = map_proj_inverse(x_prj_coor, y_prj_coor, map_structure=prj_info)
extract_lon = geo_coor[0, *].reform(ds_size)
extract_lat = geo_coor[1, *].reform(ds_size)
end
2.2 img_warp process
; +
; Function usage:
; Geometrically correct the target dataset based on the latitude and longitude dataset (reprojection-WGS84)
; Function parameters:
; target: target data set to be corrected
; lon: Longitude data set corresponding to the target data set
; lat: latitude data set corresponding to the target data set
; out_res: output resolution (°)
; target_warped: corrected target data set
; geo_info: the geographical structure corresponding to the corrected target data set
; degree<keyword argument: 5>: degree of polynomial
; interp<keyword parameter: nearest>: interpolation algorithm (including: nearest, linear, cublic)
; sub_percent<keyword parameter: 0.1>: By default, 10% of the points are used for geometric correction
;-
pro img_warp, target, lon, lat, out_res, target_warped, geo_info, degree=degree, interp=interp, $
sub_percent=sub_percent
; Get basic information
ds_size = size(target)
lon_lat_corner = hash($; considering that min()max() is the center of the corner pixel rather than the boundaries of the four corner points
'min_lon', min(lon) - out_res / 2.0d, $
'max_lon', max(lon) + out_res / 2.0d, $
'min_lat', min(lat) - out_res / 2.0d, $
'max_lat', max(lat) + out_res / 2.0d)
col_count_out = ceil((lon_lat_corner['max_lon'] - lon_lat_corner['min_lon']) / out_res)
row_count_out = ceil((lon_lat_corner['max_lat'] - lon_lat_corner['min_lat']) / out_res)
interp_code = hash($
'nearest', 0, $
'linear', 1, $
'cublic', 2)
if ~keyword_set(interp) then interp='nearest'
if ~keyword_set(degree) then degree=5
if ~keyword_set(sub_percent) then sub_percent = 0.1
; Original row and column number matrix
row_ori = make_array(ds_size[1:2], /integer)
col_ori = make_array(ds_size[1:2], /integer)
for ix=0, ds_size[1]-1 do col_ori[ix, *] = ix
for ix=0, ds_size[2]-1 do row_ori[*, ix] = ix
; Corrected row and column number matrix
col_warp = floor((lon - lon_lat_corner['min_lon']) / out_res)
row_warp = floor((lon_lat_corner['max_lat'] - lat) / out_res)
; Get the uniform sampling point index of sub_percent
all_count = ds_size[1] * ds_size[2]
sample_count = floor(all_count * sub_percent)
sample_ix = floor((findgen(sample_count) / double(sample_count)) * all_count)
polywarp, col_ori[sample_ix], row_ori[sample_ix], col_warp[sample_ix], row_warp[sample_ix], $
degree, k_col, k_row
target_warped = poly_2d(target, k_col, k_row, interp_code[interp], $
col_count_out, row_count_out, missing=!values.F_NAN)
geo_info={<!-- -->$
MODELPIXELSCALETAG: [out_res, out_res, 0.0], $ ; resolution
MODELTIEPOINTTAG: [0.0, 0.0, 0.0, lon_lat_corner['min_lon'], lon_lat_corner['max_lat'], 0.0], $ ; Corner information
GTMODELTYPEGEOKEY: 2, $ ; Set to geographic coordinate system
GTRASTERTYPEGEOKEY: 1, $ ; pixel representation type, North-Up image (North-Up)
GEOGRAPHICTYPEGEOKEY: 4326, $ ; The geographical coordinate system is WGS84
GEOGCITATIONGEOKEY: 'GCS_WGS_1984', $
GEOGANGULARUNITSGEOKEY: 9102} ; Units are degrees
end
2.3 read_h4 function
; +
; Function usage:
; Dataset and attributes used to read HDF4 files
; Function parameters:
; h4_path: The path of the HDF4 file to be read
; ds_name: name of the data set to be read
; attr_name<keyword parameter>: If the attribute name is passed in, the attribute information under the data set will be returned.
; double<keyword argument, /double>: if specified, converts the returned data set to double precision floating point type
;-
function read_h4, h4_path, ds_name, attr_name=attr_name, double=double
; This function is used to read HDF4 file data sets or attributes
h4_id = hdf_sd_start(h4_path, /read)
ds_index = hdf_sd_nametoindex(h4_id, ds_name)
ds_id = hdf_sd_select(h4_id, ds_index)
if keyword_set(attr_name) then begin
attr_index = hdf_sd_attrfind(ds_id, attr_name)
hdf_sd_attrinfo, ds_id, attr_index, data=data
endif else begin
hdf_sd_getdata, ds_id, data
if keyword_set(double) then data = double(data)
endelse
hdf_sd_endaccess, ds_id
hdf_sd_end, h4_id
return, data
end
2.4 modis_grid_warp process (main program)
; @Author : ChaoQiezi ; @Time : October 19, 2023 - 4:26:55 pm ; @Email: [email protected] ; This program is used to reproject the MODIS GRID data set pro modis_grid_warp ; Prepare in_dir = 'D:\Objects\JuniorFallTerm\IDLProgram\Experiments\ExperimentalData\chapter_3\MYD11A1\' out_dir = 'D:\Objects\JuniorFallTerm\IDLProgram\Experiments\ExperimentalData\chapter_3\IDL_out_grid_warp_me\' if ~file_test(out_dir, /directory) then file_mkdir, out_dir lst_name = 'LST_Day_1km' range_name = 'valid_range' sf_name = 'scale_factor' out_res = 0.01 ; degree h4_paths = file_search(in_dir, '*.hdf', count=h4_paths_count) foreach h4_path, h4_paths do begin start_time = systime(1) ; Get data sets and attributes, etc. lst = read_h4(h4_path, lst_name, /double) valid_range = read_h4(h4_path, lst_name, attr_name=range_name) scale_factor = read_h4(h4_path, lst_name, attr_name=sf_name) lst_size = size(lst) ; LST data preprocessing invalid_ix = where((lst lt valid_range[0]) or (lst gt valid_range[1])) lst[invalid_ix] = !values.F_NAN lst = lst * scale_factor[0] ; Get the latitude and longitude data set modis_grid_geo_transform, h4_path, lon, lat, ds_size=lst_size[1:2] ; Perform geometric correction (reprojection) img_warp, lst, lon, lat, out_res, lst_warped, geo_info ; Output out_path = out_dir + file_basename(h4_path, '.hdf') + '.tiff' write_tiff, out_path, lst_warped, geotiff=geo_info, /float print, file_basename(h4_path, '.hdf'), systime(1) - start_time, format='%s: %0.2f s' endforeach end
Bye~