forked from diegormsouza/geo-sat-python-mar-2021
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Script_16.py
131 lines (104 loc) · 5.41 KB
/
Script_16.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
# Training: Python and GOES-R Imagery: Script 16 - Level 2 Products (SST) and Average
#-----------------------------------------------------------------------------------------------------------
# Required modules
from netCDF4 import Dataset # Read / Write NetCDF4 files
import matplotlib.pyplot as plt # Plotting library
from datetime import datetime # Basic Dates and time types
import cartopy, cartopy.crs as ccrs # Plot maps
import os # Miscellaneous operating system interfaces
from osgeo import gdal # Python bindings for GDAL
import numpy as np # Scientific computing with Python
from utilities import download_PROD # Our function for download
from utilities import reproject # Our function for reproject
gdal.PushErrorHandler('CPLQuietErrorHandler') # Ignore GDAL warnings
#-----------------------------------------------------------------------------------------------------------
# Input and output directories
input = "Samples"; os.makedirs(input, exist_ok=True)
output = "Output"; os.makedirs(output, exist_ok=True)
# AMAZON repository information
# https://noaa-goes16.s3.amazonaws.com/index.html
bucket_name = 'noaa-goes16'
product_name = 'ABI-L2-SSTF'
yyyymmdd = '20210218'
# Desired extent
extent = [-60.0, -40.0, -35.0, -20.0] # Min lon, Max lon, Min lat, Max lat
########################################################################
# Sea Surface Temperature - "X" Hours
########################################################################
sum_ds = np.zeros((5424,5424))
count_ds = np.zeros((5424,5424))
#-----------------------------------------------------------------------------------------------------------
for hour in np.arange(0,23,1):
# Date structure
yyyymmddhhmn = f'{yyyymmdd}{hour:02.0f}00'
# Download the file
file_name = download_PROD(yyyymmddhhmn, product_name, input)
#-----------------------------------------------------------------------------------------------------------
# Variable
var = 'SST'
# Open the GOES-R image
file = Dataset(f'{input}/{file_name}.nc')
# Open the file
img = gdal.Open(f'NETCDF:{input}/{file_name}.nc:' + var)
# Data Quality Flag (DQF)
dqf = gdal.Open(f'NETCDF:{input}/{file_name}.nc:DQF')
# Read the header metadata
metadata = img.GetMetadata()
scale = float(metadata.get(var + '#scale_factor'))
offset = float(metadata.get(var + '#add_offset'))
undef = float(metadata.get(var + '#_FillValue'))
dtime = metadata.get('NC_GLOBAL#time_coverage_start')
# Load the data
ds = img.ReadAsArray(0, 0, img.RasterXSize, img.RasterYSize).astype(float)
ds_dqf = dqf.ReadAsArray(0, 0, dqf.RasterXSize, dqf.RasterYSize).astype(float)
# Apply the scale, offset and convert to celsius
ds = (ds * scale + offset) - 273.15
# Apply NaN's where the quality flag is greater than 1
ds[ds_dqf > 1] = np.nan
# Calculate the sum
sum_ds = np.nansum(np.dstack((sum_ds,ds)),2)
nan_ds = np.isnan(ds)
count_ds = np.nansum(np.dstack((count_ds,(ds/ds))),2)
#-----------------------------------------------------------------------------------------------------------
# Calculate the sum
ds_day = np.empty((5424,5424))
ds_day[::] = np.nan
ds_day[count_ds!=0] = sum_ds[count_ds!=0]/count_ds[count_ds!=0]
#-----------------------------------------------------------------------------------------------------------
# Reproject the file
filename_ds = f'{output}/{file_name}_ret.nc'
reproject(filename_ds, img, ds_day, extent, undef)
#-----------------------------------------------------------------------------------------------------------
# Open the reprojected GOES-R image
file = Dataset(filename_ds)
# Get the pixel values
data = file.variables['Band1'][:]
#-----------------------------------------------------------------------------------------------------------
# Choose the plot size (width x height, in inches)
plt.figure(figsize=(10,10))
# Use the Geostationary projection in cartopy
ax = plt.axes(projection=ccrs.PlateCarree())
# Define the image extent
img_extent = [extent[0], extent[2], extent[1], extent[3]]
# Plot the image
img = ax.imshow(data, vmin=15, vmax=30, cmap='jet', origin='upper', extent=img_extent)
# Add coastlines, borders and gridlines
ax.coastlines(resolution='10m', color='black', linewidth=0.8)
ax.add_feature(cartopy.feature.BORDERS, edgecolor='black', linewidth=0.5)
gl = ax.gridlines(crs=ccrs.PlateCarree(), color='gray', alpha=1.0, linestyle='--', linewidth=0.25, xlocs=np.arange(-180, 180, 5), ylocs=np.arange(-90, 90, 5), draw_labels=True)
gl.top_labels = False
gl.right_labels = False
plt.xlim(extent[0], extent[2])
plt.ylim(extent[1], extent[3])
# Add a colorbar
plt.colorbar(img, label='Brightness Temperatures (°C)', extend='both', orientation='horizontal', pad=0.05, fraction=0.05)
# Extract the date
date = (datetime.strptime(dtime, '%Y-%m-%dT%H:%M:%S.%fZ'))
# Add a title
plt.title('GOES-16 SST ' + date.strftime('%Y-%m-%d %H:%M') + ' UTC', fontweight='bold', fontsize=10, loc='left')
plt.title('Reg.: ' + str(extent) , fontsize=10, loc='right')
#-----------------------------------------------------------------------------------------------------------
# Save the image
plt.savefig(f'{output}/SST_DAY_RET_{yyyymmdd}.png', bbox_inches='tight', pad_inches=0, dpi=300)
# Show the image
plt.show()