2

I have some LSA-SAF HDF5 file data and I'm hoping to finally plot it in Python with Cartopy. I have zero experience with HDF5 files so I could be barking up the wrong tree here, but I can plot the data and the map. The major issue is that the projections don't line up. I've tried messing with the projections in both the subplot and the imshow transform argument. Since the MSG data appears to not be geolocated maybe I can't do what I was hoping easily.

My code:

FILE_NAME = 'HDF5_LSASAF_MSG_LAI_MSG-Disk_201806010000.h5' #LAI

crs = ccrs.Geostationary(central_longitude=0.0,satellite_height= 35785831)
crs2 = ccrs.PlateCarree(central_longitude=0.0) #central_longitude=0.0
fig = plt.figure(figsize=(10, 12))
ax = fig.add_subplot(1, 1, 1, projection=crs)
f = h5py.File(FILE_NAME, mode='r')
key_list = f.keys()

key_list2 = []
key_list2.append(key_list[0])

for key in key_list2:
    print(key)
    matrix = f.get(key)
    ax.add_feature(cfeature.COASTLINE.with_scale('50m'), linewidth=0.75)
    ax.add_feature(cfeature.BORDERS.with_scale('50m'), linewidth=0.5)
    ax.add_feature(cfeature.OCEAN.with_scale('50m'),alpha=0.2)
    cmap=cm.YlGn
    cmap.set_bad(alpha=0.0)

    img_extent = (-65,65,-65,65)

    ax.imshow(matrix[:], cmap=cmap, norm=colors.Normalize(vmin=-1.0,             
      vmax=7000.0), origin='upper',extent=img_extent,transform=crs2)

plt.show()

enter image description here

I had a similar problem when I was trying to plot GOES-16 data and it was resolved with satellite height calculations for lat and lon. I don't know enough about the HDF5 file hierarchy to find similar data for the MSG geostationary satellite. Any insight into whether this can be accomplished and/or HDF5 data would be very appreciated.

4
  • Please let me know if you find a solution. Struggling here too. Having difficulty with gdal too when trying to reproject the raw HDF file to GeoTiff.
    – tda
    Dec 19, 2018 at 13:54
  • So I've found a solution using gdal to get a transformed GeoTiff first. gdal_translate -a_srs "+proj=geos +h=35785831 +a=6378169 +b=6356583.8 +no_defs" -a_ullr -5568000 5568000 5568000 -5568000 HDF5:"HDF5_LSASAF_MSG_LST_MSG-Disk_201812171515.h5"://LST lst.tif then gdalwarp -t_srs EPSG:4326 -wo SOURCE_EXTRA=100 lst.tif lst2.tif. I'm looking at the MSG LST product but should be the same.
    – tda
    Dec 19, 2018 at 15:39
  • I was fortunate enough to speak with the folks at LandSAF and this is essentially the advice they gave me too. I have a bit different arguments but going through gdal was the way to go! I will post my updated code later. Dec 20, 2018 at 19:23
  • I think my parameters were based on a resolution of 3km whereas I think it's 3.1 in reality so probably a rounding error on my part - thanks for sharing though!
    – tda
    Dec 21, 2018 at 9:43

1 Answer 1

1

Just as tda mentioned I was successful with gdal as well. Here I was using the FAPAR product.

in_pathfiles = '/path/to/HDF5 files/*FAPAR*.h5' # Where .hdf5 files exist
out_pathfiles = '/path/to/new geotiff files/' # Where the new .tif file will be placed
myfiles = glob.glob(in_pathfiles) #list of all files

for f in myfiles:
    print(f),"\n"
    filename = f.split("\\")[-1]
    print "filename",out_pathfiles+filename,"\n"

    f_out = filename[:-3] + ".tif"  # splitting the .hd5 off the fileneame and making a new .tif filename
    print "f_out",out_pathfiles+f_out,"\n"

    f_rep = out_pathfiles+filename[:-3] + "_rep.tif" # create a new final .tif filename for reprojection
    print "f_rep",f_rep,"\n"

# Translating the satellite height and ellipitical values to xy values and filling the new _rep.tif file
# from the original .h5 file
os.system('gdal_translate -of GTiff -a_srs "+proj=geos +h=35785831 +a=6378169 +b=6356583.8 +no_defs"\
-a_ullr  -5568748.27576  5568748.27576 5568748.27576 -5568748.27576 "HDF5:'+ filename + '://FAPAR '+ f_out)

