Hi GEONETCasters!

Yesterday we have seen a nice presentation of the Flood Mapping Products now available on GNC-A.

Please find below preliminary Python scripts, that should get you started.

In the next version of SHOWCast, the Flood Mapping Products will be supported by default.

READING THE FLOOD MAP GEOTIFF’s AS PYTHON ARRAYS AND MAKING BASIC MANIPULATION

This opens many possibilities!

#######################################################################################################
# LICENSE
# Copyright (C) 2020 - INPE - NATIONAL INSTITUTE FOR SPACE RESEARCH - BRAZIL
# This program is free software: you can redistribute it and/or modify it under the terms of the GNU
# General Public License as published by the Free Software Foundation, either version 3 of the License,
# or (at your option) any later version.
# This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without
# even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
# Public License for more details.
# You should have received a copy of the GNU General Public License along with this program.
# If not, see http://www.gnu.org/licenses/.
#######################################################################################################
#------------------------------------------------------------------------------------------------------
#------------------------------------------------------------------------------------------------------
# Required modules
#--------------------------------
#to run in a pure text terminal:
import matplotlib
matplotlib.use('Agg')
#--------------------------------
import matplotlib.pyplot as plt                              # Plotting library
import matplotlib.colors                                     # Matplotlib colors
import numpy as np                                           # Scientific computing with Python
import cartopy, cartopy.crs as ccrs                          # Plot maps
import cartopy.io.shapereader as shpreader                   # Import shapefiles
import time as t                                             # Time access and conversion
import os                                                    # Miscellaneous operating system interfaces
import sys                                                   # Import the "system specific parameters and functions" module
from mpl_toolkits.axes_grid1.inset_locator import inset_axes # Add a child inset axes to this existing axes.
from matplotlib.colors import LinearSegmentedColormap        # Linear interpolation for color maps
from datetime import datetime, timedelta                     # Library to convert julian day to dd-mm-yyyy
from osgeo import gdal, osr, ogr                             # Import GDAL
#------------------------------------------------------------------------------------------------------
#------------------------------------------------------------------------------------------------------
# Start the counter
print('Script started.')
start = t.time()

# GeoTIFF file path
path = ("C://VLAB//Samples//Flood-Composite//RIVER-FLDglobal-composite_20200511.part032.tif")
#path = ("C://VLAB//Samples//Flood-Joint//RIVER-FLD-joint-ABI_20200511.part032.tif")
#path = ("C://VLAB//Samples//Flood-ABI//River-Flood-ABI-hourly_20200513_140000.part004.force.tif")

# Get the file name form path
file_name = os.path.splitext(os.path.basename(path))[0]
#------------------------------------------------------------------------------------------------------
#------------------------------------------------------------------------------------------------------

# Get the product type from the file name
if ("ABI" in file_name):
product = "GOES-16 ABI Hourly Composite"
hour = ''.join(filter(str.isdigit, file_name))[8:14]
print(hour)
if ("global" in file_name): product = "Suomi-NPP and NOAA-20 VIIRS 5-day Composite (daily)"
if ("joint" in file_name): product = "Joint VIIRS / ABI flood product (daily)"
print ("Product:", product)

# Get the AOI/Full Disk part, from the file name
part = ''.join(filter(str.isdigit, file_name))[-3:]
print("Part:", part)

# Get the time and date from the file name
date = ''.join(filter(str.isdigit, file_name))[0:8]
print("Date:", date)

#------------------------------------------------------------------------------------------------------
#------------------------------------------------------------------------------------------------------

# Read the GeoTIFF as array
file = gdal.Open(path)
data = file.GetRasterBand(1).ReadAsArray()
print("Size:", data.shape)

# Getting the GeoTIFF extent
geoTransform = file.GetGeoTransform()
min_lon = geoTransform[0]
max_lon = min_lon + geoTransform[1] * file.RasterXSize
max_lat = geoTransform[3]
min_lat = max_lat + geoTransform[5] * file.RasterYSize
# The full GeoTIFF extent
extent = [min_lon, min_lat, max_lon, max_lat]
print("Min Lon:", min_lon, "Min Lat:", min_lat, "Max Lon:", max_lon, "Max Lat:", max_lat)

