Converting Swath to Grid data using Python
Introduction
We are going to convert the swath data to grid data of Brightness Temperature (Tb) of Special Sensor Microwave Imager/Sounder (SSMIS) which is a satellite passive microwave radiometer aboard Defense Meteorological Satellite Program (DMSP) (F15, F16, F17 and F18) since 1987. These are near-polar orbiting satellites. SSMIS channels used to calculate brightness temperatures include 19.3 GHz vertical and horizontal, 22.2 GHz vertical, 37.0 GHz vertical and horizontal, and 91.7 GHz vertical and horizontal.
Details are available here. We are going to project the final gridded data to Polar Stereographic South.
Pre-requisite knowledge
Knowledge of numpy, matplotliband map projection in general are needed. And That’s it!!!
Data
The NetCDF file used in this excercise can be found here as “RSS_SSMIS_FCDR_V07R00_F17_D20200101_S0120_E0312_R67894.nc” file. Please refer to this article for more information about reading NetCDF file.
Necessary Libraries
NetCDF files can be read with a few different Python modules. We will use netCDF4. We will also usepyresample library. Pyresample is a python package for resampling geospatial image data. You can use pyresample to transform data in one coordinate system to another coordinate system (ex. mercator projection) using one of the available resampling algorithms. In this exercise, pyresample is used to convert swath data to grid data using resampling method. Details about this library and how to install it can be found here.
Reading netCDF file
We will use the netCDF.Dataset() attribute for reading the netCDF file. Let’s read the netCDF file and explore the variables:
#Importing the necessary libraries
from netCDF4 import Dataset
from pyresample import kd_tree, geometry, save_quicklook, SwathDefinition
#Reading the netCDF file
data = Dataset(r'\path\to\file.nc')
print (data)The output:
<class 'netCDF4._netCDF4.Dataset'>
root group (NETCDF4 data model, file format HDF5):
Conventions: CF-1.6
Metadata_Conventions: CF-1.6, Unidata Dataset Discovery v1.0, NOAA CDR v1.0, GDS v2.0
standard_name_vocabulary: CF Standard Name Table (v16, 11 October 2010)
id: RSS_SSMIS_FCDR_V07R00_F17_D20200101_S0120_E0312_R67894.nc
naming_authority: gov.noaa.ncdc
metadata_link: gov.noaa.ncdc:C00810
title: RSS Version-7 SSMIS FCDR
product_version: v07r00
date_issued: 2012-12-31
summary: Remote Sensing Systems (RSS) Version-7 Special Sensor Microwave Imager Sounder (SSMIS) Fundamental Climate Data Record (FCDR); intercalibrated and homogenized brightness temperature polar-orbiting product with quality flags
keywords: EARTH SCIENCE > SPECTRAL/ENGINEERING > MICROWAVE > BRIGHTNESS TEMPERATURE
keywords_vocabulary: NASA Global Change Master Directory (GCMD) Earth Science Keywords, Version 6.0
platform: DMSP 5D-2/F17 > Defense Meteorological Satellite Program-F17
sensor: SSMIS > Special Sensor Microwave Imager Sounder
cdm_data_type: Swath
cdr_program: Climate Data Record Program for satellites, FY 2011
cdr_variable: fcdr_brightness_temperature_19V, fcdr_brightness_temperature_19H, fcdr_brightness_temperature_22V, fcdr_brightness_temperature_37V, fcdr_brightness_temperature_37H, fcdr_brightness_temperature_92V, fcdr_brightness_temperature_92H
source: F17_r67894.dat
date_created: 20200106T210942Z
creator_name: Carl Mears
creator_url: http://www.remss.com/
creator_email: mears@remss.com
institution: Remote Sensing Systems
processing_level: NOAA Level 2
references: doi:10.1175/2007JAMC1635.1, doi:10.1126/science.1140746, SSM/I Users Interpretation Guide UG32268-900 Rev C 29 Nov 2000
history: 1) 2011-09-30, Carl Mears, Remote Sensing Systems, created initial netCDF file converted from the original RSS data format
geospatial_lat_min: -89.159996
geospatial_lat_max: 89.2
geospatial_lon_min: -180.0
geospatial_lon_max: 179.98999
geospatial_lat_units: degrees_north
geospatial_lon_units: degrees_east
spatial_resolution: 19 V/H GHz: 46.5km X 73.6km, 22 V GHz: 46.5km X 73.6km, 37 V/H GHz: 31.2km X 45.0km, 92 V/H GHz: 13.2km X 15.5km
time_coverage_start: 2020-01-01 01:20:02Z
time_coverage_end: 2020-01-01 03:12:05Z
time_coverage_duration: PT6723S
license: No restrictions on access or use
contributor_name: Frank Wentz, Carl Mears
contributor_role: Principal investigator and originator of input/source or antenna temperature data, Processor and author of entire driver routine (write_ssmi_tbs_netcdf.pro) which converts RSS native brightness temperature format to netCDF-4.
dimensions(sizes): scan_number(3542), footprint_number_hires(180), footprint_number_lores(90), eleven_flags(11), four_flags(4), numchar(20)
variables(dimensions): int32 iorbit(), float64 scan_time_hires(scan_number), float64 orbit_position(scan_number), float32 sc_lat(scan_number), float32 sc_lon(scan_number), float32 sc_alt(scan_number), int8 iscn_flag(eleven_flags, scan_number), int8 ical_flag_lores(four_flags, scan_number), int8 ical_flag_hires(four_flags, scan_number), int16 latitude_hires(scan_number, footprint_number_hires), int16 longitude_hires(scan_number, footprint_number_hires), int16 earth_incidence_angle_hires(scan_number, footprint_number_hires), int16 earth_azimuth_angle_hires(scan_number, footprint_number_hires), int16 sun_glitter_angle_hires(scan_number, footprint_number_hires), int8 land_flag_hires(scan_number, footprint_number_hires), int8 ice_flag_hires(scan_number, footprint_number_hires), float32 fcdr_brightness_temperature_92V(scan_number, footprint_number_hires), float32 fcdr_brightness_temperature_92H(scan_number, footprint_number_hires), int16 latitude_lores(scan_number, footprint_number_lores), int16 longitude_lores(scan_number, footprint_number_lores), int16 earth_incidence_angle_lores(scan_number, footprint_number_lores), int16 earth_azimuth_angle_lores(scan_number, footprint_number_lores), int16 sun_glitter_angle_lores(scan_number, footprint_number_lores), int8 land_flag_lores(scan_number, footprint_number_lores), int16 ice_flag_lores(scan_number, footprint_number_lores), float32 fcdr_brightness_temperature_19V(scan_number, footprint_number_lores), float32 fcdr_brightness_temperature_19H(scan_number, footprint_number_lores), float32 fcdr_brightness_temperature_22V(scan_number, footprint_number_lores), float32 fcdr_brightness_temperature_37V(scan_number, footprint_number_lores), float32 fcdr_brightness_temperature_37H(scan_number, footprint_number_lores)
groups: I know it is a lot of metadata. If you are new to netCDF file then I recommend you go through this article first.
Let’s read the variables:
#Reading the variables
longitudes = data ['longitude_lores'][:]
latitudes = data ['latitude_lores'][:]
tb = data['fcdr_brightness_temperature_19H'][:]The variables have the following shape:
Pyresample can be used to resample a swath dataset to a grid, a grid to a swath or a swath to another swath. Resampling can be done using nearest neighbour method, Guassian weighting, weighting with an arbitrary radial function.
Swath Definition
A swath is defined by the longitude and latitude values of the data points. We will use the SwathDefinition attribute.
#Swath Definition
swath_def = geometry.SwathDefinition(lons=longitudes, lats=latitudes)Area Definition
The cartographic definition of grid areas used by Pyresample is contained in an object of type AreaDefintion. The following arguments are needed to initialize an area:
- area_id ID of area
- name: Description
- proj_id: ID of projection
- proj_dict: Proj4 parameters as dict
- x_size: Number of grid columns
- y_size: Number of grid rows
- area_extent: (x_ll, y_ll, x_ur, y_ur)
where
- x_ll: projection x coordinate of lower left corner of lower left pixel
- y_ll: projection y coordinate of lower left corner of lower left pixel
- x_ur: projection x coordinate of upper right corner of upper right pixel
- y_ur: projection y coordinate of upper right corner of upper right pixel
Let’s implement this:
area_def = geometry.AreaDefinition('EPSG:3976 - WGS 84',
'NSIDC Sea Ice Polar Stereographic South',
'Stereographic',
'+proj=stere +lat_0=-90 +lat_ts=-70 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +type=crs',
316, 332,
[-3950000, -3950000,
3950000, 4350000])We have area_def : ‘EPSG:3976 — WGS 84’,
name: ‘NSIDC Sea Ice Polar Stereographic South’,
proj_id: ‘Stereographic’,
proj_dict: ‘+proj=stere +lat_0=-90 +lat_ts=-70 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +type=crs’,
x_size:316, y_size: 332,
area_extent: [-3950000, -3950000, 3950000, 4350000])
The above information can be found in NSIDC’s guide to Polar Stereographic Projection. The EPSG Code is 3976. The proj_dict is the PROJ.4 string of this particular EPSG. The gridded swath data will be 25 km resolution.
pyresample.kd_tree
This module contains several functions for resampling swath data. Let’s use resample_nearest attribute for resampling nearest neighbour method.
#Resampling using Nearest Neighbour Method
result = kd_tree.resample_nearest(swath_def, tb,area_def,
radius_of_influence=50000,
epsilon = 0.5)The keyword arguments:
- radius_of_influence: The radius around each grid pixel in meters to search for neighbours in the swath.
- epsilon: The distance to a found value is guaranteed to be no further than (1 + eps) times the distance to the correct neighbour. Allowing for uncertanty decreases execution time.
Displaying data quickly
Pyresample has some convenience functions for displaying data from a single channel. The function plot.save_quicklook saves the Cartopy image directly to file.
#Saving as a png file
save_quicklook('tb19h_quick_110m.png', area_def, result, label='19H Brightness Temperature (K)',
num_meridians=45,
num_parallels=10,
cmap='turbo', coast_res='10m')The keyword arguments:
- filename: path to output file. Here. “tb19h_quick_110m.png”.
- area_def: geometry.AreaDefinition object. Here, “area_def”.
- data: 2D array matching area_def. Use masked array for transparent values. Here, “result”.
- num_meridians: Number of meridians to plot on the globe.
- num_parallels: Number of parallels to plot on the globe.
- coast_res: Resolution of coastlines. “110m” or “50m” or “10m”.
The output:
Exporting to GeoTiff file
For those who wants to export as GeoTiff file, we need gdal library. For exporting files using gdal, please refer this article for more in depth information.
from osgeo import gdal, osr
#GeoReferencing information
gt = [-3950000.0, 25000.0, 0.0, 4350000.0, -0.0, -25000.0]
#Projection information
epsg = 3976
srs = osr.SpatialReference()
srs.ImportFromEPSG(epsg)
dest_wkt = srs.ExportToWkt()
#Creating output file
driver = gdal.GetDriverByName("GTiff")#Getting GeoTIFF driver
driver.Register()
outds = driver.Create ("SSMIS_Tb_19H.tif", #Name of the putput file
xsize = result.shape[1], #Setting the number of columns
ysize = result.shape[0], #Setting the number of rows
bands = 1, #Setting the number of bands
eType = gdal.GDT_Float32) #Setting the data type i.e., float 32
outds.SetGeoTransform(gt) #Setting the geoinformation
outds.SetProjection(dest_wkt) #Setting the projection information
outband = outds.GetRasterBand(1) #Setting the number of bands
outband.WriteArray(result) #Writing the 2-D array on the file
outband = None #Closing the file
outds = None #Closing the fileThe gt variable conatins georeferencing information in the form of the six rows.
The second row is the X-direction pixel size (here in meters). Third and fifth rows are the rotational components which in this case are set to zero as the image is not rotated. The sixth row is the Y-direction pixel resolution (which is almost always negative). The first and fourth rows are the coordinates (here in meters) of the upper left corner pixel for the image.
The dest_wkt variable contains the projection information which in our case is NSIDC Sea Ice Polar Stereographic South and the EPSG code is “3976”.
Let’s check if the output GeoTIFF file is properly created or not by opening the file using QGIS.
As we can see, the GeoTiff file is created. We add the OpenStreetMap and see that the output GeoTiff file is correctly georeferenced.
Conclusion
We have just learnt how to convert swath data to grid data using python. Also, we exported the gridded data as GeoTIFF file as it is much easier to deal with as it contains spatial referencing information in the form of georeferenced raster imagery.
