Hi community,

An additional nice feature to the script from the first part:

• Insert city names based on their lat / longs

This is achieved using the following code:

print("Inserting City Names...")
# Cities latitutes
city_lons = [-46.6252, -43.1971, -49.2733, -48.5482, -51.2287, -54.6478, -43.9542, -40.2958, -45.0093]
# Cities latitutes
city_lats = [-23.5337, -22.9083, -25.4284, -27.5950, -30.0277, -20.4435, -19.8157, -20.2976, -22.6650]
x,y = bmap(city_lons, city_lats)
# Plot the city markers (circles)
bmap.plot(x, y, 'bo', marker = '.', markersize=10, color = 'gold', zorder=8)
# Choose the city labels
labels = ['São Paulo', 'Rio de Janeiro', 'Curitiba', 'Florianópolis', 'Porto Alegre', 'Campo Grande', 'Belo Horizonte', 'Vitória', 'Cachoeira\n Paulista']
# Offset the labels so they do not overlap the markers
x_offsets = [10000,10000,10000,10000,10000,10000,10000,10000,0]
y_offsets = [0,0,0,0,0,0,0,0,10000]
# Plot hte labels
for label, xpt, ypt, x_offset, y_offset in zip(labels, x, y, x_offsets, y_offsets):
plt.text(xpt+x_offset, ypt+y_offset, label, fontsize=12, fontweight='bold', zorder=8, color='white')

Please find below the full example script to generate these plots:

############################################################
# LICENSE
# Copyright (C) 2019 - INPE - NATIONAL INSTITUTE FOR SPACE RESEARCH
# 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 Libraries
import matplotlib.pyplot as plt
from netCDF4 import Dataset
import numpy as np
from datetime import datetime, timedelta
import time as t
from mpl_toolkits.basemap import Basemap
from pyproj import Proj
from pyorbital import astronomy
from cpt_convert import loadCPT
from matplotlib.colors import LinearSegmentedColormap
from osgeo import gdal
from PIL import Image

# Start the time counter to measure the total processing time
start = t.time()