# The extent I would like to plot
plot_area = [-53.4, -2.04, -51.7, -0.7]
#------------------------------------------------------------------------------------------------------
#------------------------------------------------------------------------------------------------------

# Choose the plot size (width x height, in inches)
dpi = 150
fig = plt.figure(figsize=(data.shape[1]/dpi, data.shape[0]/dpi), dpi=dpi)

# Use the Geostationary projection in cartopy
ax = plt.axes([0, 0, 1, 1], projection=ccrs.PlateCarree())
# The full extent of the GeoTIFF
img_extent = [extent[0], extent[2], extent[1], extent[3]]
# The extent you want to plot (the limit is the full extent above)
ax.set_extent([plot_area[0], plot_area[2], plot_area[1], plot_area[3]], ccrs.PlateCarree())

#------------------------------------------------------------------------------------------------------
#------------------------------------------------------------------------------------------------------

# COLOR SETUP 1
# Define plot colors [VIIRS | Hourly ABI]
c_nodata = '#000000' # 14       (No Data)
c_land   = '#c4a272' # 17 | 43  (Land)
c_clouds = '#c8c8c8' # 30 | 99  (Clouds)
c_supra  = '#b400e6' # 56       (Supra Snow/Ice Water)
c_snow   = '#ffffff' # 70       (Snow)
c_ice    = '#00ffff' # 84       (Ice)
c_water  = '#0000ff' # 99 | 126 (Normal Open Water)
#Floodwater          # 100~200 | 140~254
c_floodwater = ["#3dfaab", "#2df863", "#00f700", "#c6f800", "#faf690", "#f7f701", "#f5c400", "#f79029", "#f75e00", "#ff0000"]
cmap_flood = matplotlib.colors.ListedColormap(c_floodwater)

'''
# COLOR SETUP 2
# Define plot colors [VIIRS | Hourly ABI]
c_nodata = '#000000' # 14       (No Data)
c_land   = '#282828' # 17 | 43  (Land)
c_clouds = '#383838' # 30 | 99  (Clouds)
c_supra  = '#b400e6' # 56       (Supra Snow/Ice Water)
c_snow   = '#ffffff' # 70       (Snow)
c_ice    = '#00ffff' # 84       (Ice)
c_water  = '#0000ff' # 99 | 126 (Normal Oper Water)
#Floodwater          # 100~200 | 140~254
c_floodwater = ["#3dfaab", "#2df863", "#00f700", "#c6f800", "#faf690", "#f7f701", "#f5c400", "#f79029", "#f75e00", "#ff0000"]
cmap_flood = matplotlib.colors.ListedColormap(c_floodwater)
'''
'''
# COLOR SETUP 3
# Define plot colors [VIIRS | Hourly ABI]
c_nodata = '#000000' # 14       (No Data)
c_land   = '#013220' # 17 | 43  (Land)
c_clouds = '#383838' # 30 | 99  (Clouds)
c_supra  = '#b400e6' # 56       (Supra Snow/Ice Water)
c_snow   = '#ffffff' # 70       (Snow)
c_ice    = '#00ffff' # 84       (Ice)
c_water  = '#000033' # 99 | 126 (Normal Oper Water)
#Floodwater          # 100~200 | 140~254
c_floodwater = ["#3dfaab", "#2df863", "#00f700", "#c6f800", "#faf690", "#f7f701", "#f5c400", "#f79029", "#f75e00", "#ff0000"]
cmap_flood = 'Blues'
'''

print("Plotting the Map Elements...")

