Issue
I got a HDF5 file from MOSDAC website named 3DIMG_30MAR2018_0000_L1B_STD.h5 and I'm tying to read the this file and visualize the TIR and VIS dataset in the file in map. Since I am new to programming I don't have much idea how to do this for a specific location on the global data using the Latitude and longitude values. My problem statement is to show data for Indian region.
The coding I have done till now:
#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\offic\3DIMG_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
Solution
@Black Viking, answering your question was more complicated than I initially perceived. I started from your code in the comments. After some digging, I determined the data in each 'IMG_xxx' dataset is a raw raster image (scan), and the values in the associated longitude and latitude datasets are the (lon,lat) locations for each pixel. They don't form a regular MxN grid, so can't be used with matplotlib's contourf() function. So, my previous answer is incorrect. (I will delete most of it to avoid confusing future readers.)
As I worked on this, I discovered plotting satellite raster images on a map with a different coordinate system is not a simple task (going from Geostationary to a Mercator projection). Things you have to do:
- Define the origin coordinate reference system for the image.
- When plotting an image in a Geostationary system, to display the axes you need to add
extent=toplt.imshow(). - When plotting an image in a different coordinate system, you must transform the image to the new system with
'transform='andplt.imshow().
I am going to post the code incrementally to show how to go from a simple raster plot in a Geostationary system to the same image in a Mercator projection focused on the India region. Please read to the end -- there is some additional info you might need if something doesn't work.
Plot 0: Starting point from your comments (modified slightly). This only plots the raw IMG_TIR1 image data only.
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 (w/ coastlines, image trimmed): This is a small modification to add the coastlines to the raw IMG_TIR1 image data. This was my reference to verify Mercator projections were correctly oriented in the following images. Note how it uses extent=img_extent_sat so the axes are displayed correctly. This is tricky because Geostationary Coordinates are not (longitude, latitude). ax.get_extent(projection=) returns the map extents for that projection. I had to increase them slightly to get the image to "fit" the axes.
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')
Here is the code to create plots in the Mercator projection. To avoid repetition, I posted code segments for each step. Use the first ('base") section and add sections for the plot you want.
Base Code: Gets image data and attributes along with satellite attributes from HDF5 file. I defined image = 'IMG_TIR1' so you can easily modify to get other IMG_name data. The Radiance attributes are only required if you want to calculate radiance from the image data (count values). The last line creates a masked array that maskes the_FillValue data in the corners and simplifies image plotting.
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)
Plot 1: Projects masked image data to full Mercator projection. Notice how ax1.imshow() uses transform=data_crs where data_crs references a Geostationary projection of the image data. Note: If you plot img_arr instead of img_arr_m, you will see the yellow '_FillValue' "halo" around the image. Also, I added numbers to ax# and im# to simplify plotting the following 3 figures as subplots on 1 figure.
# 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')
Plot 2: Projects masked image data to a trimmed Mercator projection. The map size is defined with the map_extent_merc parameter in ax2.set_extent(). Note that Mercator coordinates are not in degrees, so I created a function to transform the extent values (transform_extent_pts()). (For your data, it goes from map_extent_deg in PlateCarree to map_extent_merc in Mercator). Try different values of map_extent_deg to see how the map changes.
# 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')
Transform function used in Plot 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]
Plot 3: Calculates and displays radiance values. Previous plots display the raw image data (which is "count"). The paper you provided has the equation to convert "count" to radiance. I wrote a function to do that (calc_radiance()). This only requires a few minor changes from above as explained below.
First, add the radiance function.
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]
Add the following lines BEFORE: print('Start on Plot2') to calculate the radiance values. This data is used for Plot 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)
Next, modify axis references from ax2 to ax3.
Finally, modify the following 3 lines at the end:
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')
Plot 4: Adjust image data to plot cloud thresholds. This is similar to Plot 3. You can use the np.where() function to create a new image array using thresholds (less than, greater than, or both). Also, I added vmin=, vmax= parameters to match the colorbar range to your plot. Only a few minor changes are required. For brevity, only changes are explained below.
# 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)
Background and Reference info:
Note: I had problems transforming GeoStationary images with earlier versions. I had to use cartopy v0.20.2 with Proj 8.2.1 to get these images. If you have problems, run this example to test: Reprojecting images from a Geostationary projection
Also, here are links to 2 blog posts that are helpful (both include Python code):
- Blog tutorial with a detailed guide to plotting whole earth imagery in cartopy
- Another blog tutorial that shows how to plot a satellite image from a WMS with cartopy.
Finally, StackOverflow has some good Q&A on this topic. Search using tag: [cartopy] 'satellite' to find some good articles.
Answered By - kcw78
0 comments:
Post a Comment
Note: Only a member of this blog may post a comment.