读取和绘制 HDF5(.h5) 文件并显示特定纬度和经度的地图
Reading and plotting HDF5(.h5) file and showing map for a specific latitude and longitude
我从 MOSDAC 网站获得了一个名为 3DIMG_30MAR2018_0000_L1B_STD.h5 的 HDF5 文件，我想读取这个文件并在地图中可视化文件中的 TIR 和 VIS 数据集。由于我是编程新手，所以我不太清楚如何使用纬度和经度值对全球数据的特定位置执行此操作。我的问题陈述是显示印度地区的数据。
到目前为止我完成的代码：
#Importing libraries`
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import h5py
import os
import netCDF4 as nc
import seaborn as sns
import geopandas as gpd
#Reading and listing keys present in HDF
fn = '3DIMG_30MAR2018_0000_L1B_STD.h5' #filename (the ".h5" file)
f = h5py.File(fn)
list(f.keys())
#extract data from the HDF5 file
mir = f['IMG_MIR'][:] #middle infra-red
swir = f['IMG_SWIR'][:] #shortwave IR
tir1 = f['IMG_TIR1'][:] #Thermal IR1
tir2 = f['IMG_TIR2'][:] #Thermal IR2
vis = f['IMG_VIS'][:] #Visible count
wv = f['IMG_WV'][:] #Water vapor count
lat = f['Latitude'][:]
lon = f['Longitude'][:]
lat_vis = f['Latitude_VIS'][:]
lon_vis = f['Longitude_VIS'][:]
lat_wv = f['Latitude_WV'][:]
lon_wv = f['Longitude_WV'][:]
vis = f['IMG_MIR'][0,:,:]
fig,ax = plt.subplots(figsize=(10,10))
im1 = plt.imshow(vis)
plt.colorbar(im1)
tir1d = f['/IMG_TIR1'];
lat = f['/Latitude'];
lon = f['/Longitude'];
Lat = np.linspace(81.041527, -81.041527, 2805)
Lon = np.linspace(0.84329641, 163.15671, 2816)
shp = gpd.read_file(r'C:\Users\offic\shape\IND_adm0.shp')
%lat=0:gs:60;%INDIA
%lon=60:gs:100;
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
data = Dataset(r'C:\Users\officDIMG_30MAR2018_0000_L1B_STD.h5')
data
Out: <class 'netCDF4._netCDF4.Dataset'>
root group (NETCDF4 data model, file format HDF5):
conventions: CF-1.6
title: 3DIMG_30MAR2018_0000_L1B
institute: BES,SAC/ISRO,Ahmedabad,INDIA.
source: BES,SAC/ISRO,Ahmedabad,INDIA.
Unique_Id: 3DIMG_30MAR2018_0000
Satellite_Name: INSAT-3D
Sensor_Id: IMG
Sensor_Name: IMAGER
HDF_Product_File_Name: 3DIMG_30MAR2018_0000_L1B_STD.h5
Output_Format: hdf5-1.8.8
Station_Id: BES
Ground_Station: BES,SAC/ISRO,Ahmedabad,INDIA.
Product_Type: STANDARD(FULL DISK)
Processing_Level: L1B
Imaging_Mode: FULL FRAME
Acquisition_Date: 30MAR2018
Acquisition_Time_in_GMT: 0000
Acquisition_Start_Time: 30-MAR-2018T00:00:08
Acquisition_End_Time: 30-MAR-2018T00:26:58
Product_Creation_Time: 2018-03-30T06:03:39
Radiometric_Calibration_Type: LAB CALIBRATED
Nominal_Altitude(km): 36000.0
Observed_Altitude(km): 35787.8215
Field_of_View(degrees): 17.973925
Nominal_Central_Point_Coordinates(degrees)_Latitude_Longitude: [ 0. 82.]
Attitude_Source: STAR
Sun_Azimuth(Degrees): 86.314133
Sat_Azimuth(Degrees): 286.793152
Sat_Elevation(Degrees): 89.813721
Sun_Elevation(Degrees): -6.010917
FastScan_Linearity_Enabled: no
SlowScan_Linearity_Enabled: no
MCD_FS_Enabled: no
MCD_SS_Enabled: no
Yaw_Flip_Flag: Y
Software_Version: 1.0
TIR1_Gain_Mode: 2
TIR2_Gain_Mode: 2
MIR_Gain_Mode: 2
WV_Gain_Mode: 3
VIS_Gain_Mode: 2
SWIR_Gain_Mode: 3
TIR1_Acquisition_Mode: MAIN
TIR2_Acquisition_Mode: MAIN
MIR_Acquisition_Mode: MAIN
WV_Acquisition_Mode: MAIN
VIS_Acquisition_Mode: MAIN
SWIR_Acquisition_Mode: MAIN
left_longitude: 0.8432964
right_longitude: 163.15671
upper_latitude: 81.04153
lower_latitude: -81.04153
Datum: WGS84
Ellipsoid: WGS84
dimensions(sizes): GeoX(2805), GeoX1(1402), GeoX2(11220), GeoY(2816), GeoY1(1408), GeoY2(11264), GreyCount(1024), time(1)
variables(dimensions): int32 GeoX(GeoX), int32 GeoX1(GeoX1), int32 GeoX2(GeoX2), int32 GeoY(GeoY), int32 GeoY1(GeoY1), int32 GeoY2(GeoY2), int32 GreyCount(GreyCount), uint16 IMG_MIR(time, GeoY, GeoX), float32 IMG_MIR_RADIANCE(GreyCount), float32 IMG_MIR_TEMP(GreyCount), uint16 IMG_SWIR(time, GeoY2, GeoX2), float32 IMG_SWIR_RADIANCE(GreyCount), uint16 IMG_TIR1(time, GeoY, GeoX), float32 IMG_TIR1_RADIANCE(GreyCount), float32 IMG_TIR1_TEMP(GreyCount), uint16 IMG_TIR2(time, GeoY, GeoX), float32 IMG_TIR2_RADIANCE(GreyCount), float32 IMG_TIR2_TEMP(GreyCount), uint16 IMG_VIS(time, GeoY2, GeoX2), float32 IMG_VIS_ALBEDO(GreyCount), float32 IMG_VIS_RADIANCE(GreyCount), uint16 IMG_WV(time, GeoY1, GeoX1), float32 IMG_WV_RADIANCE(GreyCount), float32 IMG_WV_TEMP(GreyCount), int16 Latitude(GeoY, GeoX), int32 Latitude_VIS(GeoY2, GeoX2), int16 Latitude_WV(GeoY1, GeoX1), int16 Longitude(GeoY, GeoX), int32 Longitude_VIS(GeoY2, GeoX2), int16 Longitude_WV(GeoY1, GeoX1), <class 'str'> SCAN_LINE_TIME(GeoY1), uint16 Sat_Azimuth(time, GeoY, GeoX), int16 Sat_Elevation(time, GeoY, GeoX), uint16 Sun_Azimuth(time, GeoY, GeoX), int16 Sun_Elevation(time, GeoY, GeoX), float64 time(time)
groups:
data.variables
Out: {'GeoX': <class 'netCDF4._netCDF4.Variable'>
int32 GeoX(GeoX)
unlimited dimensions:
current shape = (2805,)
filling off,
'GeoX1': <class 'netCDF4._netCDF4.Variable'>
int32 GeoX1(GeoX1)
unlimited dimensions:
current shape = (1402,)
filling off,
'GeoX2': <class 'netCDF4._netCDF4.Variable'>
int32 GeoX2(GeoX2)
unlimited dimensions:
current shape = (11220,)
filling off,
'GeoY': <class 'netCDF4._netCDF4.Variable'>
int32 GeoY(GeoY)
unlimited dimensions:
current shape = (2816,)
filling off,
'GeoY1': <class 'netCDF4._netCDF4.Variable'>
int32 GeoY1(GeoY1)
unlimited dimensions:
current shape = (1408,)
filling off,
'GeoY2': <class 'netCDF4._netCDF4.Variable'>
int32 GeoY2(GeoY2)
unlimited dimensions:
current shape = (11264,)
filling off,
'GreyCount': <class 'netCDF4._netCDF4.Variable'>
int32 GreyCount(GreyCount)
unlimited dimensions:
current shape = (1024,)
filling off,
'IMG_MIR': <class 'netCDF4._netCDF4.Variable'>
uint16 IMG_MIR(time, GeoY, GeoX)
long_name: Middle Infrared Count
invert: true
central_wavelength: 3.9313
bandwidth: 0.2
wavelength_unit: micron
bits_per_pixel: 10
resolution: 4.0
resolution_unit: km
_FillValue: 1023
lab_radiance_scale_factor: 0.000283937
lab_radiance_add_offset: -0.00401529
lab_radiance_quad: 0.0
lab_radiance_scale_factor_gsics: 0.000366208
lab_radiance_add_offset_gsics: -0.00833466
lab_radiance_quad_gsics: 0.0
radiance_units: mW.cm-2.sr-1.micron-1
coordinates: time Latitude Longitude
unlimited dimensions:
current shape = (1, 2816, 2805)
filling on,
'IMG_MIR_RADIANCE': <class 'netCDF4._netCDF4.Variable'>
float32 IMG_MIR_RADIANCE(GreyCount)
long_name: Middle Infrared Radiance
invert: true
units: mW.cm-2.sr-1.micron-1
_FillValue: 999.0
unlimited dimensions:
current shape = (1024,)
filling on,
'IMG_MIR_TEMP': <class 'netCDF4._netCDF4.Variable'>
float32 IMG_MIR_TEMP(GreyCount)
long_name: Middle Infrared Brightness Temperature
invert: true
units: K
_FillValue: 999.0
unlimited dimensions:
current shape = (1024,)
filling on,
'IMG_SWIR': <class 'netCDF4._netCDF4.Variable'>
uint16 IMG_SWIR(time, GeoY2, GeoX2)
long_name: Shortwave Infrared Count
invert: false
central_wavelength: 1.625
bandwidth: 0.15
wavelength_unit: micron
bits_per_pixel: 10
resolution: 1.0
resolution_unit: km
_FillValue: 0
lab_radiance_scale_factor: 0.007835
lab_radiance_add_offset: -0.172166
lab_radiance_quad: 0.0
lab_radiance_scale_factor_gsics: 0.007835
lab_radiance_add_offset_gsics: -0.172166
lab_radiance_quad_gsics: 0.0
radiance_units: mW.cm-2.sr-1.micron-1
coordinates: time Latitude_VIS Longitude_VIS
unlimited dimensions:
current shape = (1, 11264, 11220)
filling on,
'IMG_SWIR_RADIANCE': <class 'netCDF4._netCDF4.Variable'>
float32 IMG_SWIR_RADIANCE(GreyCount)
long_name: Shortwave Infrared Radiance
invert: false
units: mW.cm-2.sr-1.micron-1
_FillValue: 999.0
unlimited dimensions:
current shape = (1024,)
filling on,
'IMG_TIR1': <class 'netCDF4._netCDF4.Variable'>
uint16 IMG_TIR1(time, GeoY, GeoX)
long_name: Thermal Infrared1 Count
invert: true
central_wavelength: 10.8288
bandwidth: 1.0
wavelength_unit: micron
bits_per_pixel: 10
resolution: 4.0
resolution_unit: km
_FillValue: 1023
lab_radiance_scale_factor: 0.00163731
lab_radiance_add_offset: -0.0199179
lab_radiance_quad: 0.0
lab_radiance_scale_factor_gsics: 0.00203982
lab_radiance_add_offset_gsics: -0.118067
lab_radiance_quad_gsics: 0.0
radiance_units: mW.cm-2.sr-1.micron-1
coordinates: time Latitude Longitude
unlimited dimensions:
current shape = (1, 2816, 2805)
filling on,
'IMG_TIR1_RADIANCE': <class 'netCDF4._netCDF4.Variable'>
float32 IMG_TIR1_RADIANCE(GreyCount)
long_name: Thermal Infrared1 Radiance
invert: true
units: mW.cm-2.sr-1.micron-1
_FillValue: 999.0
unlimited dimensions:
current shape = (1024,)
filling on,
'IMG_TIR1_TEMP': <class 'netCDF4._netCDF4.Variable'>
float32 IMG_TIR1_TEMP(GreyCount)
units: K
_FillValue: 999.0
long_name: Thermal Infrared1 Brightness Temperature
invert: true
unlimited dimensions:
current shape = (1024,)
filling on,
'IMG_TIR2': <class 'netCDF4._netCDF4.Variable'>
uint16 IMG_TIR2(time, GeoY, GeoX)
long_name: Thermal Infrared2 Count
invert: true
central_wavelength: 11.9593
bandwidth: 1.0
wavelength_unit: micron
bits_per_pixel: 10
resolution: 4.0
resolution_unit: km
_FillValue: 1023
lab_radiance_scale_factor: 0.00143966
lab_radiance_add_offset: -0.017032
lab_radiance_quad: 0.0
lab_radiance_scale_factor_gsics: 0.00181987
lab_radiance_add_offset_gsics: -0.090434
lab_radiance_quad_gsics: 0.0
radiance_units: mW.cm-2.sr-1.micron-1
coordinates: time Latitude Longitude
unlimited dimensions:
current shape = (1, 2816, 2805)
filling on,
'IMG_TIR2_RADIANCE': <class 'netCDF4._netCDF4.Variable'>
float32 IMG_TIR2_RADIANCE(GreyCount)
_FillValue: 999.0
long_name: Thermal Infrared2 Radiance
invert: true
units: mW.cm-2.sr-1.micron-1
unlimited dimensions:
current shape = (1024,)
filling on,
'IMG_TIR2_TEMP': <class 'netCDF4._netCDF4.Variable'>
float32 IMG_TIR2_TEMP(GreyCount)
long_name: Thermal Infrared2 Brightness Temperature
invert: true
units: K
_FillValue: 999.0
unlimited dimensions:
current shape = (1024,)
filling on,
'IMG_VIS': <class 'netCDF4._netCDF4.Variable'>
uint16 IMG_VIS(time, GeoY2, GeoX2)
long_name: Visible Count
invert: false
central_wavelength: 0.65
bandwidth: 0.25
wavelength_unit: micron
bits_per_pixel: 10
resolution: 1.0
resolution_unit: km
_FillValue: 0
lab_radiance_scale_factor: 0.0630857
lab_radiance_add_offset: -1.85891
lab_radiance_quad: 0.0
lab_radiance_scale_factor_gsics: 0.0630857
lab_radiance_add_offset_gsics: -1.85891
lab_radiance_quad_gsics: 0.0
radiance_units: mW.cm-2.sr-1.micron-1
coordinates: time Latitude_VIS Longitude_VIS
unlimited dimensions:
current shape = (1, 11264, 11220)
filling on,
'IMG_VIS_ALBEDO': <class 'netCDF4._netCDF4.Variable'>
float32 IMG_VIS_ALBEDO(GreyCount)
long_name: Visible Albedo
invert: false
units: %
unlimited dimensions:
current shape = (1024,)
filling on, default _FillValue of 9.969209968386869e+36 used,
'IMG_VIS_RADIANCE': <class 'netCDF4._netCDF4.Variable'>
float32 IMG_VIS_RADIANCE(GreyCount)
long_name: Visible Radiance
invert: false
units: mW.cm-2.sr-1.micron-1
_FillValue: 999.0
unlimited dimensions:
current shape = (1024,)
filling on,
'IMG_WV': <class 'netCDF4._netCDF4.Variable'>
uint16 IMG_WV(time, GeoY1, GeoX1)
long_name: Water Vapor Count
invert: true
central_wavelength: 6.8841
bandwidth: 0.6
wavelength_unit: micron
bits_per_pixel: 10
resolution: 8.0
resolution_unit: km
_FillValue: 1023
lab_radiance_scale_factor: 0.000858417
lab_radiance_add_offset: 0.000380381
lab_radiance_quad: 0.0
lab_radiance_scale_factor_gsics: 0.00125055
lab_radiance_add_offset_gsics: -0.0299274
lab_radiance_quad_gsics: 0.0
radiance_units: mW.cm-2.sr-1.micron-1
coordinates: time Latitude_WV Longitude_WV
unlimited dimensions:
current shape = (1, 1408, 1402)
filling on,
'IMG_WV_RADIANCE': <class 'netCDF4._netCDF4.Variable'>
float32 IMG_WV_RADIANCE(GreyCount)
long_name: Water Vapor Radiance
invert: true
units: mW.cm-2.sr-1.micron-1
_FillValue: 999.0
unlimited dimensions:
current shape = (1024,)
filling on,
'IMG_WV_TEMP': <class 'netCDF4._netCDF4.Variable'>
float32 IMG_WV_TEMP(GreyCount)
long_name: Water Vapor Brightness Temperature
invert: true
units: K
_FillValue: 999.0
unlimited dimensions:
current shape = (1024,)
filling on,
'Latitude': <class 'netCDF4._netCDF4.Variable'>
int16 Latitude(GeoY, GeoX)
long_name: latitude
add_offset: 0.0
scale_factor: 0.01
units: degrees_north
_FillValue: 32767
unlimited dimensions:
current shape = (2816, 2805)
filling on,
'Latitude_VIS': <class 'netCDF4._netCDF4.Variable'>
int32 Latitude_VIS(GeoY2, GeoX2)
long_name: latitude
add_offset: 0.0
scale_factor: 0.001
units: degrees_north
_FillValue: 327670
unlimited dimensions:
current shape = (11264, 11220)
filling on,
'Latitude_WV': <class 'netCDF4._netCDF4.Variable'>
int16 Latitude_WV(GeoY1, GeoX1)
long_name: latitude
add_offset: 0.0
scale_factor: 0.01
units: degrees_north
_FillValue: 32767
unlimited dimensions:
current shape = (1408, 1402)
filling on,
'Longitude': <class 'netCDF4._netCDF4.Variable'>
int16 Longitude(GeoY, GeoX)
long_name: longitude
units: degrees_east
add_offset: 0.0
scale_factor: 0.01
_FillValue: 32767
unlimited dimensions:
current shape = (2816, 2805)
filling on,
'Longitude_VIS': <class 'netCDF4._netCDF4.Variable'>
int32 Longitude_VIS(GeoY2, GeoX2)
add_offset: 0.0
scale_factor: 0.001
long_name: longitude
units: degrees_east
_FillValue: 327670
unlimited dimensions:
current shape = (11264, 11220)
filling on,
'Longitude_WV': <class 'netCDF4._netCDF4.Variable'>
int16 Longitude_WV(GeoY1, GeoX1)
long_name: longitude
add_offset: 0.0
scale_factor: 0.01
units: degrees_east
_FillValue: 32767
unlimited dimensions:
current shape = (1408, 1402)
filling on,
'SCAN_LINE_TIME': <class 'netCDF4._netCDF4.Variable'>
vlen SCAN_LINE_TIME(GeoY1)
long_name: Scan Time for Water Vapor Resolution
vlen data type: <class 'str'>
unlimited dimensions:
current shape = (1408,),
'Sat_Azimuth': <class 'netCDF4._netCDF4.Variable'>
uint16 Sat_Azimuth(time, GeoY, GeoX)
long_name: Satellite Azimuth
add_offset: 0.0
scale_factor: 0.01
units: degree
_FillValue: 65535
coordinates: time Latitude Longitude
unlimited dimensions:
current shape = (1, 2816, 2805)
filling on,
'Sat_Elevation': <class 'netCDF4._netCDF4.Variable'>
int16 Sat_Elevation(time, GeoY, GeoX)
long_name: Satellite Elevation
add_offset: 0.0
units: degree
scale_factor: 0.01
_FillValue: 32767
coordinates: time Latitude Longitude
unlimited dimensions:
current shape = (1, 2816, 2805)
filling on,
'Sun_Azimuth': <class 'netCDF4._netCDF4.Variable'>
uint16 Sun_Azimuth(time, GeoY, GeoX)
long_name: Sun Azimuth
add_offset: 0.0
units: degree
scale_factor: 0.01
_FillValue: 65535
coordinates: time Latitude Longitude
unlimited dimensions:
current shape = (1, 2816, 2805)
filling on,
'Sun_Elevation': <class 'netCDF4._netCDF4.Variable'>
int16 Sun_Elevation(time, GeoY, GeoX)
long_name: Sun Elevation
add_offset: 0.0
scale_factor: 0.01
units: degree
_FillValue: 32767
coordinates: time Latitude Longitude
unlimited dimensions:
current shape = (1, 2816, 2805)
filling on,
'time': <class 'netCDF4._netCDF4.Variable'>
float64 time(time)
units: minutes since 2000-01-01 00:00:00
unlimited dimensions:
current shape = (1,)
filling off}
data.variables.keys()
Out: dict_keys(['GeoX', 'GeoX1', 'GeoX2', 'GeoY', 'GeoY1', 'GeoY2', 'GreyCount', 'IMG_MIR', 'IMG_MIR_RADIANCE', 'IMG_MIR_TEMP', 'IMG_SWIR', 'IMG_SWIR_RADIANCE', 'IMG_TIR1', 'IMG_TIR1_RADIANCE', 'IMG_TIR1_TEMP', 'IMG_TIR2', 'IMG_TIR2_RADIANCE', 'IMG_TIR2_TEMP', 'IMG_VIS', 'IMG_VIS_ALBEDO', 'IMG_VIS_RADIANCE', 'IMG_WV', 'IMG_WV_RADIANCE', 'IMG_WV_TEMP', 'Latitude', 'Latitude_VIS', 'Latitude_WV', 'Longitude', 'Longitude_VIS', 'Longitude_WV', 'SCAN_LINE_TIME', 'Sat_Azimuth', 'Sat_Elevation', 'Sun_Azimuth', 'Sun_Elevation', 'time'])
lats = data.variables['Latitude'][:]
lons = data.variables['Longitude'][:]
tir = data.variables['IMG_TIR1'][:]
lats.shape
Out:(2816, 2805)
lons.shape
Out:(2816, 2805)
tir.shape
Out:(1, 2816, 2805)
` I was not so sure about exact lat and lon that I want to plot so here's two of them. Although I think the 1st one is the correct one`
#mp = Basemap(projection = 'merc',
# llcrnrlon = 68.766521,
# llcrnrlat = 3.323179,
# urcrnrlon = 100.202163,
# urcrnrlat = 37.315495,
# resolution = 'i')
mp = Basemap(projection = 'merc',
llcrnrlon = 0.8432964,
llcrnrlat = -81.04153,
urcrnrlon = 163.15671,
urcrnrlat = 81.04153,
resolution = 'i')
` I'm getting an error in this line dk what it is`
lon, lat = np.meshgrid(lons, lats)
Out: ---------------------------------------------------------------------------
MemoryError Traceback (most recent call last)
~\AppData\Local\Temp/ipykernel_12060/1726626783.py in <module>
----> 1 lon, lat = np.meshgrid(lons, lats)
~\anaconda3\lib\site-packages\numpy\core\overrides.py in meshgrid(*args, **kwargs)
~\anaconda3\lib\site-packages\numpy\lib\function_base.py in meshgrid(copy, sparse, indexing, *xi)
4947
4948 if copy:
-> 4949 output = [x.copy() for x in output]
4950
4951 return output
~\anaconda3\lib\site-packages\numpy\lib\function_base.py in <listcomp>(.0)
4947
4948 if copy:
-> 4949 output = [x.copy() for x in output]
4950
4951 return output
~\anaconda3\lib\site-packages\numpy\ma\core.py in wrapped_method(self, *args, **params)
2576 """
2577 def wrapped_method(self, *args, **params):
-> 2578 result = getattr(self._data, funcname)(*args, **params)
2579 result = result.view(type(self))
2580 result._update_from(self)
MemoryError: Unable to allocate 227. TiB for an array with shape (7898880, 7898880) and data type float32
x,y = mp(lon, lat)
Out:---------------------------------------------------------------------------
NameError Traceback (most recent call last)
~\AppData\Local\Temp/ipykernel_6908/3885003198.py in <module>
----> 1 x,y = mp(lon, lat)
NameError: name 'lon' is not defined
c_scheme = mp.pcolor(x, y, np.squeeze(tir1[0, :, :]), cmap = 'jet')
Out: ---------------------------------------------------------------------------
NameError Traceback (most recent call last)
~\AppData\Local\Temp/ipykernel_6908/254152092.py in <module>
----> 1 c_scheme = mp.pcolor(x, y, np.squeeze(tir1[0, :, :]), cmap = 'jet')
NameError: name 'x' is not defined
更新 2022-03-08：Pseudo-code 在我原来的 post 中使用 plt.contourf() 绘制数据。这不适用于保存在 'IMG_TIR1' 数据集中的栅格图像数据。相反，必须使用 plt.imshow()。有关完整过程，请参阅 2022 年 3 月 7 日编辑 post 的详细答案。
所有'IMG_name'数据集都是栅格图像数据。我修改了代码以显示一个简单示例，该示例使用 h5py 读取 HDF5 文件，然后使用 matplotlib 和 cartopy 在地球静止（原始卫星）投影上绘制图像数据。我根据 Basemap 网站上的评论选择了 cartopy：“弃用通知：不推荐使用 Basemap 以支持 Cartopy 项目。”我确定可以修改它以使用 Basemap 绘制和 NetCDF 读取 HDF5 文件。
更新代码如下：
import h5py
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
# First get data from HDF5 file with h5py:
fn = '3DIMG_30MAR2018_0000_L1B_STD.h5' #filename (the ".h5" file)
with h5py.File(fn) as data:
tir1 = data['IMG_TIR1'][:].reshape(lats.shape)
# retrieve extent of plot from file attributes:
left_lon = data.attrs['left_longitude'][0]
right_lon = data.attrs['right_longitude'][0]
lower_lat = data.attrs['lower_latitude'][0]
upper_lat = data.attrs['upper_latitude'][0]
sat_long = data.attrs['Nominal_Central_Point_Coordinates(degrees)_Latitude_Longitude'][1]
# Define extent of plot (in degrees) as:
# extent=[longitude_top_left,longitude_top_right,
# latitude_bottom_left,latitude_top_left]
# coordinates must be transformed to map coordinate system before using
img_extent_deg = (left_lon, right_lon, lower_lat, upper_lat)
# Create Geostationary plot with cartopy and matplotlib
map_proj = ccrs.Geostationary(central_longitude=sat_long)
ax = plt.axes(projection=map_proj)
ax.coastlines(color='white')
map_extend_geos = ax.get_extent(crs=map_proj)
plt.imshow(tir1, extent=map_extend_geos)
plt.colorbar(im1)
虽然我找不到数据，但我使用了 The HDF Group 的 HDFView 从数据中生成图像。这是 IMG_TIR1。
此时，我需要有关文件架构的更多信息来“查找”此数据。 MOSDAC 网站是否有任何描述文件架构的内容，或者是否有一些显示如何读取数据的代码？
@Black Viking，回答你的问题比我最初想象的要复杂。我从评论中的代码开始。经过一番挖掘，我确定每个 'IMG_xxx' 数据集中的数据是原始栅格图像（扫描），相关经度和纬度数据集中的值是每个像素的 (lon,lat) 位置。它们不形成规则的 MxN 网格，因此不能与 matplotlib 的 contourf() 函数一起使用。所以，我之前的回答是不正确的。 （我会删除其中的大部分内容，以免让未来的读者感到困惑。）
当我从事这项工作时，我发现在具有不同坐标系的地图上绘制卫星光栅图像并不是一项简单的任务（从地球静止投影到墨卡托投影）。您必须做的事情：
- 定义图像的原点坐标参考系。
- 在地球静止系统中绘制图像时，要显示轴，您需要将
extent= 添加到 plt.imshow()。
- 在不同坐标系中绘制图像时，必须使用
'transform=' 和 plt.imshow() 将图像转换到新坐标系。
我将逐步 post 代码，以展示如何从对地静止系统中的简单栅格图转到聚焦于印度地区的墨卡托投影中的同一图像。请阅读到最后 - 如果某些内容不起作用，您可能需要一些其他信息。
情节 0：起点来自您的评论（略有修改）。这仅绘制原始 IMG_TIR1 图像数据。
import h5py
import matplotlib.pyplot as plt
fn = '3DIMG_30MAR2018_0000_L1B_STD.h5' #filename (the ".h5" file)
with h5py.File(fn) as f:
img_arr = f['IMG_TIR1'][0,:,:]
fig = plt.subplots(figsize=(10,10))
plt.title('plot raw IMG_TIR1 image data (only)')
im = plt.imshow(img_arr)
plt.colorbar(im)
Plot 0（带海岸线，图像经过修剪）：这是将海岸线添加到原始 IMG_TIR1 图像数据的小修改。这是我验证墨卡托投影在以下图像中正确定向的参考。请注意它如何使用 extent=img_extent_sat 以便正确显示轴。这很棘手，因为地球静止坐标是 而不是 （经度、纬度）。 ax.get_extent(projection=) returns 该投影的地图范围。我不得不稍微增加它们以使图像“适合”轴。
import h5py
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
fn = '3DIMG_30MAR2018_0000_L1B_STD.h5' #filename (the ".h5" file)
with h5py.File(fn) as f:
img_arr = f['IMG_TIR1'][0,:,:]
map_proj = ccrs.Geostationary(central_longitude=82.0)
fig = plt.figure(figsize=(10,10))
ax = plt.axes(projection=map_proj)
ax.coastlines(color='white')
# Image extent in Geostationary coordinates:
img_extent_sat = ax.get_extent(crs=map_proj)
img_extent_sat = [1.04*x for x in img_extent_sat]
im = plt.imshow(img_arr, extent=img_extent_sat)
plt.colorbar(im)
plt.title('plot raw IMG_TIR1 data + Geostationary coastlines')
这是在墨卡托投影中创建绘图的代码。为避免重复，我 post 为每个步骤编辑了代码段。使用第一个（“基础”）部分并为您想要的情节添加部分。
基本代码：从 HDF5 文件中获取图像数据和属性以及卫星属性。我定义了 image = 'IMG_TIR1' 以便您可以轻松修改以获取其他 IMG_name 数据。仅当您要从图像数据（计数值）计算辐射率时才需要 Radiance 属性。最后一行创建了一个屏蔽数组，屏蔽了角落里的_FillValue数据并简化了图像绘制。
import h5py
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
fn = '3DIMG_30MAR2018_0000_L1B_STD.h5' #filename (the ".h5" file)
with h5py.File(fn) as f:
# retrieve image data:
image = 'IMG_TIR1'
img_arr = f[image][0,:,:]
# get _FillValue for data masking
img_arr_fill = f[image].attrs['_FillValue'][0]
# retrieve extent of plot from file attributes:
left_lon = f.attrs['left_longitude'][0]
right_lon = f.attrs['right_longitude'][0]
lower_lat = f.attrs['lower_latitude'][0]
upper_lat = f.attrs['upper_latitude'][0]
sat_long = f.attrs['Nominal_Central_Point_Coordinates(degrees)_Latitude_Longitude'][1]
sat_hght = f.attrs['Observed_Altitude(km)'][0] * 1000.0 # (for meters)
# retrieve attributes to calculate radiance from count:
Sensor_Name = f.attrs['Sensor_Name'].decode('utf-8')
img_fill = f[image].attrs['_FillValue'][0]
img_inv = f[image].attrs['invert'].decode('utf-8')
img_lrquad = f[image].attrs['lab_radiance_quad'][0]
img_lrscale = f[image].attrs['lab_radiance_scale_factor'][0]
img_lroff = f[image].attrs['lab_radiance_add_offset'][0]
print('Done reading HDF5 file')
## Use np.ma.masked_equal with integer values to
## mask '_FillValue' data in corners:
img_arr_m = np.ma.masked_equal(img_arr, img_arr_fill, copy=True)
绘图 1：将蒙版图像数据投影到全墨卡托投影。请注意 ax1.imshow() 如何使用 transform=data_crs，其中 data_crs 引用图像数据的地球静止投影。注意：如果您绘制 img_arr 而不是 img_arr_m，您将在图像周围看到黄色的“_FillValue”“光环”。此外，我在 ax# 和 im# 中添加了数字，以简化将以下 3 个图形绘制为 1 个图形的子图。
# Plot 1 uses Mercator projection
print('Start on Plot1')
map_proj = ccrs.Mercator(central_longitude=sat_long,
min_latitude=lower_lat, max_latitude=upper_lat)
data_crs = ccrs.Geostationary(central_longitude=sat_long,
satellite_height=sat_hght)
plt.figure(figsize=(10,10))
ax1 = plt.axes(projection=map_proj)
ax1.coastlines()
#ax1.add_feature(cfeature.BORDERS, edgecolor='white', linewidth=0.5)
ax1.gridlines(color='black', alpha=0.5, linestyle='--', linewidth=0.75, draw_labels=True)
map_proj_text = f'{str(type(map_proj)).split(".")[-1][:-2]}'
data_crs_text = f'{str(type(data_crs)).split(".")[-1][:-2]}'
plt.title(f'Plot1: Projection: {map_proj_text}\n' + \
f'Data Transform: {data_crs_text}\n' + \
f'\nRaster Data: {image} (masked)')
print('plotting data for Plot1 image')
im1 = ax1.imshow(img_arr_m, origin='upper', transform=data_crs)
plt.colorbar(im1)
print('done with Plot1')
绘图 2：将蒙版图像数据投影到经过修剪的墨卡托投影。地图大小由 ax2.set_extent() 中的 map_extent_merc 参数定义。请注意，墨卡托坐标不是度数，因此我创建了一个函数来转换范围值 (transform_extent_pts())。 （对于您的数据，它从 PlateCarree 中的 map_extent_deg 变为墨卡托中的 map_extent_merc）。尝试 map_extent_deg 的不同值以查看地图如何变化。
# Plot 2 on Mercator projection
print('Start on Plot2')
map_proj = ccrs.Mercator(central_longitude=sat_long,
min_latitude=lower_lat, max_latitude=upper_lat)
data_crs = ccrs.Geostationary(central_longitude=sat_long,
satellite_height=sat_hght)
plt.figure(figsize=(10,10))
ax2 = plt.axes(projection=map_proj)
ax2.coastlines()
ax2.add_feature(cfeature.BORDERS, edgecolor='white', linewidth=0.5)
ax2.gridlines(color='black', alpha=0.5, linestyle='--', linewidth=0.75, draw_labels=True)
# Focus map on Indian subcontinent:
deg_crs = ccrs.PlateCarree()
map_extent_deg = (15., 150., -67., 67.) # lon/lat focused on image
#map_extent_deg = (60.0, 97.0, -10.0, 37.5) # India
map_extent_merc = transform_extent_pts(map_extent_deg, map_proj, deg_crs)
ax2.set_extent(map_extent_merc, map_proj)
map_proj_text = f'{str(type(map_proj)).split(".")[-1][:-2]}'
data_crs_text = f'{str(type(data_crs)).split(".")[-1][:-2]}'
plt.title(f'Plot2: Projection: {map_proj_text}\n' + \
f'Data Transform: {data_crs_text}\n' + \
f"\nRaster Data: {image} (masked) ")
print('plotting data for Plot2 image')
im2 = plt.imshow(img_arr_m, origin='upper', transform=data_crs)
plt.colorbar(im2)
print('done with Plot2')
图 2 中使用的转换函数:
def transform_extent_pts(extent_pts, map_proj, pt_crs):
xul, yul = map_proj.transform_point(
x = extent_pts[0],
y = extent_pts[3],
src_crs = pt_crs)
xlr, ylr = map_proj.transform_point(
x = extent_pts[1],
y = extent_pts[2],
src_crs = pt_crs)
return [xul, xlr, ylr, yul]
绘图 3：计算并显示辐射值。以前的图显示原始图像数据（即“计数”）。您提供的论文具有将“计数”转换为辐射率的方程式。我写了一个函数来做到这一点 (calc_radiance())。这只需要从上面进行一些小的更改，如下所述。
首先，添加辐射功能。
def transform_extent_pts(extent_pts, map_proj, pt_crs):
xul, yul = map_proj.transform_point(
x = extent_pts[0],
y = extent_pts[3],
src_crs = pt_crs)
xlr, ylr = map_proj.transform_point(
x = extent_pts[1],
y = extent_pts[2],
src_crs = pt_crs)
return [xul, xlr, ylr, yul]
在之前添加以下行：print('Start on Plot2') 以计算辐射值。此数据用于图 3：
# Plot 3: Plot radiance on Mercator projection
print('Calculate radiance for Plot3')
radiance = calc_radiance(img_arr_m, img_lrquad, img_lrscale, img_lroff,
invert=img_inv, Sensor_Name=Sensor_Name)
接下来，将轴引用从 ax2 修改为 ax3。
最后修改最后3行：
plt.title(f'Plot3: Projection: {map_proj_text}\n' + \
f'Data Transform: {map_proj_text}\n' + \
f'\nRadiance from: {image} (masked)')
print('plot image for Plot3')
im3 = ax3.imshow(radiance, origin='upper', transform=data_crs)
plt.colorbar(im3)
print('done with Plot3')
图 4：调整图像数据以绘制云阈值。这类似于图 3。您可以使用 np.where() 函数创建一个使用阈值（小于、大于或两者）的新图像数组。此外，我添加了 vmin=, vmax= 参数以将颜色条范围与您的绘图相匹配。只需要做一些小的改动。为简洁起见，下面仅说明更改。
# Plot 4: Plot cloud thresholds on Mercator projection
print('Calculate cloud mask for Plot4')
# Left Plot: set values <=700 to 0
cloud1 = np.where(img_arr_m <= 700, 0, img_arr_m)
# Right Plot: set values <=700 to 0, >700 to 1
cloud2 = np.where(img_arr_m <= 700, 0, 1)
.....
.....
print('plot image for Plot4')
# to create plot on left:
im4 = ax4.imshow(cloud1, vmin=0, vmax=1000, origin='upper', transform=data_crs)
# to create plot on right:
im4 = ax4.imshow(cloud2, vmin=0, vmax=1, origin='upper', transform=data_crs)
plt.colorbar(im4)
背景和参考信息：
注意：我在使用早期版本转换 GeoStationary 图像时遇到了问题。我不得不使用带有 Proj 8.2.1 的 cartopy v0.20.2 来获取这些图像。如果有问题，运行这个例子来测试：Reprojecting images from a Geostationary projection
另外，这里有 2 个有用的博客 post 的链接（均包含 Python 代码）
- Blog tutorial 附有在 cartopy 中绘制整个地球图像的详细指南
- Another blog tutorial 展示了如何使用 cartopy 从 WMS 绘制卫星图像。
最后，Whosebug 对此主题有一些很好的问答。使用 tag: [cartopy] 'satellite' 搜索以找到一些好文章。
对于以后可能阅读这些内容的人：@kcw78 提供的代码是正确的，您只需要对图像数据的范围做一些小的改动。对于我来说，我将印度地区的 map_extent_deg = (15., 150., -67., 67.) 更改为 map_extent_deg = (65., 100., 0., 40.)。像这样你可能需要根据你的需要改变它。
对于以后可能阅读这些内容的人：@kcw78 提供的代码是正确的，您只需要对图像数据的范围做一些小的改动。对我来说，我将它从 map_extent_deg = (15., 150., -67., 67.) 更改为 map_extent_deg = (65., 100., 0., 40.)。
获得相同结果的另一个代码：
import cv2
import numpy as np
import matplotlib.pyplot as plt
fn = '3DIMG_30MAR2018_0000_L1B_STD.h5' #filename (the ".h5" file)
hf = h5py.File(fn, 'r')
hf.keys()
data = hf.get('IMG_TIR1')
data = data[0,:,:]
data=np.where(data>1022, np.nan, data)
Lat= hf.get('Latitude')
Lat = Lat[:,:]
Lat=Lat*0.0099999998
Lat=np.where(Lat>312, np.nan, Lat)
Long= hf.get('Longitude')
Long = Long[:,:]
Long=Long*0.0099999998
Long=np.where(Long>312, np.nan, Long)
X=Long[427:1421,804:1902]
Y=Lat[427:1421,804:1902]
Z=data[427:1421,804:1902]
fig = plt.figure(figsize=(10,10))
plt.pcolor(X, Y, Z)
plt.xlim(60, 100)
plt.ylim(0, 40)
fig = plt.figure(figsize=(10,10))
plt.pcolor(X, Y, Z)
plt.colorbar(label="Radiance", orientation="vertical")
plt.xlim(60, 100)
plt.ylim(0, 40)
#Threshold
Zn=np.where(Z>810, 1, 0)
fig = plt.figure(figsize=(10,10))
plt.pcolor(X, Y, Zn, cmap='jet', vmin=0, vmax=1)
plt.colorbar(label="Radiance", orientation="vertical")
plt.xlim(60, 100)
plt.ylim(0, 40)
我从 MOSDAC 网站获得了一个名为 3DIMG_30MAR2018_0000_L1B_STD.h5 的 HDF5 文件，我想读取这个文件并在地图中可视化文件中的 TIR 和 VIS 数据集。由于我是编程新手，所以我不太清楚如何使用纬度和经度值对全球数据的特定位置执行此操作。我的问题陈述是显示印度地区的数据。
到目前为止我完成的代码：
#Importing libraries`
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import h5py
import os
import netCDF4 as nc
import seaborn as sns
import geopandas as gpd
#Reading and listing keys present in HDF
fn = '3DIMG_30MAR2018_0000_L1B_STD.h5' #filename (the ".h5" file)
f = h5py.File(fn)
list(f.keys())
#extract data from the HDF5 file
mir = f['IMG_MIR'][:] #middle infra-red
swir = f['IMG_SWIR'][:] #shortwave IR
tir1 = f['IMG_TIR1'][:] #Thermal IR1
tir2 = f['IMG_TIR2'][:] #Thermal IR2
vis = f['IMG_VIS'][:] #Visible count
wv = f['IMG_WV'][:] #Water vapor count
lat = f['Latitude'][:]
lon = f['Longitude'][:]
lat_vis = f['Latitude_VIS'][:]
lon_vis = f['Longitude_VIS'][:]
lat_wv = f['Latitude_WV'][:]
lon_wv = f['Longitude_WV'][:]
vis = f['IMG_MIR'][0,:,:]
fig,ax = plt.subplots(figsize=(10,10))
im1 = plt.imshow(vis)
plt.colorbar(im1)
tir1d = f['/IMG_TIR1'];
lat = f['/Latitude'];
lon = f['/Longitude'];
Lat = np.linspace(81.041527, -81.041527, 2805)
Lon = np.linspace(0.84329641, 163.15671, 2816)
shp = gpd.read_file(r'C:\Users\offic\shape\IND_adm0.shp')
%lat=0:gs:60;%INDIA
%lon=60:gs:100;
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
data = Dataset(r'C:\Users\officDIMG_30MAR2018_0000_L1B_STD.h5')
data
Out: <class 'netCDF4._netCDF4.Dataset'>
root group (NETCDF4 data model, file format HDF5):
conventions: CF-1.6
title: 3DIMG_30MAR2018_0000_L1B
institute: BES,SAC/ISRO,Ahmedabad,INDIA.
source: BES,SAC/ISRO,Ahmedabad,INDIA.
Unique_Id: 3DIMG_30MAR2018_0000
Satellite_Name: INSAT-3D
Sensor_Id: IMG
Sensor_Name: IMAGER
HDF_Product_File_Name: 3DIMG_30MAR2018_0000_L1B_STD.h5
Output_Format: hdf5-1.8.8
Station_Id: BES
Ground_Station: BES,SAC/ISRO,Ahmedabad,INDIA.
Product_Type: STANDARD(FULL DISK)
Processing_Level: L1B
Imaging_Mode: FULL FRAME
Acquisition_Date: 30MAR2018
Acquisition_Time_in_GMT: 0000
Acquisition_Start_Time: 30-MAR-2018T00:00:08
Acquisition_End_Time: 30-MAR-2018T00:26:58
Product_Creation_Time: 2018-03-30T06:03:39
Radiometric_Calibration_Type: LAB CALIBRATED
Nominal_Altitude(km): 36000.0
Observed_Altitude(km): 35787.8215
Field_of_View(degrees): 17.973925
Nominal_Central_Point_Coordinates(degrees)_Latitude_Longitude: [ 0. 82.]
Attitude_Source: STAR
Sun_Azimuth(Degrees): 86.314133
Sat_Azimuth(Degrees): 286.793152
Sat_Elevation(Degrees): 89.813721
Sun_Elevation(Degrees): -6.010917
FastScan_Linearity_Enabled: no
SlowScan_Linearity_Enabled: no
MCD_FS_Enabled: no
MCD_SS_Enabled: no
Yaw_Flip_Flag: Y
Software_Version: 1.0
TIR1_Gain_Mode: 2
TIR2_Gain_Mode: 2
MIR_Gain_Mode: 2
WV_Gain_Mode: 3
VIS_Gain_Mode: 2
SWIR_Gain_Mode: 3
TIR1_Acquisition_Mode: MAIN
TIR2_Acquisition_Mode: MAIN
MIR_Acquisition_Mode: MAIN
WV_Acquisition_Mode: MAIN
VIS_Acquisition_Mode: MAIN
SWIR_Acquisition_Mode: MAIN
left_longitude: 0.8432964
right_longitude: 163.15671
upper_latitude: 81.04153
lower_latitude: -81.04153
Datum: WGS84
Ellipsoid: WGS84
dimensions(sizes): GeoX(2805), GeoX1(1402), GeoX2(11220), GeoY(2816), GeoY1(1408), GeoY2(11264), GreyCount(1024), time(1)
variables(dimensions): int32 GeoX(GeoX), int32 GeoX1(GeoX1), int32 GeoX2(GeoX2), int32 GeoY(GeoY), int32 GeoY1(GeoY1), int32 GeoY2(GeoY2), int32 GreyCount(GreyCount), uint16 IMG_MIR(time, GeoY, GeoX), float32 IMG_MIR_RADIANCE(GreyCount), float32 IMG_MIR_TEMP(GreyCount), uint16 IMG_SWIR(time, GeoY2, GeoX2), float32 IMG_SWIR_RADIANCE(GreyCount), uint16 IMG_TIR1(time, GeoY, GeoX), float32 IMG_TIR1_RADIANCE(GreyCount), float32 IMG_TIR1_TEMP(GreyCount), uint16 IMG_TIR2(time, GeoY, GeoX), float32 IMG_TIR2_RADIANCE(GreyCount), float32 IMG_TIR2_TEMP(GreyCount), uint16 IMG_VIS(time, GeoY2, GeoX2), float32 IMG_VIS_ALBEDO(GreyCount), float32 IMG_VIS_RADIANCE(GreyCount), uint16 IMG_WV(time, GeoY1, GeoX1), float32 IMG_WV_RADIANCE(GreyCount), float32 IMG_WV_TEMP(GreyCount), int16 Latitude(GeoY, GeoX), int32 Latitude_VIS(GeoY2, GeoX2), int16 Latitude_WV(GeoY1, GeoX1), int16 Longitude(GeoY, GeoX), int32 Longitude_VIS(GeoY2, GeoX2), int16 Longitude_WV(GeoY1, GeoX1), <class 'str'> SCAN_LINE_TIME(GeoY1), uint16 Sat_Azimuth(time, GeoY, GeoX), int16 Sat_Elevation(time, GeoY, GeoX), uint16 Sun_Azimuth(time, GeoY, GeoX), int16 Sun_Elevation(time, GeoY, GeoX), float64 time(time)
groups:
data.variables
Out: {'GeoX': <class 'netCDF4._netCDF4.Variable'>
int32 GeoX(GeoX)
unlimited dimensions:
current shape = (2805,)
filling off,
'GeoX1': <class 'netCDF4._netCDF4.Variable'>
int32 GeoX1(GeoX1)
unlimited dimensions:
current shape = (1402,)
filling off,
'GeoX2': <class 'netCDF4._netCDF4.Variable'>
int32 GeoX2(GeoX2)
unlimited dimensions:
current shape = (11220,)
filling off,
'GeoY': <class 'netCDF4._netCDF4.Variable'>
int32 GeoY(GeoY)
unlimited dimensions:
current shape = (2816,)
filling off,
'GeoY1': <class 'netCDF4._netCDF4.Variable'>
int32 GeoY1(GeoY1)
unlimited dimensions:
current shape = (1408,)
filling off,
'GeoY2': <class 'netCDF4._netCDF4.Variable'>
int32 GeoY2(GeoY2)
unlimited dimensions:
current shape = (11264,)
filling off,
'GreyCount': <class 'netCDF4._netCDF4.Variable'>
int32 GreyCount(GreyCount)
unlimited dimensions:
current shape = (1024,)
filling off,
'IMG_MIR': <class 'netCDF4._netCDF4.Variable'>
uint16 IMG_MIR(time, GeoY, GeoX)
long_name: Middle Infrared Count
invert: true
central_wavelength: 3.9313
bandwidth: 0.2
wavelength_unit: micron
bits_per_pixel: 10
resolution: 4.0
resolution_unit: km
_FillValue: 1023
lab_radiance_scale_factor: 0.000283937
lab_radiance_add_offset: -0.00401529
lab_radiance_quad: 0.0
lab_radiance_scale_factor_gsics: 0.000366208
lab_radiance_add_offset_gsics: -0.00833466
lab_radiance_quad_gsics: 0.0
radiance_units: mW.cm-2.sr-1.micron-1
coordinates: time Latitude Longitude
unlimited dimensions:
current shape = (1, 2816, 2805)
filling on,
'IMG_MIR_RADIANCE': <class 'netCDF4._netCDF4.Variable'>
float32 IMG_MIR_RADIANCE(GreyCount)
long_name: Middle Infrared Radiance
invert: true
units: mW.cm-2.sr-1.micron-1
_FillValue: 999.0
unlimited dimensions:
current shape = (1024,)
filling on,
'IMG_MIR_TEMP': <class 'netCDF4._netCDF4.Variable'>
float32 IMG_MIR_TEMP(GreyCount)
long_name: Middle Infrared Brightness Temperature
invert: true
units: K
_FillValue: 999.0
unlimited dimensions:
current shape = (1024,)
filling on,
'IMG_SWIR': <class 'netCDF4._netCDF4.Variable'>
uint16 IMG_SWIR(time, GeoY2, GeoX2)
long_name: Shortwave Infrared Count
invert: false
central_wavelength: 1.625
bandwidth: 0.15
wavelength_unit: micron
bits_per_pixel: 10
resolution: 1.0
resolution_unit: km
_FillValue: 0
lab_radiance_scale_factor: 0.007835
lab_radiance_add_offset: -0.172166
lab_radiance_quad: 0.0
lab_radiance_scale_factor_gsics: 0.007835
lab_radiance_add_offset_gsics: -0.172166
lab_radiance_quad_gsics: 0.0
radiance_units: mW.cm-2.sr-1.micron-1
coordinates: time Latitude_VIS Longitude_VIS
unlimited dimensions:
current shape = (1, 11264, 11220)
filling on,
'IMG_SWIR_RADIANCE': <class 'netCDF4._netCDF4.Variable'>
float32 IMG_SWIR_RADIANCE(GreyCount)
long_name: Shortwave Infrared Radiance
invert: false
units: mW.cm-2.sr-1.micron-1
_FillValue: 999.0
unlimited dimensions:
current shape = (1024,)
filling on,
'IMG_TIR1': <class 'netCDF4._netCDF4.Variable'>
uint16 IMG_TIR1(time, GeoY, GeoX)
long_name: Thermal Infrared1 Count
invert: true
central_wavelength: 10.8288
bandwidth: 1.0
wavelength_unit: micron
bits_per_pixel: 10
resolution: 4.0
resolution_unit: km
_FillValue: 1023
lab_radiance_scale_factor: 0.00163731
lab_radiance_add_offset: -0.0199179
lab_radiance_quad: 0.0
lab_radiance_scale_factor_gsics: 0.00203982
lab_radiance_add_offset_gsics: -0.118067
lab_radiance_quad_gsics: 0.0
radiance_units: mW.cm-2.sr-1.micron-1
coordinates: time Latitude Longitude
unlimited dimensions:
current shape = (1, 2816, 2805)
filling on,
'IMG_TIR1_RADIANCE': <class 'netCDF4._netCDF4.Variable'>
float32 IMG_TIR1_RADIANCE(GreyCount)
long_name: Thermal Infrared1 Radiance
invert: true
units: mW.cm-2.sr-1.micron-1
_FillValue: 999.0
unlimited dimensions:
current shape = (1024,)
filling on,
'IMG_TIR1_TEMP': <class 'netCDF4._netCDF4.Variable'>
float32 IMG_TIR1_TEMP(GreyCount)
units: K
_FillValue: 999.0
long_name: Thermal Infrared1 Brightness Temperature
invert: true
unlimited dimensions:
current shape = (1024,)
filling on,
'IMG_TIR2': <class 'netCDF4._netCDF4.Variable'>
uint16 IMG_TIR2(time, GeoY, GeoX)
long_name: Thermal Infrared2 Count
invert: true
central_wavelength: 11.9593
bandwidth: 1.0
wavelength_unit: micron
bits_per_pixel: 10
resolution: 4.0
resolution_unit: km
_FillValue: 1023
lab_radiance_scale_factor: 0.00143966
lab_radiance_add_offset: -0.017032
lab_radiance_quad: 0.0
lab_radiance_scale_factor_gsics: 0.00181987
lab_radiance_add_offset_gsics: -0.090434
lab_radiance_quad_gsics: 0.0
radiance_units: mW.cm-2.sr-1.micron-1
coordinates: time Latitude Longitude
unlimited dimensions:
current shape = (1, 2816, 2805)
filling on,
'IMG_TIR2_RADIANCE': <class 'netCDF4._netCDF4.Variable'>
float32 IMG_TIR2_RADIANCE(GreyCount)
_FillValue: 999.0
long_name: Thermal Infrared2 Radiance
invert: true
units: mW.cm-2.sr-1.micron-1
unlimited dimensions:
current shape = (1024,)
filling on,
'IMG_TIR2_TEMP': <class 'netCDF4._netCDF4.Variable'>
float32 IMG_TIR2_TEMP(GreyCount)
long_name: Thermal Infrared2 Brightness Temperature
invert: true
units: K
_FillValue: 999.0
unlimited dimensions:
current shape = (1024,)
filling on,
'IMG_VIS': <class 'netCDF4._netCDF4.Variable'>
uint16 IMG_VIS(time, GeoY2, GeoX2)
long_name: Visible Count
invert: false
central_wavelength: 0.65
bandwidth: 0.25
wavelength_unit: micron
bits_per_pixel: 10
resolution: 1.0
resolution_unit: km
_FillValue: 0
lab_radiance_scale_factor: 0.0630857
lab_radiance_add_offset: -1.85891
lab_radiance_quad: 0.0
lab_radiance_scale_factor_gsics: 0.0630857
lab_radiance_add_offset_gsics: -1.85891
lab_radiance_quad_gsics: 0.0
radiance_units: mW.cm-2.sr-1.micron-1
coordinates: time Latitude_VIS Longitude_VIS
unlimited dimensions:
current shape = (1, 11264, 11220)
filling on,
'IMG_VIS_ALBEDO': <class 'netCDF4._netCDF4.Variable'>
float32 IMG_VIS_ALBEDO(GreyCount)
long_name: Visible Albedo
invert: false
units: %
unlimited dimensions:
current shape = (1024,)
filling on, default _FillValue of 9.969209968386869e+36 used,
'IMG_VIS_RADIANCE': <class 'netCDF4._netCDF4.Variable'>
float32 IMG_VIS_RADIANCE(GreyCount)
long_name: Visible Radiance
invert: false
units: mW.cm-2.sr-1.micron-1
_FillValue: 999.0
unlimited dimensions:
current shape = (1024,)
filling on,
'IMG_WV': <class 'netCDF4._netCDF4.Variable'>
uint16 IMG_WV(time, GeoY1, GeoX1)
long_name: Water Vapor Count
invert: true
central_wavelength: 6.8841
bandwidth: 0.6
wavelength_unit: micron
bits_per_pixel: 10
resolution: 8.0
resolution_unit: km
_FillValue: 1023
lab_radiance_scale_factor: 0.000858417
lab_radiance_add_offset: 0.000380381
lab_radiance_quad: 0.0
lab_radiance_scale_factor_gsics: 0.00125055
lab_radiance_add_offset_gsics: -0.0299274
lab_radiance_quad_gsics: 0.0
radiance_units: mW.cm-2.sr-1.micron-1
coordinates: time Latitude_WV Longitude_WV
unlimited dimensions:
current shape = (1, 1408, 1402)
filling on,
'IMG_WV_RADIANCE': <class 'netCDF4._netCDF4.Variable'>
float32 IMG_WV_RADIANCE(GreyCount)
long_name: Water Vapor Radiance
invert: true
units: mW.cm-2.sr-1.micron-1
_FillValue: 999.0
unlimited dimensions:
current shape = (1024,)
filling on,
'IMG_WV_TEMP': <class 'netCDF4._netCDF4.Variable'>
float32 IMG_WV_TEMP(GreyCount)
long_name: Water Vapor Brightness Temperature
invert: true
units: K
_FillValue: 999.0
unlimited dimensions:
current shape = (1024,)
filling on,
'Latitude': <class 'netCDF4._netCDF4.Variable'>
int16 Latitude(GeoY, GeoX)
long_name: latitude
add_offset: 0.0
scale_factor: 0.01
units: degrees_north
_FillValue: 32767
unlimited dimensions:
current shape = (2816, 2805)
filling on,
'Latitude_VIS': <class 'netCDF4._netCDF4.Variable'>
int32 Latitude_VIS(GeoY2, GeoX2)
long_name: latitude
add_offset: 0.0
scale_factor: 0.001
units: degrees_north
_FillValue: 327670
unlimited dimensions:
current shape = (11264, 11220)
filling on,
'Latitude_WV': <class 'netCDF4._netCDF4.Variable'>
int16 Latitude_WV(GeoY1, GeoX1)
long_name: latitude
add_offset: 0.0
scale_factor: 0.01
units: degrees_north
_FillValue: 32767
unlimited dimensions:
current shape = (1408, 1402)
filling on,
'Longitude': <class 'netCDF4._netCDF4.Variable'>
int16 Longitude(GeoY, GeoX)
long_name: longitude
units: degrees_east
add_offset: 0.0
scale_factor: 0.01
_FillValue: 32767
unlimited dimensions:
current shape = (2816, 2805)
filling on,
'Longitude_VIS': <class 'netCDF4._netCDF4.Variable'>
int32 Longitude_VIS(GeoY2, GeoX2)
add_offset: 0.0
scale_factor: 0.001
long_name: longitude
units: degrees_east
_FillValue: 327670
unlimited dimensions:
current shape = (11264, 11220)
filling on,
'Longitude_WV': <class 'netCDF4._netCDF4.Variable'>
int16 Longitude_WV(GeoY1, GeoX1)
long_name: longitude
add_offset: 0.0
scale_factor: 0.01
units: degrees_east
_FillValue: 32767
unlimited dimensions:
current shape = (1408, 1402)
filling on,
'SCAN_LINE_TIME': <class 'netCDF4._netCDF4.Variable'>
vlen SCAN_LINE_TIME(GeoY1)
long_name: Scan Time for Water Vapor Resolution
vlen data type: <class 'str'>
unlimited dimensions:
current shape = (1408,),
'Sat_Azimuth': <class 'netCDF4._netCDF4.Variable'>
uint16 Sat_Azimuth(time, GeoY, GeoX)
long_name: Satellite Azimuth
add_offset: 0.0
scale_factor: 0.01
units: degree
_FillValue: 65535
coordinates: time Latitude Longitude
unlimited dimensions:
current shape = (1, 2816, 2805)
filling on,
'Sat_Elevation': <class 'netCDF4._netCDF4.Variable'>
int16 Sat_Elevation(time, GeoY, GeoX)
long_name: Satellite Elevation
add_offset: 0.0
units: degree
scale_factor: 0.01
_FillValue: 32767
coordinates: time Latitude Longitude
unlimited dimensions:
current shape = (1, 2816, 2805)
filling on,
'Sun_Azimuth': <class 'netCDF4._netCDF4.Variable'>
uint16 Sun_Azimuth(time, GeoY, GeoX)
long_name: Sun Azimuth
add_offset: 0.0
units: degree
scale_factor: 0.01
_FillValue: 65535
coordinates: time Latitude Longitude
unlimited dimensions:
current shape = (1, 2816, 2805)
filling on,
'Sun_Elevation': <class 'netCDF4._netCDF4.Variable'>
int16 Sun_Elevation(time, GeoY, GeoX)
long_name: Sun Elevation
add_offset: 0.0
scale_factor: 0.01
units: degree
_FillValue: 32767
coordinates: time Latitude Longitude
unlimited dimensions:
current shape = (1, 2816, 2805)
filling on,
'time': <class 'netCDF4._netCDF4.Variable'>
float64 time(time)
units: minutes since 2000-01-01 00:00:00
unlimited dimensions:
current shape = (1,)
filling off}
data.variables.keys()
Out: dict_keys(['GeoX', 'GeoX1', 'GeoX2', 'GeoY', 'GeoY1', 'GeoY2', 'GreyCount', 'IMG_MIR', 'IMG_MIR_RADIANCE', 'IMG_MIR_TEMP', 'IMG_SWIR', 'IMG_SWIR_RADIANCE', 'IMG_TIR1', 'IMG_TIR1_RADIANCE', 'IMG_TIR1_TEMP', 'IMG_TIR2', 'IMG_TIR2_RADIANCE', 'IMG_TIR2_TEMP', 'IMG_VIS', 'IMG_VIS_ALBEDO', 'IMG_VIS_RADIANCE', 'IMG_WV', 'IMG_WV_RADIANCE', 'IMG_WV_TEMP', 'Latitude', 'Latitude_VIS', 'Latitude_WV', 'Longitude', 'Longitude_VIS', 'Longitude_WV', 'SCAN_LINE_TIME', 'Sat_Azimuth', 'Sat_Elevation', 'Sun_Azimuth', 'Sun_Elevation', 'time'])
lats = data.variables['Latitude'][:]
lons = data.variables['Longitude'][:]
tir = data.variables['IMG_TIR1'][:]
lats.shape
Out:(2816, 2805)
lons.shape
Out:(2816, 2805)
tir.shape
Out:(1, 2816, 2805)
` I was not so sure about exact lat and lon that I want to plot so here's two of them. Although I think the 1st one is the correct one`
#mp = Basemap(projection = 'merc',
# llcrnrlon = 68.766521,
# llcrnrlat = 3.323179,
# urcrnrlon = 100.202163,
# urcrnrlat = 37.315495,
# resolution = 'i')
mp = Basemap(projection = 'merc',
llcrnrlon = 0.8432964,
llcrnrlat = -81.04153,
urcrnrlon = 163.15671,
urcrnrlat = 81.04153,
resolution = 'i')
` I'm getting an error in this line dk what it is`
lon, lat = np.meshgrid(lons, lats)
Out: ---------------------------------------------------------------------------
MemoryError Traceback (most recent call last)
~\AppData\Local\Temp/ipykernel_12060/1726626783.py in <module>
----> 1 lon, lat = np.meshgrid(lons, lats)
~\anaconda3\lib\site-packages\numpy\core\overrides.py in meshgrid(*args, **kwargs)
~\anaconda3\lib\site-packages\numpy\lib\function_base.py in meshgrid(copy, sparse, indexing, *xi)
4947
4948 if copy:
-> 4949 output = [x.copy() for x in output]
4950
4951 return output
~\anaconda3\lib\site-packages\numpy\lib\function_base.py in <listcomp>(.0)
4947
4948 if copy:
-> 4949 output = [x.copy() for x in output]
4950
4951 return output
~\anaconda3\lib\site-packages\numpy\ma\core.py in wrapped_method(self, *args, **params)
2576 """
2577 def wrapped_method(self, *args, **params):
-> 2578 result = getattr(self._data, funcname)(*args, **params)
2579 result = result.view(type(self))
2580 result._update_from(self)
MemoryError: Unable to allocate 227. TiB for an array with shape (7898880, 7898880) and data type float32
x,y = mp(lon, lat)
Out:---------------------------------------------------------------------------
NameError Traceback (most recent call last)
~\AppData\Local\Temp/ipykernel_6908/3885003198.py in <module>
----> 1 x,y = mp(lon, lat)
NameError: name 'lon' is not defined
c_scheme = mp.pcolor(x, y, np.squeeze(tir1[0, :, :]), cmap = 'jet')
Out: ---------------------------------------------------------------------------
NameError Traceback (most recent call last)
~\AppData\Local\Temp/ipykernel_6908/254152092.py in <module>
----> 1 c_scheme = mp.pcolor(x, y, np.squeeze(tir1[0, :, :]), cmap = 'jet')
NameError: name 'x' is not defined
更新 2022-03-08：Pseudo-code 在我原来的 post 中使用 plt.contourf() 绘制数据。这不适用于保存在 'IMG_TIR1' 数据集中的栅格图像数据。相反，必须使用 plt.imshow()。有关完整过程，请参阅 2022 年 3 月 7 日编辑 post 的详细答案。
所有'IMG_name'数据集都是栅格图像数据。我修改了代码以显示一个简单示例，该示例使用 h5py 读取 HDF5 文件，然后使用 matplotlib 和 cartopy 在地球静止（原始卫星）投影上绘制图像数据。我根据 Basemap 网站上的评论选择了 cartopy：“弃用通知：不推荐使用 Basemap 以支持 Cartopy 项目。”我确定可以修改它以使用 Basemap 绘制和 NetCDF 读取 HDF5 文件。
更新代码如下：
import h5py
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
# First get data from HDF5 file with h5py:
fn = '3DIMG_30MAR2018_0000_L1B_STD.h5' #filename (the ".h5" file)
with h5py.File(fn) as data:
tir1 = data['IMG_TIR1'][:].reshape(lats.shape)
# retrieve extent of plot from file attributes:
left_lon = data.attrs['left_longitude'][0]
right_lon = data.attrs['right_longitude'][0]
lower_lat = data.attrs['lower_latitude'][0]
upper_lat = data.attrs['upper_latitude'][0]
sat_long = data.attrs['Nominal_Central_Point_Coordinates(degrees)_Latitude_Longitude'][1]
# Define extent of plot (in degrees) as:
# extent=[longitude_top_left,longitude_top_right,
# latitude_bottom_left,latitude_top_left]
# coordinates must be transformed to map coordinate system before using
img_extent_deg = (left_lon, right_lon, lower_lat, upper_lat)
# Create Geostationary plot with cartopy and matplotlib
map_proj = ccrs.Geostationary(central_longitude=sat_long)
ax = plt.axes(projection=map_proj)
ax.coastlines(color='white')
map_extend_geos = ax.get_extent(crs=map_proj)
plt.imshow(tir1, extent=map_extend_geos)
plt.colorbar(im1)
虽然我找不到数据，但我使用了 The HDF Group 的 HDFView 从数据中生成图像。这是 IMG_TIR1。
此时，我需要有关文件架构的更多信息来“查找”此数据。 MOSDAC 网站是否有任何描述文件架构的内容，或者是否有一些显示如何读取数据的代码？
@Black Viking，回答你的问题比我最初想象的要复杂。我从评论中的代码开始。经过一番挖掘，我确定每个 'IMG_xxx' 数据集中的数据是原始栅格图像（扫描），相关经度和纬度数据集中的值是每个像素的 (lon,lat) 位置。它们不形成规则的 MxN 网格，因此不能与 matplotlib 的 contourf() 函数一起使用。所以，我之前的回答是不正确的。 （我会删除其中的大部分内容，以免让未来的读者感到困惑。）
当我从事这项工作时，我发现在具有不同坐标系的地图上绘制卫星光栅图像并不是一项简单的任务（从地球静止投影到墨卡托投影）。您必须做的事情：
- 定义图像的原点坐标参考系。
- 在地球静止系统中绘制图像时，要显示轴，您需要将
extent=添加到plt.imshow()。 - 在不同坐标系中绘制图像时，必须使用
'transform='和plt.imshow()将图像转换到新坐标系。
我将逐步 post 代码，以展示如何从对地静止系统中的简单栅格图转到聚焦于印度地区的墨卡托投影中的同一图像。请阅读到最后 - 如果某些内容不起作用，您可能需要一些其他信息。
情节 0：起点来自您的评论（略有修改）。这仅绘制原始 IMG_TIR1 图像数据。
import h5py
import matplotlib.pyplot as plt
fn = '3DIMG_30MAR2018_0000_L1B_STD.h5' #filename (the ".h5" file)
with h5py.File(fn) as f:
img_arr = f['IMG_TIR1'][0,:,:]
fig = plt.subplots(figsize=(10,10))
plt.title('plot raw IMG_TIR1 image data (only)')
im = plt.imshow(img_arr)
plt.colorbar(im)
Plot 0（带海岸线，图像经过修剪）：这是将海岸线添加到原始 IMG_TIR1 图像数据的小修改。这是我验证墨卡托投影在以下图像中正确定向的参考。请注意它如何使用 extent=img_extent_sat 以便正确显示轴。这很棘手，因为地球静止坐标是 而不是 （经度、纬度）。 ax.get_extent(projection=) returns 该投影的地图范围。我不得不稍微增加它们以使图像“适合”轴。
import h5py
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
fn = '3DIMG_30MAR2018_0000_L1B_STD.h5' #filename (the ".h5" file)
with h5py.File(fn) as f:
img_arr = f['IMG_TIR1'][0,:,:]
map_proj = ccrs.Geostationary(central_longitude=82.0)
fig = plt.figure(figsize=(10,10))
ax = plt.axes(projection=map_proj)
ax.coastlines(color='white')
# Image extent in Geostationary coordinates:
img_extent_sat = ax.get_extent(crs=map_proj)
img_extent_sat = [1.04*x for x in img_extent_sat]
im = plt.imshow(img_arr, extent=img_extent_sat)
plt.colorbar(im)
plt.title('plot raw IMG_TIR1 data + Geostationary coastlines')
这是在墨卡托投影中创建绘图的代码。为避免重复，我 post 为每个步骤编辑了代码段。使用第一个（“基础”）部分并为您想要的情节添加部分。
基本代码：从 HDF5 文件中获取图像数据和属性以及卫星属性。我定义了 image = 'IMG_TIR1' 以便您可以轻松修改以获取其他 IMG_name 数据。仅当您要从图像数据（计数值）计算辐射率时才需要 Radiance 属性。最后一行创建了一个屏蔽数组，屏蔽了角落里的_FillValue数据并简化了图像绘制。
import h5py
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
fn = '3DIMG_30MAR2018_0000_L1B_STD.h5' #filename (the ".h5" file)
with h5py.File(fn) as f:
# retrieve image data:
image = 'IMG_TIR1'
img_arr = f[image][0,:,:]
# get _FillValue for data masking
img_arr_fill = f[image].attrs['_FillValue'][0]
# retrieve extent of plot from file attributes:
left_lon = f.attrs['left_longitude'][0]
right_lon = f.attrs['right_longitude'][0]
lower_lat = f.attrs['lower_latitude'][0]
upper_lat = f.attrs['upper_latitude'][0]
sat_long = f.attrs['Nominal_Central_Point_Coordinates(degrees)_Latitude_Longitude'][1]
sat_hght = f.attrs['Observed_Altitude(km)'][0] * 1000.0 # (for meters)
# retrieve attributes to calculate radiance from count:
Sensor_Name = f.attrs['Sensor_Name'].decode('utf-8')
img_fill = f[image].attrs['_FillValue'][0]
img_inv = f[image].attrs['invert'].decode('utf-8')
img_lrquad = f[image].attrs['lab_radiance_quad'][0]
img_lrscale = f[image].attrs['lab_radiance_scale_factor'][0]
img_lroff = f[image].attrs['lab_radiance_add_offset'][0]
print('Done reading HDF5 file')
## Use np.ma.masked_equal with integer values to
## mask '_FillValue' data in corners:
img_arr_m = np.ma.masked_equal(img_arr, img_arr_fill, copy=True)
绘图 1：将蒙版图像数据投影到全墨卡托投影。请注意 ax1.imshow() 如何使用 transform=data_crs，其中 data_crs 引用图像数据的地球静止投影。注意：如果您绘制 img_arr 而不是 img_arr_m，您将在图像周围看到黄色的“_FillValue”“光环”。此外，我在 ax# 和 im# 中添加了数字，以简化将以下 3 个图形绘制为 1 个图形的子图。
# Plot 1 uses Mercator projection
print('Start on Plot1')
map_proj = ccrs.Mercator(central_longitude=sat_long,
min_latitude=lower_lat, max_latitude=upper_lat)
data_crs = ccrs.Geostationary(central_longitude=sat_long,
satellite_height=sat_hght)
plt.figure(figsize=(10,10))
ax1 = plt.axes(projection=map_proj)
ax1.coastlines()
#ax1.add_feature(cfeature.BORDERS, edgecolor='white', linewidth=0.5)
ax1.gridlines(color='black', alpha=0.5, linestyle='--', linewidth=0.75, draw_labels=True)
map_proj_text = f'{str(type(map_proj)).split(".")[-1][:-2]}'
data_crs_text = f'{str(type(data_crs)).split(".")[-1][:-2]}'
plt.title(f'Plot1: Projection: {map_proj_text}\n' + \
f'Data Transform: {data_crs_text}\n' + \
f'\nRaster Data: {image} (masked)')
print('plotting data for Plot1 image')
im1 = ax1.imshow(img_arr_m, origin='upper', transform=data_crs)
plt.colorbar(im1)
print('done with Plot1')
绘图 2：将蒙版图像数据投影到经过修剪的墨卡托投影。地图大小由 ax2.set_extent() 中的 map_extent_merc 参数定义。请注意，墨卡托坐标不是度数，因此我创建了一个函数来转换范围值 (transform_extent_pts())。 （对于您的数据，它从 PlateCarree 中的 map_extent_deg 变为墨卡托中的 map_extent_merc）。尝试 map_extent_deg 的不同值以查看地图如何变化。
# Plot 2 on Mercator projection
print('Start on Plot2')
map_proj = ccrs.Mercator(central_longitude=sat_long,
min_latitude=lower_lat, max_latitude=upper_lat)
data_crs = ccrs.Geostationary(central_longitude=sat_long,
satellite_height=sat_hght)
plt.figure(figsize=(10,10))
ax2 = plt.axes(projection=map_proj)
ax2.coastlines()
ax2.add_feature(cfeature.BORDERS, edgecolor='white', linewidth=0.5)
ax2.gridlines(color='black', alpha=0.5, linestyle='--', linewidth=0.75, draw_labels=True)
# Focus map on Indian subcontinent:
deg_crs = ccrs.PlateCarree()
map_extent_deg = (15., 150., -67., 67.) # lon/lat focused on image
#map_extent_deg = (60.0, 97.0, -10.0, 37.5) # India
map_extent_merc = transform_extent_pts(map_extent_deg, map_proj, deg_crs)
ax2.set_extent(map_extent_merc, map_proj)
map_proj_text = f'{str(type(map_proj)).split(".")[-1][:-2]}'
data_crs_text = f'{str(type(data_crs)).split(".")[-1][:-2]}'
plt.title(f'Plot2: Projection: {map_proj_text}\n' + \
f'Data Transform: {data_crs_text}\n' + \
f"\nRaster Data: {image} (masked) ")
print('plotting data for Plot2 image')
im2 = plt.imshow(img_arr_m, origin='upper', transform=data_crs)
plt.colorbar(im2)
print('done with Plot2')
图 2 中使用的转换函数:
def transform_extent_pts(extent_pts, map_proj, pt_crs):
xul, yul = map_proj.transform_point(
x = extent_pts[0],
y = extent_pts[3],
src_crs = pt_crs)
xlr, ylr = map_proj.transform_point(
x = extent_pts[1],
y = extent_pts[2],
src_crs = pt_crs)
return [xul, xlr, ylr, yul]
绘图 3：计算并显示辐射值。以前的图显示原始图像数据（即“计数”）。您提供的论文具有将“计数”转换为辐射率的方程式。我写了一个函数来做到这一点 (calc_radiance())。这只需要从上面进行一些小的更改，如下所述。
首先，添加辐射功能。
def transform_extent_pts(extent_pts, map_proj, pt_crs):
xul, yul = map_proj.transform_point(
x = extent_pts[0],
y = extent_pts[3],
src_crs = pt_crs)
xlr, ylr = map_proj.transform_point(
x = extent_pts[1],
y = extent_pts[2],
src_crs = pt_crs)
return [xul, xlr, ylr, yul]
在之前添加以下行：print('Start on Plot2') 以计算辐射值。此数据用于图 3：
# Plot 3: Plot radiance on Mercator projection
print('Calculate radiance for Plot3')
radiance = calc_radiance(img_arr_m, img_lrquad, img_lrscale, img_lroff,
invert=img_inv, Sensor_Name=Sensor_Name)
接下来，将轴引用从 ax2 修改为 ax3。
最后修改最后3行：
plt.title(f'Plot3: Projection: {map_proj_text}\n' + \
f'Data Transform: {map_proj_text}\n' + \
f'\nRadiance from: {image} (masked)')
print('plot image for Plot3')
im3 = ax3.imshow(radiance, origin='upper', transform=data_crs)
plt.colorbar(im3)
print('done with Plot3')
图 4：调整图像数据以绘制云阈值。这类似于图 3。您可以使用 np.where() 函数创建一个使用阈值（小于、大于或两者）的新图像数组。此外，我添加了 vmin=, vmax= 参数以将颜色条范围与您的绘图相匹配。只需要做一些小的改动。为简洁起见，下面仅说明更改。
# Plot 4: Plot cloud thresholds on Mercator projection
print('Calculate cloud mask for Plot4')
# Left Plot: set values <=700 to 0
cloud1 = np.where(img_arr_m <= 700, 0, img_arr_m)
# Right Plot: set values <=700 to 0, >700 to 1
cloud2 = np.where(img_arr_m <= 700, 0, 1)
.....
.....
print('plot image for Plot4')
# to create plot on left:
im4 = ax4.imshow(cloud1, vmin=0, vmax=1000, origin='upper', transform=data_crs)
# to create plot on right:
im4 = ax4.imshow(cloud2, vmin=0, vmax=1, origin='upper', transform=data_crs)
plt.colorbar(im4)
背景和参考信息：
注意：我在使用早期版本转换 GeoStationary 图像时遇到了问题。我不得不使用带有 Proj 8.2.1 的 cartopy v0.20.2 来获取这些图像。如果有问题，运行这个例子来测试：Reprojecting images from a Geostationary projection
另外，这里有 2 个有用的博客 post 的链接（均包含 Python 代码）
- Blog tutorial 附有在 cartopy 中绘制整个地球图像的详细指南
- Another blog tutorial 展示了如何使用 cartopy 从 WMS 绘制卫星图像。
最后，Whosebug 对此主题有一些很好的问答。使用 tag: [cartopy] 'satellite' 搜索以找到一些好文章。
对于以后可能阅读这些内容的人：@kcw78 提供的代码是正确的，您只需要对图像数据的范围做一些小的改动。对于我来说，我将印度地区的 map_extent_deg = (15., 150., -67., 67.) 更改为 map_extent_deg = (65., 100., 0., 40.)。像这样你可能需要根据你的需要改变它。
对于以后可能阅读这些内容的人：@kcw78 提供的代码是正确的，您只需要对图像数据的范围做一些小的改动。对我来说，我将它从 map_extent_deg = (15., 150., -67., 67.) 更改为 map_extent_deg = (65., 100., 0., 40.)。
获得相同结果的另一个代码：
import cv2
import numpy as np
import matplotlib.pyplot as plt
fn = '3DIMG_30MAR2018_0000_L1B_STD.h5' #filename (the ".h5" file)
hf = h5py.File(fn, 'r')
hf.keys()
data = hf.get('IMG_TIR1')
data = data[0,:,:]
data=np.where(data>1022, np.nan, data)
Lat= hf.get('Latitude')
Lat = Lat[:,:]
Lat=Lat*0.0099999998
Lat=np.where(Lat>312, np.nan, Lat)
Long= hf.get('Longitude')
Long = Long[:,:]
Long=Long*0.0099999998
Long=np.where(Long>312, np.nan, Long)
X=Long[427:1421,804:1902]
Y=Lat[427:1421,804:1902]
Z=data[427:1421,804:1902]
fig = plt.figure(figsize=(10,10))
plt.pcolor(X, Y, Z)
plt.xlim(60, 100)
plt.ylim(0, 40)
fig = plt.figure(figsize=(10,10))
plt.pcolor(X, Y, Z)
plt.colorbar(label="Radiance", orientation="vertical")
plt.xlim(60, 100)
plt.ylim(0, 40)
#Threshold
Zn=np.where(Z>810, 1, 0)
fig = plt.figure(figsize=(10,10))
plt.pcolor(X, Y, Zn, cmap='jet', vmin=0, vmax=1)
plt.colorbar(label="Radiance", orientation="vertical")
plt.xlim(60, 100)
plt.ylim(0, 40)