# Plot the map elements (land, clouds, snow, ice, etc)
if (product == "GOES-16 ABI Hourly Composite"):
cmap = matplotlib.colors.ListedColormap([c_nodata, c_land, c_supra, c_snow, c_ice, c_clouds, c_water])
boundaries = [0, 15, 44, 57, 71, 85, 100, 127]
# Get only the flood from the full data
data_flood = data.astype(np.float64)
data_flood[data_flood < 142] = np.nan
data_flood = data_flood - 140
vmin = 0
vmax = 104
else:
cmap = matplotlib.colors.ListedColormap([c_nodata, c_land, c_clouds, c_supra, c_snow, c_ice, c_water])
boundaries = [0, 15, 18, 31, 57, 71, 85, 100]
# Get only the flood from the full data
data_flood = data.astype(np.float64)
data_flood[data_flood < 100] = np.nan
data_flood = data_flood - 100
vmin = 0
vmax = 100

norm = matplotlib.colors.BoundaryNorm(boundaries, cmap.N, clip=True)
img1 = ax.imshow(data, origin='upper', extent=img_extent, cmap=cmap, norm=norm, zorder=1)

print("Plotting the Floodwater fraction...")

# Plot the floodwater fraction (%)
img2 = ax.imshow(data_flood, vmin=vmin, vmax=vmax, origin='upper', extent=img_extent, cmap=cmap_flood, zorder=2)

#------------------------------------------------------------------------------------------------------
#------------------------------------------------------------------------------------------------------

print("Plotting other elements...")

# Required to put the colorbar inside plot
axins1 = inset_axes(ax, width="100%", height="1%", loc='lower center', borderpad=0.0)

# Add states and provinces
shapefile = list(shpreader.Reader('..//Shapefiles//ne_10m_admin_1_states_provinces.shp').geometries())
ax.add_geometries(shapefile, ccrs.PlateCarree(), edgecolor='white',facecolor='none', linewidth=2.00, zorder=3)
# Add countries
shapefile = list(shpreader.Reader('..//Shapefiles//ne_50m_admin_0_countries.shp').geometries())
ax.add_geometries(shapefile, ccrs.PlateCarree(), edgecolor='black',facecolor='none', linewidth=2.25, zorder=4)
# Add continents
shapefile = list(shpreader.Reader('..//Shapefiles//ne_10m_coastline.shp').geometries())
ax.add_geometries(shapefile, ccrs.PlateCarree(), edgecolor='black',facecolor='none', linewidth=2.50, zorder=5)
# Add coastlines, borders and gridlines
ax.gridlines(color='white', alpha=0.5, linestyle='--', linewidth=2.0, xlocs=np.arange(-180, 180, 0.5), ylocs=np.arange(-180, 180, 0.5), draw_labels=False, zorder=6)

# Put the floodwater fraction (%) colorbar inside plot
ticks = np.arange(0, 100, 10).tolist()
ticks = 10 * np.round(np.true_divide(ticks,10))
ticks = ticks[1:]
cb = fig.colorbar(img2, cax=axins1, orientation="horizontal", ticks=ticks)
cb.outline.set_visible(False)
cb.ax.tick_params(width = 0)
cb.ax.xaxis.set_tick_params(pad=0)
cb.ax.xaxis.set_ticks_position('top')
cb.ax.tick_params(axis='x', colors='black', labelsize=int(data.shape[0] * 0.005))

# Add a title
#plt.annotate(product + " " + date , xy=(int(data.shape[1] * 0.01), data.shape[0] - int(data.shape[0] * 0.016)), xycoords='figure pixels', fontsize=int(data.shape[1] * 0.005), fontweight='bold', color='white', bbox=dict(boxstyle="round",fc=(0.0, 0.0, 0.0), ec=(1., 1., 1.)), zorder=8)

#------------------------------------------------------------------------------------------------------
#------------------------------------------------------------------------------------------------------

print("Saving...")

# Save the figure
plt.savefig('flood_map_example.png', bbox_inches='tight', pad_inches=0, facecolor='white')

# Finish the counter
print('Finished. Total processing time:', round((t.time() - start),2), 'seconds.')

This plot was produced with the “color setup 2” from the example script. It uses dark gray tones for land and clouds, and the same colors for rivers and floodwater fractions and the previous example. It is easy to change the plot colors on the script!

This plot was produced with the “color setup 3” from the example script. It uses dark green for lan, dark gray for clouds, dark blue for rivers and bluish tones for floodwater fraction. It is easy to change the plot colors on the script!