# Load the Channel 02 contrast curve
contrast_curve = [0.00000, 0.02576, 0.05148, 0.07712, 0.10264, 0.12799, 0.15313, 0.17803, 0.20264, 0.22692, 0.25083, 0.27432, 0.29737, 0.31991, 0.34193, 0.36336,
0.38418, 0.40433, 0.42379, 0.44250, 0.46043, 0.47754, 0.49378, 0.50911, 0.52350, 0.53690, 0.54926, 0.56055, 0.57073, 0.57976, 0.58984, 0.59659,
0.60321, 0.60969, 0.61604, 0.62226, 0.62835, 0.63432, 0.64016, 0.64588, 0.65147, 0.65694, 0.66230, 0.66754, 0.67267, 0.67768, 0.68258, 0.68738,
0.69206, 0.69664, 0.70112, 0.70549, 0.70976, 0.71394, 0.71802, 0.72200, 0.72589, 0.72968, 0.73339, 0.73701, 0.74055, 0.74399, 0.74736, 0.75065,
0.75385, 0.75698, 0.76003, 0.76301, 0.76592, 0.76875, 0.77152, 0.77422, 0.77686, 0.77943, 0.78194, 0.78439, 0.78679, 0.78912, 0.79140, 0.79363,
0.79581, 0.79794, 0.80002, 0.80206, 0.80405, 0.80600, 0.80791, 0.80978, 0.81162, 0.81342, 0.81518, 0.81692, 0.81862, 0.82030, 0.82195, 0.82358,
0.82518, 0.82676, 0.82833, 0.82987, 0.83140, 0.83292, 0.83442, 0.83592, 0.83740, 0.83888, 0.84036, 0.84183, 0.84329, 0.84476, 0.84623, 0.84771,
0.84919, 0.85068, 0.85217, 0.85368, 0.85520, 0.85674, 0.85829, 0.85986, 0.86145, 0.86306, 0.86469, 0.86635, 0.86803, 0.86974, 0.87149, 0.87326,
0.87500, 0.87681, 0.87861, 0.88038, 0.88214, 0.88388, 0.88560, 0.88730, 0.88898, 0.89064, 0.89228, 0.89391, 0.89552, 0.89711, 0.89868, 0.90023,
0.90177, 0.90329, 0.90479, 0.90627, 0.90774, 0.90919, 0.91063, 0.91205, 0.91345, 0.91483, 0.91620, 0.91756, 0.91890, 0.92022, 0.92153, 0.92282,
0.92410, 0.92536, 0.92661, 0.92784, 0.92906, 0.93027, 0.93146, 0.93263, 0.93380, 0.93495, 0.93608, 0.93720, 0.93831, 0.93941, 0.94050, 0.94157,
0.94263, 0.94367, 0.94471, 0.94573, 0.94674, 0.94774, 0.94872, 0.94970, 0.95066, 0.95162, 0.95256, 0.95349, 0.95441, 0.95532, 0.95622, 0.95711,
0.95799, 0.95886, 0.95973, 0.96058, 0.96142, 0.96225, 0.96307, 0.96389, 0.96469, 0.96549, 0.96628, 0.96706, 0.96783, 0.96860, 0.96936, 0.97010,
0.97085, 0.97158, 0.97231, 0.97303, 0.97374, 0.97445, 0.97515, 0.97584, 0.97653, 0.97721, 0.97789, 0.97856, 0.97922, 0.97988, 0.98053, 0.98118,
0.98182, 0.98246, 0.98309, 0.98372, 0.98435, 0.98497, 0.98559, 0.98620, 0.98681, 0.98741, 0.98802, 0.98862, 0.98921, 0.98980, 0.99039, 0.99098,
0.99157, 0.99215, 0.99273, 0.99331, 0.99389, 0.99446, 0.99503, 0.99561, 0.99618, 0.99675, 0.99732, 0.99788, 0.99845, 0.99902, 0.99959, 1.000000]

# Convert the contrast curve to a numpy array
curve = np.array(contrast_curve)

# Visualization extent [min_lon, min_lat, max_lon, max_lat]
extent = [-55, -28.0, -40.0, -18.0]

print("Converting G16 coordinates to lons and lats...")

# Band 02 file to read (Note: using 1 km Band 02 from GNC-A. If using 500 m, you should read every 2 pixels)
image1 = "FalseColor\OR_ABI-L2-CMIPF-M6C13_G16_s20190931000228_e20190931009547_c20190931010017.nc"
# Read the Band 02 file using the NetCDF library
file1 = Dataset(image1)
# Satellite height
sat_h = file1.variables['goes_imager_projection'].perspective_point_height
# Satellite longitude
sat_lon = file1.variables['goes_imager_projection'].longitude_of_projection_origin
# Satellite sweep
sat_sweep = file1.variables['goes_imager_projection'].sweep_angle_axis
# The projection x and y coordinates equals
# the scanning angle (in radians) multiplied by the satellite height (http://proj4.org/projections/geos.html)
X = file1.variables['x'][:][::1] * sat_h
Y = file1.variables['y'][:][::1] * sat_h
# map object with pyproj
p = Proj(proj='geos', h=sat_h, lon_0=sat_lon, sweep=sat_sweep, a=6378137.0)
# Convert map points to latitude and longitude with the magic provided by Pyproj
XX, YY = np.meshgrid(X, Y)
lons, lats = p(XX, YY, inverse=True)

print("Converting min lon to G16 coordinates...")
value = extent[0]
X = np.abs( lons - value)
idx_min_lon = np.where( X == X.min() )
min_lon = np.round(idx_min_lon[1]/2.)*2
#print(idx_min_lon[1])
#print(min_lon)

