Plotting ILWIS raster and vector data in conjunction with Cartopy

Notebook prepared by Ben Maathuis and Bas Retsios. ITC-University of Twente, Enschede. The Netherlands

At the end of your raster analysis, you may want to present the results using proper cartographic techniques. This notebook demonstrates how to create a basemap for a specified region of interest, which can then serve as the background for visualizing the results of your analysis.

In this example, the Nile Basin is used as the area of interest. NDVI raster data (sourced from NOAA STAR - https://www.star.nesdis.noaa.gov/smcd/emb/vci/VH/index.php) processed using ILWIS, but saved in a geotif format along with vector data in shapefile format are combined using Cartopy to generate the maps. The resulting plots are saved as PNG files, allowing you to easily include them in a Word document or other reports.

We continue with a gauge corrected radar rainfall time series obtained from the National Meteorological Service of the Netherlands (KNMI), consiting of 11 hours of rainfall maps with a temporal interval of 5 minutes. The data is imported, resampled and animated.

Some additional installion requirements for the creation of Animations

Before you continue: Installation of FFmpeg (see also: https://ffmpeg.org/download.html).

FFmpeg is a universal media converter. It can read a wide variety of inputs - including live grabbing/recording devices - filter, and transcode them into various output formats.

A step-by-step installation guide:

1. Download FFmpeg.zip file (see FFmpeg download link above, but also available from https://filetransfer.itc.nl/pub/52n/ilwis_py/temporary_tools/ffmpeg.zip). Extract the downloaded .ZIP archive. Rename the extracted folder to something simple like "FFmpeg" for easier access.
2. Move the folder: Move the FFmpeg folder to a convenient location, such as the root of your C: drive (e.g., C:\FFmpeg).
3. Update the PATH environment variable: Search for "Edit the system environment variables" in the Windows search bar and open it. Go to Advanced -> Environment Variables. Under "System variables", find and select "Path", then click "Edit". Click "New" and add the path to the bin folder inside your FFmpeg folder (e.g., C:\FFmpeg\bin). Click "OK" on all open windows to save the changes.
4. Verify the installation: Open Command Prompt or PowerShell. Type ffmpeg -version and press Enter. If FFmpeg is installed correctly, you should see information about the FFmpeg version and build. By following these steps, you'll be able to use FFmpeg commands from any command prompt window on your system.
5. Eventually reboot your system

#load packages
import os
import ilwis
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import matplotlib.animation as animation
from cartopy import config
import cartopy
from cartopy import crs as ccrs, feature as cfeature
import cartopy.io.img_tiles as cimgt
from cartopy.io.shapereader import Reader
import numpy as np
import h5py
from datetime import datetime, timedelta
from IPython.display import HTML

import warnings
warnings.filterwarnings("ignore")

#sample data expected in folder 'nile_data' under this notebook folder
work_dir = os.getcwd()+r'/Carto_data'

print("current dir is: %s" % (os.getcwd()))
print("current working directory is:",work_dir) 

if os.path.isdir(work_dir):
    print("Folder exists")
else:
    print("Folder doesn't exists")
    os.makedirs(work_dir, exist_ok=True)

current dir is: d:\jupyter\notebook_scripts\ilwispy_tutorial\tutorial_release
current working directory is: d:\jupyter\notebook_scripts\ilwispy_tutorial\tutorial_release/Carto_data
Folder exists

#set the working directory for ILWISPy
ilwis.setWorkingCatalog(work_dir)
print(work_dir)

d:\jupyter\notebook_scripts\ilwispy_tutorial\tutorial_release/Carto_data

Read background base map

#Retrieve background image - in ilwis format
Input_file = 'ne_nile.mpl'
print(Input_file)

# open a raster coverage object:
raster = ilwis.RasterCoverage(Input_file)

# probe its size
print(raster.size().xsize)
print(raster.size().ysize)
print(raster.size().zsize)

ne_nile.mpl
1140
2220
3

# inspect raster georeference and coordinate system
roi = raster.envelope()
grf = raster.geoReference()
csy = grf.coordinateSystem()
print()
print('region of interest = ', roi)
print()
print('proj4-info = ', csy.toProj4())
print()
print('wks-info = ', csy.toWKT())

region of interest =  21.983333 -4.983333 40.983333 32.016667

proj4-info =  +proj=longlat +a=6378137.000000000000 +b=6356752.314245179296

wks-info =  PROJCS["ne_nile.csy",GEOCS["ne_nile.csy",DATUM[" WGS 84",[DWGS84],ELLIPSOID["WGS 84",6378137.000000000000,298.257223563000],PRIMEM["Greenwich",0, AUTHORITY["EPSG",8901"]]],PROJECTION[""],,UNIT[degrees,1.0]]

b1 = ilwis.do('selection',raster,"rasterbands(0)")
b2 = ilwis.do('selection',raster,"rasterbands(1)")
b3 = ilwis.do('selection',raster,"rasterbands(2)")

#read raster as np-array and transform - reshape from 1D to 2D
b1_im = np.fromiter(iter(b1), np.ubyte, b1.size().linearSize())
b1_im = b1_im.reshape((b1.size().ysize, b1.size().xsize ))

#read raster as np-array and transform - reshape from 1D to 2D
b2_im = np.fromiter(iter(b2), np.ubyte, b1.size().linearSize())
b2_im = b2_im.reshape((b1.size().ysize, b1.size().xsize ))

#read raster as np-array and transform - reshape from 1D to 2D
b3_im = np.fromiter(iter(b3), np.ubyte, b1.size().linearSize())
b3_im = b3_im.reshape((b1.size().ysize, b1.size().xsize ))

# create true color composite
nc = np.dstack((b1_im, b2_im, b3_im))

# display the background image
plt.figure(figsize = (4,8))
plt.imshow(nc, interpolation='none')
plt.title('Colour shaded relief background')
plt.axis('off')
plt.show();

projection = csy.toProj4()
geotransform = raster.envelope()
print(projection)
print(geotransform)

+proj=longlat +a=6378137.000000000000 +b=6356752.314245179296
21.983333 -4.983333 40.983333 32.016667

Add countries and sub basin boundary shape file

basin_shape = work_dir+r'/Sub_basin.shp'
lakes_nile = work_dir+r'/lakes.shp'
disp_bound = work_dir+'/disputed_NBI.shp'
no_disp_bound = work_dir+'/countries_seg_no_disputed.shp'

Create base map template

img_extent = (22, 41, -5, 32)
fig = plt.figure(figsize=(10, 14))
ax = plt.axes(projection=ccrs.PlateCarree())
plt.title('Nile Basin basemap')

# set a margin around the data
ax.set_xmargin(0.05)
ax.set_ymargin(0.10)

#add land mask Natural earth
#LAND_10m = cartopy.feature.NaturalEarthFeature('physical', 'land', '10m',edgecolor='face', facecolor=cartopy.feature.COLORS['land'])
#ax.add_feature(LAND_10m)

#add Natural earth shape files
ax.coastlines(resolution='50m', color='white', linewidth=1)
#ax.add_feature(cfeature.BORDERS, linewidth=0.5, color='red')
ax.add_feature(cfeature.COASTLINE, linewidth=0.5, color='red')
ax.add_feature(cfeature.RIVERS, linewidth=0.5, color='blue')
ax.add_feature(cfeature.LAKES, color='blue')

ax.add_geometries(Reader(basin_shape).geometries(),ccrs.PlateCarree(),facecolor='none', edgecolor='white', linewidth=1.50,zorder=3) 
ax.add_geometries(Reader(disp_bound).geometries(),ccrs.PlateCarree(),facecolor='none', edgecolor='red', linewidth=0.50,linestyle='dashed',zorder=2) 
ax.add_geometries(Reader(no_disp_bound).geometries(),ccrs.PlateCarree(),facecolor='none', edgecolor='red', linewidth=0.50,zorder=2) 

#background map
ax.imshow(nc, origin='upper',extent=img_extent, transform=ccrs.PlateCarree(),zorder=1)


# mark  known places to help us geo-locate ourselves - coordinates in lon, lat
#Capitals
ax.plot(31.32, 29.94, 'rs', markersize=4, zorder = 10, transform=ccrs.Geodetic()) #marker is red, e.g. bo for blue bs is square marker in blue
ax.text(31.32, 30.34,'Cairo', color= 'red', zorder = 10,transform=ccrs.Geodetic())

ax.plot(38.96, 15.26, 'rs', markersize=4, zorder = 10,transform=ccrs.Geodetic()) 
ax.text(37.96, 15.66, 'Asmara', color= 'red', zorder = 10, transform=ccrs.Geodetic())

ax.plot(32.58, 15.53, 'rs', markersize=4, zorder = 10,transform=ccrs.Geodetic()) 
ax.text(32.58, 15.93, 'Khartoum', color= 'red', zorder = 10, transform=ccrs.Geodetic())

ax.plot(38.81, 8.92, 'rs', markersize=4, zorder = 10,transform=ccrs.Geodetic()) 
ax.text(37.11, 9.32, 'Addis Ababa', color= 'red', zorder = 10, transform=ccrs.Geodetic())

ax.plot(31.62, 4.78, 'rs', markersize=4, zorder = 10,transform=ccrs.Geodetic()) 
ax.text(31.62, 5.18, 'Juba', color= 'red', zorder = 10, transform=ccrs.Geodetic())

ax.plot(32.62, 0.24, 'rs', markersize=4, zorder = 10,transform=ccrs.Geodetic()) 
ax.text(32.62, 0.64, 'Kampala', color= 'red', zorder = 10, transform=ccrs.Geodetic())

ax.plot(36.85, -1.37, 'rs', markersize=4, zorder = 10,transform=ccrs.Geodetic()) 
ax.text(36.85, -0.97, 'Nairobi', color= 'red', zorder = 10, transform=ccrs.Geodetic())

ax.plot(30.06, -1.95, 'rs', markersize=4, zorder = 10,transform=ccrs.Geodetic()) 
ax.text(29.36, -1.65, 'Kigali', color= 'red', zorder = 10, transform=ccrs.Geodetic())

ax.plot(29.38, -3.40, 'rs', markersize=4, zorder = 10,transform=ccrs.Geodetic()) 
ax.text(29.78, -3.60, 'Bujumbura', color= 'red', zorder = 10, transform=ccrs.Geodetic())

#Sub-Basin names
ax.text(28.01, 16.45, 'Main Nile', color= 'black', zorder = 10, transform=ccrs.Geodetic())
ax.text(30.5, 12.62, 'White Nile', color= 'black', zorder = 10, transform=ccrs.Geodetic())
ax.text(32.07, -1.65, 'Lake Victoria', color= 'black', zorder = 10, transform=ccrs.Geodetic())
ax.text(27.94, 0.27, 'Lake Albert', color= 'black', zorder = 10, transform=ccrs.Geodetic())
ax.text(31.85, 1.77, 'Victoria Nile', color= 'black', zorder = 10, transform=ccrs.Geodetic())
ax.text(24.83, 11.47, 'Bahr el Ghazal', color= 'black', zorder = 10, transform=ccrs.Geodetic())
ax.text(28.95, 6.28, 'Bahr el Jebel', color= 'black', zorder = 10, transform=ccrs.Geodetic())
ax.text(31.25,7.16, 'Baro Akobbo Sobat', color= 'black', zorder = 10, transform=ccrs.Geodetic())
ax.text(36,10.16, 'Blue Nile', color= 'black', zorder = 10, transform=ccrs.Geodetic())
ax.text(35.25,14.55, 'Tekeze Atbara', color= 'black', zorder = 10, transform=ccrs.Geodetic())


gl = ax.gridlines(draw_labels=True,color = 'grey', linestyle = '--', linewidth = 0.5)
gl.top_labels = False
gl.right_labels = False
gl.xlocator = mticker.FixedLocator([25, 30, 35, 40])


plt.show()
fig.savefig(work_dir+'/BaseMap_plot.png')

Create some new plots with additional raster data for the Nile Basin

Using as example NDVI data retrieved from NOAA STAR (see: https://www.star.nesdis.noaa.gov/smcd/emb/vci/VH/vh_ftp.php)

Month1 = ilwis.RasterCoverage('NB_smn202401_b.tif')
Month2 = ilwis.RasterCoverage('NB_smn202402_b.tif')
Month3 = ilwis.RasterCoverage('NB_smn202403_b.tif')

#read raster as np-array and transform - reshape from 1D to 2D
Month1_im = np.fromiter(iter(Month1), np.float64, Month1.size().linearSize())
Month1_im = Month1_im.reshape((Month1.size().ysize, Month1.size().xsize ))

#read raster as np-array and transform - reshape from 1D to 2D
Month2_im = np.fromiter(iter(Month2), np.float64, Month1.size().linearSize())
Month2_im = Month2_im.reshape((Month1.size().ysize, Month1.size().xsize ))

#read raster as np-array and transform - reshape from 1D to 2D
Month3_im = np.fromiter(iter(Month3), np.float64, Month1.size().linearSize())
Month3_im = Month3_im.reshape((Month1.size().ysize, Month1.size().xsize ))

#values outside Nile basin have been assigned -999 (as nodata)
#nodata is assigned to nan - masked smn
smn1 = np.ma.masked_where(Month1_im == -999, Month1_im)
smn2 = np.ma.masked_where(Month2_im == -999, Month2_im)
smn3 = np.ma.masked_where(Month3_im == -999, Month3_im)

# retrieve lowest and highest map values values
min_all = (np.nanmin(smn1)+np.nanmin(smn2)+np.nanmin(smn3))/3
max_all = (np.nanmax(smn1)+np.nanmax(smn2)+np.nanmax(smn3))/3
print(max_all)

0.5802000562349955

img_extent = (22, 41, -5, 32)
fig = plt.figure(figsize=(8, 12))
ax = plt.axes(projection=ccrs.PlateCarree())
plt.title('Monthly NDVI 202401')


#background map
ax.imshow(nc, origin='upper',extent=img_extent, transform=ccrs.PlateCarree(),zorder=1)

#can be changed with land mask Natural earth
#LAND_10m = cartopy.feature.NaturalEarthFeature('physical', 'land', '10m',edgecolor='face', facecolor=cartopy.feature.COLORS['land'])
#ax.add_feature(LAND_10m)

#add Natural earth shape files
ax.coastlines(resolution='50m', color='white', linewidth=1, zorder=3)
#ax.add_feature(cfeature.BORDERS, linewidth=0.5, color='red', zorder=3)
ax.add_feature(cfeature.COASTLINE, linewidth=0.5, color='red', zorder=3)
ax.add_feature(cfeature.RIVERS, linewidth=0.5, color='blue', zorder=3)
#ax.add_feature(cfeature.LAKES, color='blue',zorder=3)

#add custom shape files: ax.add_geometries(Reader(fname).geometries(),ccrs.PlateCarree(),facecolor='white', hatch='xxxx')
ax.add_geometries(Reader(basin_shape).geometries(),ccrs.PlateCarree(),facecolor='none', edgecolor='black', linewidth=0.50,zorder=3) 
ax.add_geometries(Reader(lakes_nile).geometries(),ccrs.PlateCarree(),facecolor='blue', edgecolor='blue', linewidth=0.50,zorder=3) 
ax.add_geometries(Reader(disp_bound).geometries(),ccrs.PlateCarree(),facecolor='none', edgecolor='red', linewidth=0.50,linestyle='dashed',zorder=3) 
ax.add_geometries(Reader(no_disp_bound).geometries(),ccrs.PlateCarree(),facecolor='none', edgecolor='red', linewidth=0.50,zorder=3) 

#data raster - here ndvi - see title
ax.imshow(smn1, origin='upper', extent=img_extent, transform=ccrs.PlateCarree(), vmin=min_all, vmax=max_all, cmap='YlGn', interpolation='None', zorder=2)

#Overall annotation
#Capitals
ax.plot(31.32, 29.94, 'rs', markersize=4, zorder = 10, transform=ccrs.Geodetic()) #marker is red, e.g. bo for blue bs is square marker in blue
ax.text(31.32, 30.34,'Cairo', color= 'red', zorder = 10,transform=ccrs.Geodetic())

ax.plot(38.96, 15.26, 'rs', markersize=4, zorder = 10,transform=ccrs.Geodetic()) 
ax.text(37.96, 15.66, 'Asmara', color= 'red', zorder = 10, transform=ccrs.Geodetic())

ax.plot(32.58, 15.53, 'rs', markersize=4, zorder = 10,transform=ccrs.Geodetic()) 
ax.text(32.58, 15.93, 'Khartoum', color= 'red', zorder = 10, transform=ccrs.Geodetic())

ax.plot(38.81, 8.92, 'rs', markersize=4, zorder = 10,transform=ccrs.Geodetic()) 
ax.text(37.11, 9.32, 'Addis Ababa', color= 'red', zorder = 10, transform=ccrs.Geodetic())

ax.plot(31.62, 4.78, 'rs', markersize=4, zorder = 10,transform=ccrs.Geodetic()) 
ax.text(31.62, 5.18, 'Juba', color= 'red', zorder = 10, transform=ccrs.Geodetic())

ax.plot(32.62, 0.24, 'rs', markersize=4, zorder = 10,transform=ccrs.Geodetic()) 
ax.text(32.62, 0.64, 'Kampala', color= 'red', zorder = 10, transform=ccrs.Geodetic())

ax.plot(36.85, -1.37, 'rs', markersize=4, zorder = 10,transform=ccrs.Geodetic()) 
ax.text(36.85, -0.97, 'Nairobi', color= 'red', zorder = 10, transform=ccrs.Geodetic())

ax.plot(30.06, -1.95, 'rs', markersize=4, zorder = 10,transform=ccrs.Geodetic()) 
ax.text(29.36, -1.65, 'Kigali', color= 'red', zorder = 10, transform=ccrs.Geodetic())

ax.plot(29.38, -3.40, 'rs', markersize=4, zorder = 10,transform=ccrs.Geodetic()) 
ax.text(29.78, -3.60, 'Bujumbura', color= 'red', zorder = 10, transform=ccrs.Geodetic())

#Sub-Basins
ax.text(28.01, 16.45, 'Main Nile', color= 'black', zorder = 10, transform=ccrs.Geodetic())
ax.text(30.5, 12.62, 'White Nile', color= 'black', zorder = 10, transform=ccrs.Geodetic())
ax.text(32.07, -1.65, 'Lake Victoria', color= 'black', zorder = 10, transform=ccrs.Geodetic())
ax.text(27.94, 0.27, 'Lake Albert', color= 'black', zorder = 10, transform=ccrs.Geodetic())
ax.text(31.85, 1.77, 'Victoria Nile', color= 'black', zorder = 10, transform=ccrs.Geodetic())
ax.text(24.83, 11.47, 'Bahr el Ghazal', color= 'black', zorder = 10, transform=ccrs.Geodetic())
ax.text(28.95, 6.28, 'Bahr el Jebel', color= 'black', zorder = 10, transform=ccrs.Geodetic())
ax.text(31.25,7.16, 'Baro Akobbo Sobat', color= 'black', zorder = 10, transform=ccrs.Geodetic())
ax.text(36,10.16, 'Blue Nile', color= 'black', zorder = 10, transform=ccrs.Geodetic())
ax.text(35.25,14.55, 'Tekeze Atbara', color= 'black', zorder = 10, transform=ccrs.Geodetic())

gl = ax.gridlines(draw_labels=True,color = 'grey', linestyle = '--', linewidth = 0.5)
gl.top_labels = False
gl.right_labels = False
gl.xlocator = mticker.FixedLocator([25, 30, 35, 40])

#legend
sm = plt.cm.ScalarMappable(cmap='YlGn',norm=plt.Normalize(min_all, max_all))
sm._A = []
cbar = plt.colorbar(sm,ax=ax, shrink=0.45, pad = 0.05)
cbar.set_label('NDVI index')

plt.show()
fig.savefig(work_dir+'/SMN_202401.png')

Template to plot 3 months

# Use subplot_kw to declare the projection
fig, (ax1, ax2, ax3) = plt.subplots(1,3, figsize=(15, 8), subplot_kw={"projection": ccrs.PlateCarree()})
fig.suptitle('Normalized Difference Vegetation Index (NDVI)')

#Month1
ax1.set_title('NDVI_202401')
img_extent = (22, 41, -5, 32)

#background map
ax1.imshow(nc, origin='upper',extent=img_extent, transform=ccrs.PlateCarree(),zorder=1)
# LAND_10m = cartopy.feature.NaturalEarthFeature('physical', 'land', '10m',edgecolor='face', facecolor=cartopy.feature.COLORS['land'])
# ax1.add_feature(LAND_10m)

#add Natural earth shape files
ax1.coastlines(resolution='50m', color='white', linewidth=1, zorder=3)
#ax1.add_feature(cfeature.BORDERS, linewidth=0.5, color='red', zorder=3)
ax1.add_feature(cfeature.COASTLINE, linewidth=0.5, color='red', zorder=3)
ax1.add_feature(cfeature.RIVERS, linewidth=0.5, color='blue', zorder=3)

#add custom shape files: ax.add_geometries(Reader(fname).geometries(),ccrs.PlateCarree(),facecolor='white', hatch='xxxx')
ax1.add_geometries(Reader(basin_shape).geometries(),ccrs.PlateCarree(),facecolor='none', edgecolor='black', linewidth=0.50,zorder=3) 
ax1.add_geometries(Reader(lakes_nile).geometries(),ccrs.PlateCarree(),facecolor='blue', edgecolor='blue', linewidth=0.50,zorder=3) 
ax1.add_geometries(Reader(disp_bound).geometries(),ccrs.PlateCarree(),facecolor='none', edgecolor='red', linewidth=0.50,linestyle='dashed',zorder=3) 
ax1.add_geometries(Reader(no_disp_bound).geometries(),ccrs.PlateCarree(),facecolor='none', edgecolor='red', linewidth=0.50,zorder=3) 

#data raster - see title
ax1.imshow(smn1, origin='upper', extent=img_extent, transform=ccrs.PlateCarree(), vmin=min_all, vmax=max_all, cmap='YlGn', interpolation='None', zorder=2)
gl = ax1.gridlines(draw_labels=True,color = 'grey', linestyle = '--', linewidth = 0.5)
gl.top_labels = False
gl.right_labels = False
gl.xlocator = mticker.FixedLocator([25, 30, 35, 40])


#Month2
ax2.set_title('NDVI_202402')
img_extent = (22, 41, -5, 32)

#background map
ax2.imshow(nc, origin='upper',extent=img_extent, transform=ccrs.PlateCarree(),zorder=1)
# LAND_10m = cartopy.feature.NaturalEarthFeature('physical', 'land', '10m',edgecolor='face', facecolor=cartopy.feature.COLORS['land'])
# ax2.add_feature(LAND_10m)
#add Natural earth shape files
ax2.coastlines(resolution='50m', color='white', linewidth=1, zorder=3)
#ax2.add_feature(cfeature.BORDERS, linewidth=0.5, color='red', zorder=3)
ax2.add_feature(cfeature.COASTLINE, linewidth=0.5, color='red', zorder=3)
ax2.add_feature(cfeature.RIVERS, linewidth=0.5, color='blue', zorder=3)

#add custom shape files: ax.add_geometries(Reader(fname).geometries(),ccrs.PlateCarree(),facecolor='white', hatch='xxxx')
ax2.add_geometries(Reader(basin_shape).geometries(),ccrs.PlateCarree(),facecolor='none', edgecolor='black', linewidth=0.50,zorder=3) 
ax2.add_geometries(Reader(lakes_nile).geometries(),ccrs.PlateCarree(),facecolor='blue', edgecolor='blue', linewidth=0.50,zorder=3) 
ax2.add_geometries(Reader(disp_bound).geometries(),ccrs.PlateCarree(),facecolor='none', edgecolor='red', linewidth=0.50,linestyle='dashed',zorder=3) 
ax2.add_geometries(Reader(no_disp_bound).geometries(),ccrs.PlateCarree(),facecolor='none', edgecolor='red', linewidth=0.50,zorder=3) 

#data raster - see title
ax2.imshow(smn2, origin='upper', extent=img_extent, transform=ccrs.PlateCarree(), vmin=min_all, vmax=max_all, cmap='YlGn', interpolation='None', zorder=2)
gl = ax2.gridlines(draw_labels=True,color = 'grey', linestyle = '--', linewidth = 0.5)
gl.top_labels = False
gl.right_labels = False
gl.xlocator = mticker.FixedLocator([25, 30, 35, 40])

#month3
ax3.set_title('NDVI_202403')
img_extent = (22, 41, -5, 32)

#background map
ax3.imshow(nc, origin='upper',extent=img_extent, transform=ccrs.PlateCarree(),zorder=1)
# LAND_10m = cartopy.feature.NaturalEarthFeature('physical', 'land', '10m',edgecolor='face', facecolor=cartopy.feature.COLORS['land'])
# ax3.add_feature(LAND_10m)
#add Natural earth shape files
ax3.coastlines(resolution='50m', color='white', linewidth=1, zorder=3)
#ax3.add_feature(cfeature.BORDERS, linewidth=0.5, color='red', zorder=3)
ax3.add_feature(cfeature.COASTLINE, linewidth=0.5, color='red', zorder=3)
ax3.add_feature(cfeature.RIVERS, linewidth=0.5, color='blue', zorder=3)

#add custom shape files: ax.add_geometries(Reader(fname).geometries(),ccrs.PlateCarree(),facecolor='white', hatch='xxxx')
ax3.add_geometries(Reader(basin_shape).geometries(),ccrs.PlateCarree(),facecolor='none', edgecolor='black', linewidth=0.50,zorder=3) 
ax3.add_geometries(Reader(lakes_nile).geometries(),ccrs.PlateCarree(),facecolor='blue', edgecolor='blue', linewidth=0.50,zorder=3) 
ax3.add_geometries(Reader(disp_bound).geometries(),ccrs.PlateCarree(),facecolor='none', edgecolor='red', linewidth=0.50,linestyle='dashed',zorder=3) 
ax3.add_geometries(Reader(no_disp_bound).geometries(),ccrs.PlateCarree(),facecolor='none', edgecolor='red', linewidth=0.50,zorder=3) 

#data raster - see title
cs = ax3.imshow(smn3, origin='upper', extent=img_extent, transform=ccrs.PlateCarree(), vmin=min_all, vmax=max_all, cmap='YlGn', interpolation='None', zorder=2)
gl = ax3.gridlines(draw_labels=True,color = 'grey', linestyle = '--', linewidth = 0.5)
gl.top_labels = False
gl.right_labels = False
gl.xlocator = mticker.FixedLocator([25, 30, 35, 40])

plt.subplots_adjust(left=0.1, bottom=0.1, right=0.9,top=0.9, wspace=0.3, hspace=0.1)

# Add a colorbar axis at the bottom of the graph
cbar_ax = fig.add_axes([0.35, 0.02, 0.3, 0.02])

# Draw the colorbar
cbar=fig.colorbar(cs, cax=cbar_ax,orientation='horizontal')
cbar.set_label('NDVI index')

plt.show();
fig.savefig(work_dir+'/SMN_202401 _202403.png',bbox_inches = 'tight')

Create Animations

Using precipitation - radar/gauge 5 minute real-time accumulations over the Netherlands

The archive of gridded files of radar-derived 5 minute precipitation accumulations is corrected by rain gauge data. Radar data from 1500 m over the Netherlands and surrounding area measured by Dutch, Belgian, and German radars are corrected by available data from automatic rain gauges. Time interval is 5 minutes.

Starting with data from 21 July 2024 - 10.00 UTC onwards (for the next 3 hours), this dataset is created using improved algorithms. This includes correction for signal attenuation, correction for vertical variation of precipitation, correction for fast-moving showers and use of uncertainty information in merging data from multiple radars.

Data set name: nl_rdr_data_rtcor_5m_tar

Data access: https://dataplatform.knmi.nl/dataset/access/nl-rdr-data-rtcor-5m-tar-1-0

Review Buienradar archive: https://www.buienradar.nl/nederland/neerslag/buienradar-terugkijken/archief/202407211000 (local time)

Data has been downloaded and is locally available.

# create a list of files containing the rainfall data
rain_data = []
# Iterate directory
for file in os.listdir(work_dir):
    # check only h5 files
    if file.endswith('.h5'):
        rain_data.append(file)
print(rain_data)
print('Number of files = ',len( rain_data))

['RAD_NL25_RAC_RT_202407211000.h5', 'RAD_NL25_RAC_RT_202407211005.h5', 'RAD_NL25_RAC_RT_202407211010.h5', 'RAD_NL25_RAC_RT_202407211015.h5', 'RAD_NL25_RAC_RT_202407211020.h5', 'RAD_NL25_RAC_RT_202407211025.h5', 'RAD_NL25_RAC_RT_202407211030.h5', 'RAD_NL25_RAC_RT_202407211035.h5', 'RAD_NL25_RAC_RT_202407211040.h5', 'RAD_NL25_RAC_RT_202407211045.h5', 'RAD_NL25_RAC_RT_202407211050.h5', 'RAD_NL25_RAC_RT_202407211055.h5', 'RAD_NL25_RAC_RT_202407211100.h5', 'RAD_NL25_RAC_RT_202407211105.h5', 'RAD_NL25_RAC_RT_202407211110.h5', 'RAD_NL25_RAC_RT_202407211115.h5', 'RAD_NL25_RAC_RT_202407211120.h5', 'RAD_NL25_RAC_RT_202407211125.h5', 'RAD_NL25_RAC_RT_202407211130.h5', 'RAD_NL25_RAC_RT_202407211135.h5', 'RAD_NL25_RAC_RT_202407211140.h5', 'RAD_NL25_RAC_RT_202407211145.h5', 'RAD_NL25_RAC_RT_202407211150.h5', 'RAD_NL25_RAC_RT_202407211155.h5', 'RAD_NL25_RAC_RT_202407211200.h5', 'RAD_NL25_RAC_RT_202407211205.h5', 'RAD_NL25_RAC_RT_202407211210.h5', 'RAD_NL25_RAC_RT_202407211215.h5', 'RAD_NL25_RAC_RT_202407211220.h5', 'RAD_NL25_RAC_RT_202407211225.h5', 'RAD_NL25_RAC_RT_202407211230.h5', 'RAD_NL25_RAC_RT_202407211235.h5', 'RAD_NL25_RAC_RT_202407211240.h5', 'RAD_NL25_RAC_RT_202407211245.h5', 'RAD_NL25_RAC_RT_202407211250.h5', 'RAD_NL25_RAC_RT_202407211255.h5', 'RAD_NL25_RAC_RT_202407211300.h5', 'RAD_NL25_RAC_RT_202407211305.h5', 'RAD_NL25_RAC_RT_202407211310.h5', 'RAD_NL25_RAC_RT_202407211315.h5', 'RAD_NL25_RAC_RT_202407211320.h5', 'RAD_NL25_RAC_RT_202407211325.h5', 'RAD_NL25_RAC_RT_202407211330.h5', 'RAD_NL25_RAC_RT_202407211335.h5', 'RAD_NL25_RAC_RT_202407211340.h5', 'RAD_NL25_RAC_RT_202407211345.h5', 'RAD_NL25_RAC_RT_202407211350.h5', 'RAD_NL25_RAC_RT_202407211355.h5', 'RAD_NL25_RAC_RT_202407211400.h5', 'RAD_NL25_RAC_RT_202407211405.h5', 'RAD_NL25_RAC_RT_202407211410.h5', 'RAD_NL25_RAC_RT_202407211415.h5', 'RAD_NL25_RAC_RT_202407211420.h5', 'RAD_NL25_RAC_RT_202407211425.h5', 'RAD_NL25_RAC_RT_202407211430.h5', 'RAD_NL25_RAC_RT_202407211435.h5', 'RAD_NL25_RAC_RT_202407211440.h5', 'RAD_NL25_RAC_RT_202407211445.h5', 'RAD_NL25_RAC_RT_202407211450.h5', 'RAD_NL25_RAC_RT_202407211455.h5', 'RAD_NL25_RAC_RT_202407211500.h5', 'RAD_NL25_RAC_RT_202407211505.h5', 'RAD_NL25_RAC_RT_202407211510.h5', 'RAD_NL25_RAC_RT_202407211515.h5', 'RAD_NL25_RAC_RT_202407211520.h5', 'RAD_NL25_RAC_RT_202407211525.h5', 'RAD_NL25_RAC_RT_202407211530.h5', 'RAD_NL25_RAC_RT_202407211535.h5', 'RAD_NL25_RAC_RT_202407211540.h5', 'RAD_NL25_RAC_RT_202407211545.h5', 'RAD_NL25_RAC_RT_202407211550.h5', 'RAD_NL25_RAC_RT_202407211555.h5', 'RAD_NL25_RAC_RT_202407211600.h5', 'RAD_NL25_RAC_RT_202407211605.h5', 'RAD_NL25_RAC_RT_202407211610.h5', 'RAD_NL25_RAC_RT_202407211615.h5', 'RAD_NL25_RAC_RT_202407211620.h5', 'RAD_NL25_RAC_RT_202407211625.h5', 'RAD_NL25_RAC_RT_202407211630.h5', 'RAD_NL25_RAC_RT_202407211635.h5', 'RAD_NL25_RAC_RT_202407211640.h5', 'RAD_NL25_RAC_RT_202407211645.h5', 'RAD_NL25_RAC_RT_202407211650.h5', 'RAD_NL25_RAC_RT_202407211655.h5', 'RAD_NL25_RAC_RT_202407211700.h5', 'RAD_NL25_RAC_RT_202407211705.h5', 'RAD_NL25_RAC_RT_202407211710.h5', 'RAD_NL25_RAC_RT_202407211715.h5', 'RAD_NL25_RAC_RT_202407211720.h5', 'RAD_NL25_RAC_RT_202407211725.h5', 'RAD_NL25_RAC_RT_202407211730.h5', 'RAD_NL25_RAC_RT_202407211735.h5', 'RAD_NL25_RAC_RT_202407211740.h5', 'RAD_NL25_RAC_RT_202407211745.h5', 'RAD_NL25_RAC_RT_202407211750.h5', 'RAD_NL25_RAC_RT_202407211755.h5', 'RAD_NL25_RAC_RT_202407211800.h5', 'RAD_NL25_RAC_RT_202407211805.h5', 'RAD_NL25_RAC_RT_202407211810.h5', 'RAD_NL25_RAC_RT_202407211815.h5', 'RAD_NL25_RAC_RT_202407211820.h5', 'RAD_NL25_RAC_RT_202407211825.h5', 'RAD_NL25_RAC_RT_202407211830.h5', 'RAD_NL25_RAC_RT_202407211835.h5', 'RAD_NL25_RAC_RT_202407211840.h5', 'RAD_NL25_RAC_RT_202407211845.h5', 'RAD_NL25_RAC_RT_202407211850.h5', 'RAD_NL25_RAC_RT_202407211855.h5', 'RAD_NL25_RAC_RT_202407211900.h5', 'RAD_NL25_RAC_RT_202407211905.h5', 'RAD_NL25_RAC_RT_202407211910.h5', 'RAD_NL25_RAC_RT_202407211915.h5', 'RAD_NL25_RAC_RT_202407211920.h5', 'RAD_NL25_RAC_RT_202407211925.h5', 'RAD_NL25_RAC_RT_202407211930.h5', 'RAD_NL25_RAC_RT_202407211935.h5', 'RAD_NL25_RAC_RT_202407211940.h5', 'RAD_NL25_RAC_RT_202407211945.h5', 'RAD_NL25_RAC_RT_202407211950.h5', 'RAD_NL25_RAC_RT_202407211955.h5', 'RAD_NL25_RAC_RT_202407212000.h5', 'RAD_NL25_RAC_RT_202407212005.h5', 'RAD_NL25_RAC_RT_202407212010.h5', 'RAD_NL25_RAC_RT_202407212015.h5', 'RAD_NL25_RAC_RT_202407212020.h5', 'RAD_NL25_RAC_RT_202407212025.h5', 'RAD_NL25_RAC_RT_202407212030.h5', 'RAD_NL25_RAC_RT_202407212035.h5', 'RAD_NL25_RAC_RT_202407212040.h5', 'RAD_NL25_RAC_RT_202407212045.h5', 'RAD_NL25_RAC_RT_202407212050.h5', 'RAD_NL25_RAC_RT_202407212055.h5']
Number of files =  132

#create definition to open and import the *.h5 files
def openh5(filename):
    try:
        f = h5py.File(filename, 'r')
        f.close
        return f
    except:
        print('File {:s} cannot be opened'.format(filename))
        return []

image_size = None
band_data = []

for f in (rain_data):
    print ('now processed = ',f)
    f = openh5(work_dir + '/' + f)
    image = f['image1']['image_data'][()]
    
    band_data.append(image)
    if image_size is None:
        # image_size contains the number of pixels of the images
        image_size = [f['geographic'].attrs['geo_number_columns'][0], f['geographic'].attrs['geo_number_rows'][0]]

print(image_size)

now processed =  RAD_NL25_RAC_RT_202407211000.h5
now processed =  RAD_NL25_RAC_RT_202407211005.h5
now processed =  RAD_NL25_RAC_RT_202407211010.h5
now processed =  RAD_NL25_RAC_RT_202407211015.h5
now processed =  RAD_NL25_RAC_RT_202407211020.h5
now processed =  RAD_NL25_RAC_RT_202407211025.h5
now processed =  RAD_NL25_RAC_RT_202407211030.h5
now processed =  RAD_NL25_RAC_RT_202407211035.h5
now processed =  RAD_NL25_RAC_RT_202407211040.h5
now processed =  RAD_NL25_RAC_RT_202407211045.h5
now processed =  RAD_NL25_RAC_RT_202407211050.h5
now processed =  RAD_NL25_RAC_RT_202407211055.h5
now processed =  RAD_NL25_RAC_RT_202407211100.h5
now processed =  RAD_NL25_RAC_RT_202407211105.h5
now processed =  RAD_NL25_RAC_RT_202407211110.h5
now processed =  RAD_NL25_RAC_RT_202407211115.h5
now processed =  RAD_NL25_RAC_RT_202407211120.h5
now processed =  RAD_NL25_RAC_RT_202407211125.h5
now processed =  RAD_NL25_RAC_RT_202407211130.h5
now processed =  RAD_NL25_RAC_RT_202407211135.h5
now processed =  RAD_NL25_RAC_RT_202407211140.h5
now processed =  RAD_NL25_RAC_RT_202407211145.h5
now processed =  RAD_NL25_RAC_RT_202407211150.h5
now processed =  RAD_NL25_RAC_RT_202407211155.h5
now processed =  RAD_NL25_RAC_RT_202407211200.h5
now processed =  RAD_NL25_RAC_RT_202407211205.h5
now processed =  RAD_NL25_RAC_RT_202407211210.h5
now processed =  RAD_NL25_RAC_RT_202407211215.h5
now processed =  RAD_NL25_RAC_RT_202407211220.h5
now processed =  RAD_NL25_RAC_RT_202407211225.h5
now processed =  RAD_NL25_RAC_RT_202407211230.h5
now processed =  RAD_NL25_RAC_RT_202407211235.h5
now processed =  RAD_NL25_RAC_RT_202407211240.h5
now processed =  RAD_NL25_RAC_RT_202407211245.h5
now processed =  RAD_NL25_RAC_RT_202407211250.h5
now processed =  RAD_NL25_RAC_RT_202407211255.h5
now processed =  RAD_NL25_RAC_RT_202407211300.h5
now processed =  RAD_NL25_RAC_RT_202407211305.h5
now processed =  RAD_NL25_RAC_RT_202407211310.h5
now processed =  RAD_NL25_RAC_RT_202407211315.h5
now processed =  RAD_NL25_RAC_RT_202407211320.h5
now processed =  RAD_NL25_RAC_RT_202407211325.h5
now processed =  RAD_NL25_RAC_RT_202407211330.h5
now processed =  RAD_NL25_RAC_RT_202407211335.h5
now processed =  RAD_NL25_RAC_RT_202407211340.h5
now processed =  RAD_NL25_RAC_RT_202407211345.h5
now processed =  RAD_NL25_RAC_RT_202407211350.h5
now processed =  RAD_NL25_RAC_RT_202407211355.h5
now processed =  RAD_NL25_RAC_RT_202407211400.h5
now processed =  RAD_NL25_RAC_RT_202407211405.h5
now processed =  RAD_NL25_RAC_RT_202407211410.h5
now processed =  RAD_NL25_RAC_RT_202407211415.h5
now processed =  RAD_NL25_RAC_RT_202407211420.h5
now processed =  RAD_NL25_RAC_RT_202407211425.h5
now processed =  RAD_NL25_RAC_RT_202407211430.h5
now processed =  RAD_NL25_RAC_RT_202407211435.h5
now processed =  RAD_NL25_RAC_RT_202407211440.h5
now processed =  RAD_NL25_RAC_RT_202407211445.h5
now processed =  RAD_NL25_RAC_RT_202407211450.h5
now processed =  RAD_NL25_RAC_RT_202407211455.h5
now processed =  RAD_NL25_RAC_RT_202407211500.h5
now processed =  RAD_NL25_RAC_RT_202407211505.h5
now processed =  RAD_NL25_RAC_RT_202407211510.h5
now processed =  RAD_NL25_RAC_RT_202407211515.h5
now processed =  RAD_NL25_RAC_RT_202407211520.h5
now processed =  RAD_NL25_RAC_RT_202407211525.h5
now processed =  RAD_NL25_RAC_RT_202407211530.h5
now processed =  RAD_NL25_RAC_RT_202407211535.h5
now processed =  RAD_NL25_RAC_RT_202407211540.h5
now processed =  RAD_NL25_RAC_RT_202407211545.h5
now processed =  RAD_NL25_RAC_RT_202407211550.h5
now processed =  RAD_NL25_RAC_RT_202407211555.h5
now processed =  RAD_NL25_RAC_RT_202407211600.h5
now processed =  RAD_NL25_RAC_RT_202407211605.h5
now processed =  RAD_NL25_RAC_RT_202407211610.h5
now processed =  RAD_NL25_RAC_RT_202407211615.h5
now processed =  RAD_NL25_RAC_RT_202407211620.h5
now processed =  RAD_NL25_RAC_RT_202407211625.h5
now processed =  RAD_NL25_RAC_RT_202407211630.h5
now processed =  RAD_NL25_RAC_RT_202407211635.h5
now processed =  RAD_NL25_RAC_RT_202407211640.h5
now processed =  RAD_NL25_RAC_RT_202407211645.h5
now processed =  RAD_NL25_RAC_RT_202407211650.h5
now processed =  RAD_NL25_RAC_RT_202407211655.h5
now processed =  RAD_NL25_RAC_RT_202407211700.h5
now processed =  RAD_NL25_RAC_RT_202407211705.h5
now processed =  RAD_NL25_RAC_RT_202407211710.h5
now processed =  RAD_NL25_RAC_RT_202407211715.h5
now processed =  RAD_NL25_RAC_RT_202407211720.h5
now processed =  RAD_NL25_RAC_RT_202407211725.h5
now processed =  RAD_NL25_RAC_RT_202407211730.h5
now processed =  RAD_NL25_RAC_RT_202407211735.h5
now processed =  RAD_NL25_RAC_RT_202407211740.h5
now processed =  RAD_NL25_RAC_RT_202407211745.h5
now processed =  RAD_NL25_RAC_RT_202407211750.h5
now processed =  RAD_NL25_RAC_RT_202407211755.h5
now processed =  RAD_NL25_RAC_RT_202407211800.h5
now processed =  RAD_NL25_RAC_RT_202407211805.h5
now processed =  RAD_NL25_RAC_RT_202407211810.h5
now processed =  RAD_NL25_RAC_RT_202407211815.h5
now processed =  RAD_NL25_RAC_RT_202407211820.h5
now processed =  RAD_NL25_RAC_RT_202407211825.h5
now processed =  RAD_NL25_RAC_RT_202407211830.h5
now processed =  RAD_NL25_RAC_RT_202407211835.h5
now processed =  RAD_NL25_RAC_RT_202407211840.h5
now processed =  RAD_NL25_RAC_RT_202407211845.h5
now processed =  RAD_NL25_RAC_RT_202407211850.h5
now processed =  RAD_NL25_RAC_RT_202407211855.h5
now processed =  RAD_NL25_RAC_RT_202407211900.h5
now processed =  RAD_NL25_RAC_RT_202407211905.h5
now processed =  RAD_NL25_RAC_RT_202407211910.h5
now processed =  RAD_NL25_RAC_RT_202407211915.h5
now processed =  RAD_NL25_RAC_RT_202407211920.h5
now processed =  RAD_NL25_RAC_RT_202407211925.h5
now processed =  RAD_NL25_RAC_RT_202407211930.h5
now processed =  RAD_NL25_RAC_RT_202407211935.h5
now processed =  RAD_NL25_RAC_RT_202407211940.h5
now processed =  RAD_NL25_RAC_RT_202407211945.h5
now processed =  RAD_NL25_RAC_RT_202407211950.h5
now processed =  RAD_NL25_RAC_RT_202407211955.h5
now processed =  RAD_NL25_RAC_RT_202407212000.h5
now processed =  RAD_NL25_RAC_RT_202407212005.h5
now processed =  RAD_NL25_RAC_RT_202407212010.h5
now processed =  RAD_NL25_RAC_RT_202407212015.h5
now processed =  RAD_NL25_RAC_RT_202407212020.h5
now processed =  RAD_NL25_RAC_RT_202407212025.h5
now processed =  RAD_NL25_RAC_RT_202407212030.h5
now processed =  RAD_NL25_RAC_RT_202407212035.h5
now processed =  RAD_NL25_RAC_RT_202407212040.h5
now processed =  RAD_NL25_RAC_RT_202407212045.h5
now processed =  RAD_NL25_RAC_RT_202407212050.h5
now processed =  RAD_NL25_RAC_RT_202407212055.h5
[np.int32(700), np.int32(765)]

#check results of data import
print(image_size[0], image_size[1])
print(len(band_data))

700 765
132

#create animation of raw imported data

imageList = band_data
def getImageFromList(x):
    return imageList[x]

fig = plt.figure(figsize=(8, 8))
ims = []
for i in range(len(imageList)):
    im = plt.imshow(getImageFromList(i), vmin='0', vmax='155', cmap='jet', animated=True)
    ims.append([im])
    
ani = animation.ArtistAnimation(fig, ims, interval=400, blit=True, repeat=True)

plt.xticks(color='w')
plt.yticks(color='w')

plt.grid(True, alpha = 0.5)
plt.close()

# Show the animation inline (below)
HTML(ani.to_html5_video())

# Uncomment line below to save animation as video 
#ani.save(work_dir+'/ani_raw.mp4', fps=2)

Create a georeference and assign this to an empty raster coverage and fill it with the data imported before and review the projection and dimensions (x, y, z) of the new ilwis map list (data stack) created

# from meta data information provided: +proj=stere +lat_0=90 +lon_0=0 +lat_ts=60 +a=6378.14 +b=6356.75 +x_0=0 y_0=0
grf = ilwis.GeoReference('code=georef:type=corners,csy=proj4:+proj=stere +lat_0=60 +lon_0=0 +lat_ts=60 +a=6378137 +b=6356750 +x_0=0 +y_0=90 +datum=WGS84 +units=m, envelope=0 -448497.662 686614.850 -1196791.122, gridsize=700 765,cornerofcorners=yes')

#create empty ilwis raster 
pcpAll = ilwis.RasterCoverage()
defNumr = ilwis.DataDefinition(ilwis.NumericDomain('code=value'), ilwis.NumericRange(0, 1000, 0))
pcpAll.setDataDef(defNumr)
pcpAll.setSize(ilwis.Size(int(image_size[0]), int(image_size[1]), len(rain_data)))
pcpAll.setGeoReference(grf)

pcpAll.array2raster(np.array(band_data).flatten())
print(pcpAll.size())

Size(700, 765, 132)

#data check
print(pcpAll.size())
print(pcpAll.envelope())
pcpAllcoordSys = pcpAll.coordinateSystem()
print(pcpAllcoordSys.toWKT())
print(pcpAllcoordSys.toProj4())

Size(700, 765, 132)
0.000000 -1196791.122000 686614.850000 -448497.662000
PROJCS["_ANONYMOUS_1005288",GEOCS["_ANONYMOUS_1005288",DATUM[" WGS 84",[DWGS84],ELLIPSOID["User Defined",6378137.000000000000,298.224949735821],PRIMEM["Greenwich",0, AUTHORITY["EPSG",8901"]]],PROJECTION["Polar_Stereographic"],PARAMETER[false_easting,0],PARAMETER[false_northing,90],PARAMETER[central_meridian,0],PARAMETER[latitude_of_origin,60],UNIT[meter,1.0]]
+proj=stere +a=6378137 +b=6356750 +datum=WGS84 +lat_0=60 +lon_0=0 +x_0=0 +y_0=90

#uncomment line below if you want to store results of the raw data as an ILWIS maplist
#pcpAll.store('pcp_raw.mpl')

#apply calibration formulas = GEO = 0.010000*PV+0.000000
pcp_acc5m = ilwis.do('mapcalc', '0.01 * @1', pcpAll)

#aggregation to get total over time period
R_all = ilwis.do('aggregaterasterstatistics', pcp_acc5m, 'sum') # add all 5 minutes accumumulated timesteps
R_all_fin = ilwis.do('mapcalc', 'iff(@1<655, @1, -999)', R_all) # set background to -999

R_all_fin.store('pcp_accum_gauge_cor_stgr.mpr')

Temporary issue

To display the map in ILWIS386, the coordinate system created 'pcp_accum_gauge_cor_stgr.csy' has to be adapted. The code field below changes 'Projection=?' into 'Projection=StereoGraphic'. Now you can display the map without an error message.

#input file
fin = open(work_dir+"/pcp_accum_gauge_cor_stgr.csy", "rt")
#output file to write the result to
fout = open(work_dir+"/temp.csy", "wt")
#for each line in the input file
for line in fin:
    #read replace the string and write to output file
    fout.write(line.replace('Projection=?', 'Projection=StereoGraphic'))
#close input and output files
fin.close()
fout.close()

#remove file
if os.path.exists(work_dir+"/pcp_accum_gauge_cor_stgr.csy"):
  os.remove(work_dir+"/pcp_accum_gauge_cor_stgr.csy")
else:
  print("The file does not exist")


old_name = work_dir+"/temp.csy"
new_name = work_dir+"/pcp_accum_gauge_cor_stgr.csy"
# Renaming the file
os.rename(old_name, new_name)

#read raster as np-array and transform - reshape from 1D to 2D
R_all_fin_im = np.fromiter(iter(R_all_fin), np.float64, R_all_fin.size().linearSize())
R_all_fin_im = R_all_fin_im.reshape((R_all_fin.size().ysize, R_all_fin.size().xsize ))

#nodata is assigned to nan - masked 
R_all_fin_im_mask = np.ma.masked_where(R_all_fin_im == -999, R_all_fin_im)

# display the accumulated precipitation map
plt.figure(figsize = (4,8))
plt.imshow(R_all_fin_im_mask, vmin=0, vmax =100, interpolation='none', cmap = 'jet')
plt.title('3 hr accumulated precipitation')
plt.axis('off')
plt.show();

Resample the map to UTM - zone 32N

#transform from stereographic to UTM projection, here zone 32 -'EPSG:32632' 
#create target georeference
grf_utm = ilwis.GeoReference('code=georef:type=corners, csy=epsg:32632, envelope= 84512.281 5591939.637 401512.281 5954939.637, gridsize=317 363, cornerofcorners=yes')

print(grf_utm.coordinateSystem())

WGS 84 / UTM zone 32N

#conduct the resampling 
R_all_res = ilwis.do('resample', pcp_acc5m, grf_utm, 'nearestneighbour')

#save the results
R_all_res.store(work_dir+'/gauge_cor_radar_sum_utm.mpl') #ilwis maplist
R_all_res.store(work_dir+'/gauge_cor_radar_sum_utm.tif', 'GTiff', 'gdal')#for use in QGIS - now as geotif

Retrieve precipitation time series values for certain location

First we specify a location (using line and column position) and extract the time series values, then we create the time annotation to be used on the x-axis of the plot. Finally the plot is created

R_ts5 = ilwis.do('mapcalc', 'iff(@1<655, @1, ?)',R_all_res) # set background to undefined (representing no data)

#city centre of Enschede - the Netherlands
l = 168
c = 272

z = (R_ts5.size().zsize)
print(z)

PCP_values = []
for n in range(0,(z)):
    point_value = (R_ts5.pix2value(ilwis.Pixel(c,l,n)))
    PCP_values.append(point_value)

print('Values extracted for selected location:', PCP_values)

132
Values extracted for selected location: [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.06, 0.0, 0.0, 0.0, 0.04, 0.02, 0.02, 0.0, 0.0, 0.0, 0.1, 0.18, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.02, 0.03, 0.07, 0.04, 0.07, 0.05, 0.01, 0.02, 0.04, 0.06, 0.07, 0.07, 0.09, 0.13, 0.19, 0.58, 0.54, 0.39, 0.32, 0.32, 0.37, 0.58, 1.29, 3.17, 3.92, 3.36, 3.15, 2.36, 1.3, 0.66, 0.34, 0.28, 0.42, 0.68, 0.72, 0.52, 0.36, 0.26, 0.18, 0.09, 0.06, 0.05, 0.03, 0.02, 0.02, 0.0, 0.0, 0.0, 0.0, 0.02, 0.02, 0.02, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.02, 0.03, 0.03, 0.03, 0.02, 0.0, 0.0, 0.0, 0.0, 0.0, 0.02, 0.02, 0.02, 0.02, 0.04, 0.06, 0.04, 0.08, 0.1, 0.16, 0.13, 0.1, 0.08, 0.03, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]

# initializing time range for X-Axis
plot_date = datetime.strptime("10:00", "%H:%M")

z = (R_ts5.size().zsize)

current_time = plot_date
end_time = plot_date + timedelta(hours=11) 
interval = timedelta(minutes=5)

date_ts = []
while current_time < end_time:
    time_ts = (current_time.strftime("%H:%M"))
    current_time += interval
    date_ts.append(time_ts)

print(date_ts)

['10:00', '10:05', '10:10', '10:15', '10:20', '10:25', '10:30', '10:35', '10:40', '10:45', '10:50', '10:55', '11:00', '11:05', '11:10', '11:15', '11:20', '11:25', '11:30', '11:35', '11:40', '11:45', '11:50', '11:55', '12:00', '12:05', '12:10', '12:15', '12:20', '12:25', '12:30', '12:35', '12:40', '12:45', '12:50', '12:55', '13:00', '13:05', '13:10', '13:15', '13:20', '13:25', '13:30', '13:35', '13:40', '13:45', '13:50', '13:55', '14:00', '14:05', '14:10', '14:15', '14:20', '14:25', '14:30', '14:35', '14:40', '14:45', '14:50', '14:55', '15:00', '15:05', '15:10', '15:15', '15:20', '15:25', '15:30', '15:35', '15:40', '15:45', '15:50', '15:55', '16:00', '16:05', '16:10', '16:15', '16:20', '16:25', '16:30', '16:35', '16:40', '16:45', '16:50', '16:55', '17:00', '17:05', '17:10', '17:15', '17:20', '17:25', '17:30', '17:35', '17:40', '17:45', '17:50', '17:55', '18:00', '18:05', '18:10', '18:15', '18:20', '18:25', '18:30', '18:35', '18:40', '18:45', '18:50', '18:55', '19:00', '19:05', '19:10', '19:15', '19:20', '19:25', '19:30', '19:35', '19:40', '19:45', '19:50', '19:55', '20:00', '20:05', '20:10', '20:15', '20:20', '20:25', '20:30', '20:35', '20:40', '20:45', '20:50', '20:55']

fig = plt.figure(figsize =(8, 5))
plt.plot(date_ts, PCP_values, label='Rainfall')
plt.ylim(0, 4)#set value range for y-axis
frequency = 12
plt.xticks(date_ts[::frequency])#plot only 15 minutes time step on x-axis
plt.xlabel('Time period selected')
plt.ylabel('Accumulated rainfall (mm/5-min)')
plt.title('Rainfall extracted for Enschede city centre')
plt.legend();

Create animation of resulting utm resampled data

#transform to numpy array
R_ts5_2np = np.fromiter(iter(R_ts5), np.float64, R_ts5.size().linearSize()) 
R_ts5_2np = R_ts5_2np.reshape((R_ts5.size().zsize, R_ts5.size().ysize, R_ts5.size().xsize))

#check the data
R_ts5_2np.shape

(132, 363, 317)

R_ts5_2np = np.ma.masked_where(R_ts5_2np < 0.001, R_ts5_2np) #set background to white for animation

imageList = R_ts5_2np
def getImageFromList(x):
    return imageList[x]

fig = plt.figure(figsize=(8,4))
ims = []
for i in range(len(imageList)):
    im = plt.imshow(getImageFromList(i), vmin='0', vmax='4', cmap='jet', animated=True)
    ims.append([im])
    
ani = animation.ArtistAnimation(fig, ims, interval=400, blit=True, repeat=True)

plt.xticks(color='w')
plt.yticks(color='w')
cbar = plt.colorbar(shrink=0.65)
cbar.set_label('pcp (mm/5 min)')

plt.grid(True, alpha = 0.5)
plt.close()

# Uncomment line below to save animation as video 
#ani.save(work_dir+'/ani.mp4', fps=2)

# Show the animation inline (below)
HTML(ani.to_html5_video())

Now also using cartopy

Display first layer of time series Using Cartopy we set the projection to UTM, zone 32, with the plot extent of our map. Retrieve OpenStreetMap image tiles and then the coastline and country borders are added as well as the first layer of our array 'R_ts5_2np[0]' and this is plotted using 'imshow'

warnings.filterwarnings('ignore')

fig = plt.figure(figsize=(10, 6))
subplot_extent = [84512, 401512, 5591939, 5954939]

request = cimgt.OSM()
#request = cimgt.QuadtreeTiles()
#request = cimgt.GoogleTiles(style="satellite") #(style="satellite") (style="streets") (style="only_streets") (style="terrain")
projection = ccrs.UTM(32)

ax = fig.add_subplot(projection=projection)
ax.set_extent(subplot_extent, projection)


ax.coastlines(linewidth=0.5, color='blue')
ax.add_feature(cfeature.BORDERS, linewidth=0.5, color='red')
plt.imshow(R_ts5_2np[0], extent=subplot_extent,  interpolation = 'None', vmin='0', vmax='4', cmap='jet', alpha=1) #first rainfall event
ax.add_image(request, 8, alpha=0.4)
ax = plt.subplot(1, 1, 1, projection=ccrs.UTM(32))

cbar = plt.colorbar(shrink=0.65)
cbar.set_label('PCP in mm accumulated', labelpad=5, y=0.5, rotation=90)
ax.set_title("Rainfall timestep-1, UTM projection");

Final animation

With the projection and additional layers displayed correctly we can now create the animation of the rainfall time series, we also include the time stamp for every time step. Use is made of 'datetime' to define a list of time steps which can be included in the annotation when the animation is created

start_time = datetime(2024, 7, 21, 10, 00) # Replace this with your desired start time
end_time = start_time + timedelta(hours=11)  # Replace this with your desired end time

current_time = start_time
interval = timedelta(minutes=5)

ts = []
while current_time < end_time:
    time_ani = (current_time.strftime("%H:%M"))
    current_time += interval
    ts.append(time_ani)
print(ts)

['10:00', '10:05', '10:10', '10:15', '10:20', '10:25', '10:30', '10:35', '10:40', '10:45', '10:50', '10:55', '11:00', '11:05', '11:10', '11:15', '11:20', '11:25', '11:30', '11:35', '11:40', '11:45', '11:50', '11:55', '12:00', '12:05', '12:10', '12:15', '12:20', '12:25', '12:30', '12:35', '12:40', '12:45', '12:50', '12:55', '13:00', '13:05', '13:10', '13:15', '13:20', '13:25', '13:30', '13:35', '13:40', '13:45', '13:50', '13:55', '14:00', '14:05', '14:10', '14:15', '14:20', '14:25', '14:30', '14:35', '14:40', '14:45', '14:50', '14:55', '15:00', '15:05', '15:10', '15:15', '15:20', '15:25', '15:30', '15:35', '15:40', '15:45', '15:50', '15:55', '16:00', '16:05', '16:10', '16:15', '16:20', '16:25', '16:30', '16:35', '16:40', '16:45', '16:50', '16:55', '17:00', '17:05', '17:10', '17:15', '17:20', '17:25', '17:30', '17:35', '17:40', '17:45', '17:50', '17:55', '18:00', '18:05', '18:10', '18:15', '18:20', '18:25', '18:30', '18:35', '18:40', '18:45', '18:50', '18:55', '19:00', '19:05', '19:10', '19:15', '19:20', '19:25', '19:30', '19:35', '19:40', '19:45', '19:50', '19:55', '20:00', '20:05', '20:10', '20:15', '20:20', '20:25', '20:30', '20:35', '20:40', '20:45', '20:50', '20:55']

warnings.filterwarnings('ignore')

imageList = R_ts5_2np

fig = plt.figure(figsize=(10, 6))
subplot_extent = [84512.281, 401512.281, 5591939.637, 5954939.637]

#retrieve openstreetmap using request
request = cimgt.OSM()
projection = ccrs.UTM(32)


ax = plt.axes(projection=projection)
ax.set_extent(subplot_extent, projection)


ax.coastlines(color='blue', linewidth=0.5, alpha=0.3)
#ax.add_feature(cfeature.BORDERS, linewidth=0.5, color='red', alpha=0.3)
ax.add_image(request, 8, alpha=0.4)

ims = []
for i in range(len(imageList)):
    im = plt.imshow(imageList[i], extent=subplot_extent, interpolation = 'None', vmin='0', vmax='4', cmap='jet', animated=True, alpha=1)
    #ann = ax.annotate(f"Time step {i:.0f}.", (1.0, 1.01), xycoords="axes fraction", ha="right")
    ann = ax.annotate (f'20240721 - {ts[i]}',(1.0, 1.01), xycoords="axes fraction", ha="right")
    ims.append([im, ann])
    
ani_fin = animation.ArtistAnimation(fig, ims, interval=400, blit=True, repeat=True)


ax = plt.subplot(1, 1, 1, projection=projection)
ax.set_title("Final rainfall animation", loc='left');
cbar = plt.colorbar(shrink=0.65)
cbar.set_label('PCP in mm accumulated', labelpad=5, y=0.5, rotation=90)
#clb.ax.set_title('Accum-PCP',x=1.07) 
#clb.ax.set_ylabel('(mm/5min)', labelpad=20, rotation=90) 

plt.close()

# Show the animation inline (below)
HTML(ani_fin.to_html5_video())

# Uncomment line below to save animation as video 
#ani_fin.save(work_dir+'/ani_fin.mp4', fps=2)