Please download the script by clicking at the image below.

CONVERTING TO KMZ, MOSAICKING AND SUBSECTING SUB REGIONS

#######################################################################################################
# LICENSE
# Copyright (C) 2020 - INPE - NATIONAL INSTITUTE FOR SPACE RESEARCH - BRAZIL
# This program is free software: you can redistribute it and/or modify it under the terms of the GNU
# General Public License as published by the Free Software Foundation, either version 3 of the License,
# or (at your option) any later version.
# This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without
# even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
# Public License for more details.
# You should have received a copy of the GNU General Public License along with this program.
# If not, see http://www.gnu.org/licenses/.
#######################################################################################################
#------------------------------------------------------------------------------------------------------
#------------------------------------------------------------------------------------------------------
# Required modules
#--------------------------------
from osgeo import gdal, osr, ogr                             # Import GDAL
from datetime import datetime, timedelta                     # Library to convert julian day to dd-mm-yyyy
import time as t                                             # Time access and conversion
import os                                                    # Miscellaneous operating system interfaces
import sys                                                   # Import the "system specific parameters and functions" module
import glob                                                  # Unix style pathname pattern expansion
#------------------------------------------------------------------------------------------------------
#------------------------------------------------------------------------------------------------------
# Start the time counter
print('Script started.')
start = t.time()

#------------------------------------------------------------------------------------------------------
#------------------------------------------------------------------------------------------------------
# CONVERTION TO KMZ
print("Example 1: Converting an AOI to KMZ...")

# GeoTIFF file path
path = ("C://VLAB//Samples//Flood-Composite//RIVER-FLDglobal-composite_20200511.part032.tif")

# Converting from GeoTIFF to KMZ (Google Earth)
translate_options = gdal.TranslateOptions(rgbExpand='RGB')
gdal.Translate(path + '.kmz', path, format='KMLSUPEROVERLAY', options=translate_options)
#------------------------------------------------------------------------------------------------------
#------------------------------------------------------------------------------------------------------
# MOSAICKING AOI's AND CONVERTING TO KMZ

print("Example 2: Mosaicking AOI's and converting to KMZ...")

# Mosaicking diferent AOI's (in this example parts 20 to 29 from May 11th):
gdal.BuildVRT('mosaic_1.vrt', sorted(glob.glob( 'C://VLAB//Samples//Flood-Composite//*20200511*part02*.tif'), key=os.path.getmtime))
gdal.Translate('mosaic_1.tif',  'mosaic_1.vrt')

# Converting the Mosaic from GeoTIFF to KMZ (Google Earth)
translate_options = gdal.TranslateOptions(rgbExpand='RGB')
gdal.Translate('mosaic_1.kmz', 'mosaic_1.tif', format='KMLSUPEROVERLAY', options=translate_options)
#------------------------------------------------------------------------------------------------------
#------------------------------------------------------------------------------------------------------
# MOSAICKING AOI's AND CUTTING A REGION

print("Example 3: Mosaicking AOI's and cutting a region...")

# Mosaicking diferent AOI's (in this example parts 20 to 29 from May 11th):
gdal.BuildVRT('mosaic_2.vrt', sorted(glob.glob( 'C://VLAB//Samples//Flood-Composite//*20200511*part02*.tif'), key=os.path.getmtime))
#translate_options = gdal.TranslateOptions(projWin = [-85.0, 16.0, -64.0, 24.0])
translate_options = gdal.TranslateOptions(projWin = [-85.0, 24.0, -64.0, 16.0])
gdal.Translate('mosaic_2.tif',  'mosaic_2.vrt', options=translate_options)
#------------------------------------------------------------------------------------------------------
#------------------------------------------------------------------------------------------------------

# Finish the counter
print('Finished. Total processing time:', round((t.time() - start),2), 'seconds.')

A single Flood Map Area of Interest (AOI’s) converted to KMZ using Python

Areas of Interest (AOI’s) number 21 to 29 mosaicked and converted to KMZ using Python

After mosaicking different AOI’s, the GeoTIFF has been cut to a custom region using the Python script above

Please download the script by clicking at the image below.