The Advected Layer Precipitable Water product from CIRA will be added to the GNC-A broadcast in the near future.
The preliminary tests with the HDF files look very promising. Please find below some preliminary plots.
Received the following question:
I attended the Virtual Training on GEO Satellite Imagery Processing and I would like to combine the custom color scale from Script 08 with the reprojected imagery from Script 14. How can I do that?
Please find an example script below, with both concepts combined:
| # Training: Python and GOES-R Imagery: Script 14 - Reprojection with GDAL | |
| #----------------------------------------------------------------------------------------------------------- | |
| # 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 osr # Python bindings for GDAL | |
| from osgeo import gdal # Python bindings for GDAL | |
| import numpy as np # Scientific computing with Python | |
| from utilities import download_CMI # Our function for download | |
| 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) | |
| # Desired extent | |
| extent = [-64.0, -35.0, -35.0, -15.0] # Min lon, Max lon, Min lat, Max lat | |
| # AMAZON repository information | |
| # https://noaa-goes16.s3.amazonaws.com/index.html | |
| bucket_name = 'noaa-goes16' | |
| product_name = 'ABI-L2-CMIPF' | |
| yyyymmddhhmn = '202102181800' | |
| band = '13' | |
| # Download the file | |
| file_name = download_CMI(yyyymmddhhmn, band, input) | |
| #----------------------------------------------------------------------------------------------------------- | |
| # CUSTOM COLORMAP EXAMPLE 2 | |
| # Creating the INPE DISSM IR colormap | |
| # Online color picker: https://imagecolorpicker.com/ | |
| from matplotlib import cm # Colormap handling utilities | |
| vmin2 = -80 # Min. value | |
| vmax2 = 40 # Max. value | |
| gray_cmap = cm.get_cmap('gray_r', 120) # Read the reversed 'gray' cmap | |
| gray_cmap = gray_cmap(np.linspace(0, 1, 120)) # Create the array | |
| colors = ["#ffa0ff", "#0806ff", "#3bcfff", "#feff65", "#ff7516"] # Custom colors | |
| my_colors = cm.colors.ListedColormap(colors) # Create a custom colormap | |
| my_colors = my_colors(np.linspace(0, 1, 50)) # Create the array | |
| gray_cmap[:50, :] = my_colors # Join both cmaps arrays | |
| my_cmap2 = cm.colors.ListedColormap(gray_cmap) # Create the custom colormap | |
| #----------------------------------------------------------------------------------------------------------- | |
| # Variable | |
| var = 'CMI' | |
| # Open the file | |
| img = gdal.Open(f'NETCDF:{input}/{file_name}.nc:' + var) | |
| # 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) | |
| # Apply the scale, offset and convert to celsius | |
| ds = (ds * scale + offset) - 273.15 | |
| # Read the original file projection and configure the output projection | |
| source_prj = osr.SpatialReference() | |
| source_prj.ImportFromProj4(img.GetProjectionRef()) | |
| target_prj = osr.SpatialReference() | |
| target_prj.ImportFromProj4("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs") | |
| # Reproject the data | |
| GeoT = img.GetGeoTransform() | |
| driver = gdal.GetDriverByName('MEM') | |
| raw = driver.Create('raw', ds.shape[0], ds.shape[1], 1, gdal.GDT_Float32) | |
| raw.SetGeoTransform(GeoT) | |
| raw.GetRasterBand(1).WriteArray(ds) | |
| # Define the parameters of the output file | |
| kwargs = {'format': 'netCDF', \ | |
| 'srcSRS': source_prj, \ | |
| 'dstSRS': target_prj, \ | |
| 'outputBounds': (extent[0], extent[3], extent[2], extent[1]), \ | |
| 'outputBoundsSRS': target_prj, \ | |
| 'outputType': gdal.GDT_Float32, \ | |
| 'srcNodata': undef, \ | |
| 'dstNodata': 'nan', \ | |
| 'xRes': 0.02, \ | |
| 'yRes': 0.02, \ | |
| 'resampleAlg': gdal.GRA_NearestNeighbour} | |
| # Write the reprojected file on disk | |
| gdal.Warp(f'{output}/{file_name}_ret.nc', raw, **kwargs) | |
| #----------------------------------------------------------------------------------------------------------- | |
| # Open the reprojected GOES-R image | |
| file = Dataset(f'{output}/{file_name}_ret.nc') | |
| # 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]] | |
| # Define the color scale based on the channel | |
| colormap = "gray_r" # White to black for IR channels | |
| # Plot the image | |
| img = ax.imshow(data, origin='upper', extent=img_extent, vmin=vmin2, vmax=vmax2, cmap=my_cmap2) | |
| # Add coastlines, borders and gridlines | |
| ax.coastlines(resolution='10m', color='white', linewidth=0.8) | |
| ax.add_feature(cartopy.feature.BORDERS, edgecolor='white', 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 | |
| gl.xpadding = -100 | |
| gl.ypadding = -100 | |
| gl.xlabel_style = {'color': 'red'} | |
| gl.ylabel_style = {'color': 'red'} | |
| # 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 Band 13 ' + 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}/{file_name}_ret.png', bbox_inches='tight', pad_inches=0, dpi=300) | |
| # Show the image | |
| plt.show() |
Received the following question:
I would like to crop the GOES-16 DSIF files for an specific region (based on lat lon coordinates), but I do not know how to determine the x and y intervals. Do you have an example script for this?
Please find below an example script. The conversion between lat / lons and GOES-R x / y (seen on the utilities.py script) may be found on the GOES-R Series PUG, page 21.
| # Python and GOES-R Imagery Exampe: Cropping the DSIF imagery | |
| #----------------------------------------------------------------------------------------------------------- | |
| # 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 utilities import download_CMI # Our own utilities | |
| from utilities import geo2grid, convertExtent2GOESProjection # Our own utilities | |
| #----------------------------------------------------------------------------------------------------------- | |
| # Input and output directories | |
| input = "Samples"; os.makedirs(input, exist_ok=True) | |
| output = "Output"; os.makedirs(output, exist_ok=True) | |
| # Desired extent | |
| extent = [-64.0, -36.0, -40.0, -15.0] # Min lon, Max lon, Min lat, Max lat | |
| #----------------------------------------------------------------------------------------------------------- | |
| # Open the GOES-R image | |
| file = Dataset('OR_ABI-L2-DSIF-M6_G16_s20210632010164_e20210632019472_c20210632020547.nc') | |
| # Convert lat/lon to grid-coordinates | |
| lly, llx = geo2grid(extent[1], extent[0], file) | |
| ury, urx = geo2grid(extent[3], extent[2], file) | |
| # Get the pixel values | |
| data = file.variables['CAPE'][ury:lly, llx:urx] | |
| #----------------------------------------------------------------------------------------------------------- | |
| # Compute data-extent in GOES projection-coordinates | |
| img_extent = convertExtent2GOESProjection(extent) | |
| #----------------------------------------------------------------------------------------------------------- | |
| # 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.Geostationary(central_longitude=-75.0, satellite_height=35786023.0)) | |
| # Define the color scale based on the channel | |
| colormap = "jet" # White to black for IR channels | |
| # Plot the image | |
| img = ax.imshow(data, origin='upper', extent=img_extent, cmap=colormap) | |
| # 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) | |
| ax.gridlines(color='black', alpha=0.5, linestyle='--', linewidth=0.5) | |
| # Add a colorbar | |
| plt.colorbar(img, label='CAPE', extend='both', orientation='horizontal', pad=0.05, fraction=0.05) | |
| # Extract the date | |
| date = (datetime.strptime(file.time_coverage_start, '%Y-%m-%dT%H:%M:%S.%fZ')) | |
| # Add a title | |
| plt.title('GOES-16 CAPE ' + 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'OR_ABI-L2-DSIF-M6_G16_s20210632010164_e20210632019472_c20210632020547.png', bbox_inches='tight', pad_inches=0, dpi=300) | |
| # Show the image | |
| plt.show() |
You may find the utilities.py script, used in this example, at this link.
Received the following question:
I’m trying to process MODIS / TERRA “MOD14” data (Thermal Anomalies / Fires) downloaded from: https://ladsweb.modaps.eosdis.nasa.gov/search/
These images are in the HDF format. Do you have any examples to read it with Python and PyHDF?
Yes! The former Blended TPW files from GNC-A were HDF’s and I had to use PyHDF to read them. Please find an example code below, which worked with a MODIS sample downloaded from the provided link:
| from pyhdf.SD import SD, SDC # Import the HDF library | |
| # Reading the file: | |
| path = "MOD14.A2021085.0315.006.2021085095713.hdf" | |
| file = SD(path, SDC.READ) | |
| # Printing all the datasets names | |
| datasets_dic = file.datasets() | |
| for idx,sds in enumerate(datasets_dic.keys()): | |
| print (idx,sds) | |
| # Reading the dataset | |
| sds_obj = file.select('fire mask') # select sds | |
| data = sds_obj.get() # get sds data | |
| # Check the dimensions | |
| print(data.shape) |
Hint: If you want to quickly check the content of HDF, NetCDF or GRIB files without programming, I recommend using NASA’s Panoply:
Hi GEONETCasters!
Please find below the presentations from the 9th GEONETCast-Americas User Group Webinar (click on the images to download):
NOAA UPDATES
INPE UPDATES
GNC-A USER CASE STUDIES
MARN – MINISTRY OF ENVIRONMENT AND NATURAL RESOURCES – EL SALVADOR
Thank you to all the participants!
Please find the presentations from the other GEONETCast-Americas User Group Webinars, at the following links:
GEONETCasters;
Save the date! We hope to see you there.
We have scheduled the next GNC-A User Group Meeting Webex/Telecon to be held on March 25th, 2021 at 11:00 EDT / 15:00 UTC to 12:30 EDT / 16:30 UTC.
Webex Link:
Meeting Number: 198 124 3991
Meeting Password: gnc1
Join by Phone Info: +1-415-527-5035 US Toll
Access code: 198 124 3991
We ask that attendees use the Webex VOIP or the above number to join the webex.
The slides will be sent out shortly before the meeting.
Tentative Agenda:
1. Opening remarks/role call
2. NOAA updates
3. GNC-A Broadcast Status
4. INPE updates
5. User case study
6. User Outreach
Please download SHOWCast v.2.3.0 at the link below:
Changes on this version:
In this initial set, the following plots / scripts are provided: 2 m Temperature, Total Precipitation, PSML + Winds (850 hPa), Galvez Davison Index (GDI), Relative Vorticity + Geopotential Height (500 hPa), PSML + Geopotential Difference (1000 – 850 hPa) + Winds (925 hPa) + Mean RH (1000~400 hPa), Precipitable Water + CAPE, Low-Middle-High Clouds, Wind Speed and Direction (500 to 850 hPa).
By default, only the analysis is plotted, and users may change the configuration on the gfs scripts. In the example below, the GFS plots will be made from 00 h (analysis) until 90 h, with 3 h intervals. In order to plot only the analysis, just leave “hour_end” and “hour_inc” with the same value.
Note: In the showcast_config.py please select as “True” OR South America OR Central America + Caribbean.
Note: By default, for the density and heatmap products the GLM data is accumulated every 10 minutes. Users may change that if needed, in the “acum_minutes” variable:
In this initial set, at the showcast_config.py script, users may select to process the following ISCS data: SYNOP, DRIFTING BUOYS, METAR, SPECI, TAF, SIGMETS, AIRMETS, VOLCANIC ASH and TSUNAMI Warnings. At the “identifier” variable for each product, users may select the WMO Country Territory Designators: https://www.weather.gov/tg/tablec1
Please find below the procedure to install SHOWCast v 2:
SHOWCast 2 Installation Procedure
Please find below previous posts about SHOWCast:
There will be a GOES outage today at 1400Z that will last for 6 hours until 2000Z. This will affect all CMI and non-CMI full disk imagery for both GOES-16 and GOES-17, as well as the 20-second GLM data provided by NESDIS. The 5-minute aggregated INPE provided GLM product will still be available as well as all the other non-NOAA data products.
During the outage, users have the option to use the SHOWCast Cloud module, to download data automatically from AWS or UNIDATA THREDDS.
Stay tuned for news!
GNC-A Blog 