# Mapping the new values and filling the new _rep.tif file
os.system('gdalwarp -ot Float32 -s_srs "+proj=geos +h=35785831 +a=6378169 +b=6356583.8 +no_defs"\
-t_srs EPSG:4326 -r near -of GTiff ' + f_out + ' ' + f_rep)

Plot:

# enable gdal exceptions (instead of the silent failure which is gdal default)
gdal.UseExceptions()

fname = "/path/to/rep.tif file/"
ds = gdal.Open(fname)

print( "[ RASTER BAND COUNT ]: ", ds.RasterCount)
cols = ds.RasterXSize
print('cols = ',cols)
rows = ds.RasterYSize
print(' rows = ', rows)
bands = ds.RasterCount
print('bands = ', bands)
driver = ds.GetDriver().LongName
print('driver =', driver)

print('MetaData = ',ds.GetMetadata())

Meta = ds.GetMetadata()
#print Meta.values()
Product = Meta.values()[3]
print Product

# print various metadata for the image
geotransform = ds.GetGeoTransform()
if not geotransform is None:
    print ('Origin = (',geotransform[0], ',',geotransform[3],')')
    print ('Pixel Size = (',geotransform[1], ',',geotransform[5],')')
proj = ds.GetProjection()
print proj
inproj = osr.SpatialReference()
inproj.ImportFromWkt(proj)

print('inproj = \n', inproj)
data = ds.ReadAsArray()
crs = ccrs.Geostationary(central_longitude=0.0)
crs2 = ccrs.PlateCarree(central_longitude=0.0)
fig = plt.figure(figsize=(10, 12))
ax = fig.add_subplot(1, 1, 1, projection=crs)

ax.add_feature(cfeature.COASTLINE.with_scale('50m'), linewidth=0.75)
ax.add_feature(cfeature.BORDERS.with_scale('50m'), linewidth=0.5)
ax.add_feature(cfeature.OCEAN.with_scale('50m'),alpha=0.2)

cmap=cm.YlGn
cmap.set_bad(alpha=0.0)
#ax.set_extent([-60,60,-60,60])
img_extent = (-81.26765645410755,81.26765645410755,-74.11423113858775,74.11423113858775)
ax.imshow(data, cmap=cmap, norm=colors.Normalize(vmin=-1.0, vmax=7000.0), origin='upper'
         ,extent=img_extent,transform=crs2) 

plt.show()

enter image description here

Plot a new region and masked array for where data is 0. This allowed me to display the oceans and other areas where the data wasn't relevant:

fig = plt.figure(figsize=(10, 12))

 # enable gdal exceptions (instead of the silent failure which is gdal default)
gdal.UseExceptions()

fname = "/path/to/rep.tif file/"
ds = gdal.Open(fname)

Meta = ds.GetMetadata()

Product = Meta.values()[3]
#print Product

Date = Meta.values()[38]
Date_End = Date[:8]

geotransform = ds.GetGeoTransform()
data = ds.ReadAsArray()
data = np.ma.masked_where(data <= -1, data)

crs = ccrs.Geostationary(central_longitude=0.0)
crs2 = ccrs.PlateCarree(central_longitude=0.0)

ax = fig.add_subplot(1, 1, 1, projection=crs2)
gl = ax.gridlines(crs=crs2, draw_labels=True,
    linewidth=2, color='gray', alpha=0.5, linestyle='--')
gl.xlabels_top = False
gl.ylabels_left = False

ax.add_feature(cfeature.COASTLINE.with_scale('50m'), linewidth=0.75)
ax.add_feature(cfeature.BORDERS.with_scale('50m'), linewidth=0.5)

cmap=cm.YlGn
cmap.set_bad(alpha=0.0)

ax.set_extent([5,40,-10,8]) # Congo

img_extent = (-81.26765645410755,81.26765645410755,-74.11423113858775,74.11423113858775)

cf = ax.imshow(data, cmap="RdYlGn", origin='upper'
    ,extent=img_extent,transform=crs2)

cbar = plt.colorbar(cf, orientation='horizontal')
ax.add_feature(cfeature.OCEAN.with_scale('50m'),alpha=0.5)

plt.show()

enter image description here

2
  • Are you using jupyter notebok? gdal_translate is not creating a Tiff file
    – eLg
    Jun 2, 2021 at 4:44
  • What's the role of 'inproj' in the code?
    – vicemagui
    Feb 1, 2022 at 16:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.