print("Converting max lon to G16 coordinates...")
value = value = extent[2]
X = np.abs( lons - value)
idx_max_lon = np.where( X == X.min() )
max_lon = np.round(idx_max_lon[1]/2.)*2
#print(idx_max_lon[1])
#print(max_lon)

print("Converting min lat to G16 coordinates...")
value = extent[3]
X = np.abs( lats - value)
idx_min_lat = np.where( X == X.min() )
min_lat = np.round(idx_min_lat[0]/2.)*2
#print(idx_min_lat[0])
#print(min_lat)

print("Converting max lat to G16 coordinates...")
value = extent[1]
X = np.abs( lats - value)
idx_max_lat = np.where( X == X.min() )
max_lat = np.round(idx_max_lat[0]/2.)*2

#print(idx_max_lat[0])
#print(max_lat)

XX = XX[int(min_lat):int(max_lat),int(min_lon):int(max_lon)][::1,::1]
YY = YY[int(min_lat):int(max_lat),int(min_lon):int(max_lon)][::1,::1]

#print(XX[0,0])
#print(XX[-1,-1])
#print(YY[0,0])
#print(YY[-1,-1])

print("Reading Band 02...")

# Band 02 file to read (Note: using 1 km Band 02 from GNC-A. If using 500 m, you should read every 2 pixels)
image1 = "FalseColor\OR_ABI-L2-CMIPF-M6C02_G16_s20190931000228_e20190931009536_c20190931010012-114300_0.nc"
# Read the Band 02 file using the NetCDF library
file1 = Dataset(image1)
# Extract the reflactance factor, every two pixels
data1 = file1.variables['CMI'][int(min_lat):int(max_lat),int(min_lon):int(max_lon)][::1,::1]
# Minimum and maximum reflectances
VISmin = 0.0
VISmax = 1.0
# Anything that is below the min or greater than max, keep min and max
data1[data1VISmax] = VISmax
# Convert to 0 - 255
data1 = ((data1 - VISmin) / (VISmax - VISmin)) * 255
# Convert to int
data1 = data1.astype(int)
# Apply the contrast curve
data1 = curve[data1] * 255
# Convert to int
data1 = data1.astype(int)
#print(data1)
#print(np.shape(data1))

print("Reading Band 13...")

# Band 13 file to read
image2 = "FalseColor\OR_ABI-L2-CMIPF-M6C13_G16_s20190931000228_e20190931009547_c20190931010017.nc"
# Read the Band 13 file using the NetCDF library
file2 = Dataset(image2)
# Extract the brightness temperatures, every pixel
data2 = file2.variables['CMI'][int(min_lat/2):int(max_lat/2),int(min_lon/2):int(max_lon/2)][::1,::1]
# Double the size of the Band 13 array, so it may be used with an 1 km image
data3 = np.repeat(data2, 2, axis=0)
data3 = np.repeat(data3, 2, axis=1)
# Convert Kelvin to Celsius (will be used later in the code)
data2a = data3
data2b = data3 - 273.15
# Minimum and maximum brightness temperatures
IRmin = 89.62
IRmax = 341.27
# Anything that is below the min or greater than max, keep min and max
data3[data3IRmax] = IRmax
# Convert to 0 - 255
data3 = ((data3 - IRmin) / (IRmax - IRmin)) * 255
# Convert to int
data2 = data3.astype(int)
#print(np.shape(data2))

print("Reading the Time and Date...")
# Getting the file date
add_seconds = int(file2.variables['time_bounds'][0])
date = datetime(2000,1,1,12) + timedelta(seconds=add_seconds)
date_format = date.strftime('%Y%m%d%H%M')
year = date.strftime('%Y')
month = date.strftime('%m')
day = date.strftime('%d')
hour = date.strftime('%H')
minutes = date.strftime('%M')

print("Reading the False Color Look Up Table...")

