Search code examples
pythoncartopymap-projections

How to offset the cartopy coastline to fit the actual data?

I am trying to plot a simple satellite image from an HDF5 file using imshow() with cartopy. However my map doesn’t exactly overlap with the correct region. Following is the sample code -

import matplotlib.pyplot as plt    
import cartopy.crs as ccrs
import cartopy.feature as cfeature

import numpy as np

# First get data from HDF5 file with h5py:
fn = '/home/swadhin/project/insat/3DIMG_15JUN2022_0930_L1C_ASIA_MER.h5' #filename (the ".h5" file)
with h5py.File(fn) as data:      
    # retrieve image data:
    image = 'IMG_TIR1'
    img_arr = data[image][0,:,:] 
    # get _FillValue for data masking
    img_arr_fill = data[image].attrs['_FillValue'][0] 


# 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]
img_extent_deg = (left_lon, right_lon, lower_lat, upper_lat)
print(img_extent_deg)
img_arr_m = np.ma.masked_equal(img_arr, img_arr_fill, copy=True)   
map_proj = ccrs.Mercator(central_longitude=sat_long)
ax = plt.axes(extent = img_extent_deg,projection=map_proj)
ax.coastlines(color='white')
ax.add_feature(cfeature.BORDERS, edgecolor='white', linewidth=0.5)
ax.add_feature(cfeature.STATES,edgecolor = 'red',linewidth = 0.5)
ax.gridlines(color='black', alpha=0.5, linestyle='--', linewidth=0.75, draw_labels=True)
map_extend_mer = ax.get_extent(crs=map_proj)
plt.imshow(img_arr_m, extent=map_extend_mer, cmap = 'gray')
plt.colorbar()

How ever the plot generated ranges slightly below from 10S and hence, the features are not correctly displayed. For example the clouds which should be over Gujurat region are now on Pakistan border. enter image description here

Update 1: The problem was resolved after I added the latitudinal extent in the map projection.


map_proj = ccrs.Mercator(central_longitude=sat_long, min_latitude= lower_lat,max_latitude=upper_lat)
Solution

  • The problem was resolved after I added the latitudinal extent in the map projection.

    map_proj = ccrs.Mercator(central_longitude=sat_long, min_latitude= lower_lat,max_latitude=upper_lat)