[FIXED] How to plot the map correctly over the SST data in cartopy?

Issue

I am trying to plot L2 Sea Surface Temperature data and I want to plot it over the globe in a geostationary projection. I tried the following code:

import h5py
import sys
import numpy as np
import matplotlib.pyplot as plt    
import cartopy.crs as ccrs
import cartopy.feature as cfeature

# First get data from HDF5 file with h5py:
fn = '/home/swadhin/project/insat/data/3RIMG_30MAR2018_0014_L2B_SST_V01R00.h5'

with h5py.File(fn) as f: 
    print(list(f.keys()))
    image = 'SST'
    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['Nominal_Altitude(km)'][0] * 1000.0 # (for meters)
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)
print(img_arr_fill)
print(np.max(img_arr_m))
print(np.min(img_arr_m))
#print(np.shape(img_arr_m))


# # Create Geostationary plot with cartopy and matplotlib  


map_proj = ccrs.Geostationary(central_longitude=sat_long,satellite_height=sat_hght)
ax = plt.axes(projection=map_proj)


ax.coastlines(color='black',linewidth = 0.5)
#ax.add_feature(cfeature.BORDERS, edgecolor='white', linewidth=0.25)
#ax.add_feature(cfeature.STATES,edgecolor = 'red',linewidth = 0.5)
ax.gridlines(color='black', alpha=0.5, linestyle='--', linewidth=0.75, draw_labels=True)
#ax.add_geometries(ind_shapes,crs = map_proj, edgecolor = 'black', alpha = 0.5)

map_extend_geos = ax.get_extent(crs=map_proj)
plt.imshow(img_arr_m, interpolation='none',origin='upper',extent=map_extend_geos, cmap = 'jet')
plt.colorbar()
#plt.clim(-10,5)
plt.savefig('/home/swadhin/project/insat/data/l2_sst.png',format = 'png', dpi=1000)

The output I got is not very accurate. There are some SST values over some of the land areas which should not be the case.

I am adding the data here for people who wanna give it a try.

https://drive.google.com/file/d/126oW36JXua-zz3XMUcyZxwPj8UISDgUM/view?usp=sharing

Solution

I have checked your HDF5 file, and there are Longitude and Latitude variables in the file. So I think these WGS84 coordinates should be used.

First, the imshow method needs the image boundary information that cannot be obtained.

I also tried the pcolormesh method, but this method can not accept lon/lat array with NaN value.

In conclusion, the contourf seems to be the best choice, but this method still has the disadvantage that it is time-consuming to run.

import h5py
import sys
import numpy as np
import matplotlib.pyplot as plt    
import cartopy.crs as ccrs
import cartopy.feature as cfeature

fn ='3RIMG_30MAR2018_0014_L2B_SST_V01R00.h5'

with h5py.File(fn) as f: 
    print(list(f.keys()))
    image = 'SST'
    img_arr = f[image][0,:,:]
    
    lon = f['Longitude'][:]*0.01
    lat = f['Latitude'][:]*0.01
#     # 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['Nominal_Altitude(km)'][0] * 1000.0 # (for meters)
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)
print(img_arr_fill)
print(np.max(img_arr_m))
print(np.min(img_arr_m))

lon_m = np.ma.masked_equal(lon, 327.67)
lat_m = np.ma.masked_equal(lat, 327.67)

# # Create Geostationary plot with cartopy and matplotlib  
map_proj = ccrs.Geostationary(central_longitude=sat_long,satellite_height=sat_hght)
# or map_proj = ccrs.PlateCarree()
ax = plt.axes(projection=map_proj)

ax.set_global()

ax.coastlines(color='black',linewidth = 0.5)
ax.add_feature(cfeature.BORDERS, edgecolor='white', linewidth=0.25)
ax.add_feature(cfeature.STATES,edgecolor = 'red',linewidth = 0.5)
ax.gridlines(color='black', alpha=0.5, linestyle='--', linewidth=0.75, draw_labels=True)

cb = ax.contourf(lon_m,lat_m,img_arr_m, cmap = 'jet',transform = ccrs.PlateCarree())
plt.colorbar(cb)
plt.savefig('l2_sst1.png',format = 'png', dpi=300)

Here is the output figure.

or using a lon-lat projection.

Answered By - Li Yupeng