import matplotlib.image as mpimg
# Open the False Color Look Up Table
img = mpimg.imread('wx-star.com_GOES-R_ABI_False-Color-LUT_sat20.png')
# Flip the array (for some reason, is necessary to flip the LUT horizontally)
img = np.fliplr(img)
# Apply the look up table based on the Band 02 reflectances and Band 13 Brightness Temperatures
final_image = img[data1,data2,0:3]

print("Calculating the sun zenith angle...")
utc_time = datetime(int(year), int(month), int(day), int(hour), int(minutes))
sun_zenith = np.zeros((data1.shape[0], data1.shape[1]))
sun_zenith = astronomy.sun_zenith_angle(utc_time, lons[int(min_lat):int(max_lat),int(min_lon):int(max_lon)][::1,::1], lats[int(min_lat):int(max_lat),int(min_lon):int(max_lon)][::1,::1])
#print(np.shape(sun_zenith))

print("Putting zeroes in the areas without sunlight")
final_image[sun_zenith > 80] = 0
mask = (final_image == [0.0,0.0,0.0]).all(axis=2)
#apply the mask to overwrite the pixels
final_image[mask] = [0,0,0]
alpha = ~np.all(final_image == 0, axis=2)
final_image_alpha = np.dstack((final_image, alpha))
#print(np.shape(final_image))
#print(np.shape(final_image_alpha))

print("Reading the night lights...")

raster = gdal.Open('BlackMarble_2016_6km_geo.tif')
ulx, xres, xskew, uly, yskew, yres = raster.GetGeoTransform()
lrx = ulx + (raster.RasterXSize * xres)
lry = uly + (raster.RasterYSize * yres)
corners = [ulx, lry, lrx, uly]
min_lon = extent[0]; max_lon = extent[2]; min_lat = extent[1]; max_lat = extent[3]
raster = gdal.Translate('teste.tif', raster, projWin = [min_lon, max_lat, max_lon, min_lat])
array2 = raster.ReadAsArray()
r = array2[0,:,:].astype(float)
g = array2[1,:,:].astype(float)
b = array2[2,:,:].astype(float)
r[r==4] = 0
g[g==5] = 0
b[b==15] = 0
geo = raster.GetGeoTransform()
xres = geo[1]
yres = geo[5]
xmin = geo[0]
xmax = geo[0] + (xres * raster.RasterXSize)
ymin = geo[3] + (yres * raster.RasterYSize)
ymax = geo[3]
lons_n = np.arange(xmin,xmax,xres)
lats_n = np.arange(ymax,ymin,yres)
lons_n, lats_n = np.meshgrid(lons_n,lats_n)
color_tuples = (np.array([r[:-1,:-1].flatten(), g[:-1,:-1].flatten(), b[:-1,:-1].flatten()]).transpose())/255
#print(np.shape(color_tuples))

print("Creating the plot...")

# Define the size of the saved picture=========================================
DPI = 150
fig = plt.figure(figsize=(data1.shape[1]/float(DPI), data1.shape[0]/float(DPI)), frameon=False, dpi=DPI)
ax = plt.Axes(fig, [0., 0., 1., 1.])
ax.set_axis_off()
fig.add_axes(ax)
ax = plt.axis('off')
#==============================================================================

# Create the basemap for geostationary view
bmap = Basemap(projection='geos', rsphere=(6378140.0,6356750.0), resolution='c', area_thresh=0, lon_0=-75.0, satellite_height=35786000, llcrnrx= XX[0,0], llcrnry= YY[-1,-1] , urcrnrx= XX[-1,-1], urcrnry= YY[0,0])

print("First layer...")
# FIRST LAYER
# Converts a CPT file to be used in Python
cpt = loadCPT('IR4AVHRR6.cpt')
IR4AVHRR6 = LinearSegmentedColormap('cpt', cpt)
# Everything under 88° is a nan
data2a[sun_zenith < 80] = np.nan
# Apply range limits for clean IR channel
data2a = np.maximum(data2a, 90)
data2a = np.minimum(data2a, 313)
# Normalize the channel between a range
data2a = (data2a-90)/(313-90)
# Invert colors
data2a = 1 - data2a
bmap.imshow(data2a, cmap='gray', vmin=0.1, vmax=0.25, alpha = 0.4, origin='upper', zorder=4)

print("Second layer...")
# SECOND LAYER
data2a[data2a < 0.20] = np.nan
bmap.imshow(data2a, cmap='gray', vmin=0.15, vmax=0.30, alpha = 1.0, origin='upper', zorder=6)

print("Third layer...")
# THIRD LAYER
data2b[sun_zenith < 80] = np.nan
data2b[np.logical_or(data2b  -28)] = np.nan
bmap.imshow(data2b, cmap=IR4AVHRR6, vmin=-103, vmax=84, alpha=1.0, origin='upper', zorder=7)

print("Fourth layer...")
# FOURTH LAYER
# Show th night lights
bmap.pcolormesh(lons_n, lats_n, r, latlon=True, color=color_tuples, zorder=1)

print("Fifth layer...")
# FIFTH LAYER
# Show the image
bmap.imshow(final_image_alpha, origin='upper', zorder=2)

print("Configuring the map...")
# Insert shapefiles
bmap.readshapefile('BRMUE250GC_SIR', 'BRMUE250GC_SIR', linewidth=0.1, color='white', zorder=6)
bmap.readshapefile('estados_2010', 'estados_2010', linewidth=1.0, color='yellow', zorder=6)
bmap.readshapefile('Americas', 'Americas', linewidth=2.0, color='orange', zorder=6)
bmap.drawparallels(np.arange(-90.0, 90.0, 5), dashes = [4 ,4], linewidth=0.5, color='white', zorder=7)
bmap.drawmeridians(np.arange(0.0, 360.0, 5), dashes = [4,4], linewidth=0.5, color='white', zorder=7)

print("Inserting Title...")
plt.annotate('GOES-16 RGB ' + date.strftime('%B %d %Y - %H:%M UTC (INPE/CPTEC)'), xy=(20, 20), xycoords='figure pixels', fontsize=13, fontweight='bold', zorder=8, color='whitesmoke')

print("Inserting City Names...")
# Cities latitutes
city_lons = [-46.6252, -43.1971, -49.2733, -48.5482, -51.2287, -54.6478, -43.9542, -40.2958, -45.0093]
# Cities latitutes
city_lats = [-23.5337, -22.9083, -25.4284, -27.5950, -30.0277, -20.4435, -19.8157, -20.2976, -22.6650]
x,y = bmap(city_lons, city_lats)
# Plot the city markers (circles)
bmap.plot(x, y, 'bo', marker = '.', markersize=10, color = 'gold', zorder=8)
# Choose the city labels
labels = ['São Paulo', 'Rio de Janeiro', 'Curitiba', 'Florianópolis', 'Porto Alegre', 'Campo Grande', 'Belo Horizonte', 'Vitória', 'Cachoeira\n Paulista']
# Offset the labels so they do not overlap the markers

x_offsets = [10000,10000,10000,10000,10000,10000,10000,10000,0]
y_offsets = [0,0,0,0,0,0,0,0,10000]
# Plot hte labels
for label, xpt, ypt, x_offset, y_offset in zip(labels, x, y, x_offsets, y_offsets):
plt.text(xpt+x_offset, ypt+y_offset, label, fontsize=12, fontweight='bold', zorder=8, color='white')

print("Saving Image...")
# Save the image
plt.savefig('G16_GNCOLOR_' + date_format + '.png', dpi=DPI, facecolor="black", pad_inches=0)

print ('- finished! Time:', t.time() - start, 'seconds')

Please download at this link, the additional files required to tun this example script (shapefiles, night lights GeoTIFF, color tables). You should extract the zip in the script’s directory. Stay tuned for news.