This module shows the structure of VIIRS Level 1B
data and what information of the data files can be used in order to load, browse and visualize the data.
According to NASA, "VIIRS is aboard two satellites: Suomi National Polar-orbiting Partnership (Suomi NPP) and National Oceanic and Atmospheric Administration-20 (NOAA-20).
VIIRS observes the entire Earth’s surface twice each day crossing equator approximately 14 times each day per satellite. Suomi NPP crosses the equator about 1:30 a.m. and 1:30 p.m. and NOAA-20 crossing at 2:20 am and 2:20 pm, local time; from NPP’s polar orbit 824 kilometers (km) (512 miles) above the Earth’s surface. [...] The VIIRS instrument provides 22 spectral bands from 412 nanometers (nm) to 12 micrometers (µm) at two spatial resolutions, 375 meters (m) and 750 m, which are resampled to 500 m, 1 km, and 0.05 degrees in the NASA produced data products to promote consistency with the MODIS heritage."
Spatial resolution:
375m or 750m at nadir
Spatial coverage:Global
Revisit time:Daily
Data availability:since 2012
This notebook uses VIIRS Level 1B data from the SNPP platform. This data can be ordered via the LAADS DAAC and are distributed in netcdf
format.
Note that you need to also search for and download the 6 data files in order to get all the data you need to geolocate the data and create the full set of composites:
The filenames for the same granule will have a common middle portion which reflects the date and time. For example, A2021037.1306
. This means the file date and time is 2021-02-06 13:06:00
.
You need to register for an Earthdata account in order to be able to download data. Read this webpage for more information on downloading data programmatically using wget
.
import glob
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors
from matplotlib.axes import Axes
import satpy
from satpy.scene import Scene
from satpy import find_files_and_readers
import pyresample as prs
import cartopy.crs as ccrs
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
from cartopy.mpl.geoaxes import GeoAxes
GeoAxes._pcolormesh_patched = Axes.pcolormesh
import warnings
warnings.filterwarnings('ignore')
warnings.simplefilter(action = "ignore", category = RuntimeWarning)
%run ../functions.ipynb
From the LAADS DAAC, we downloaded one tile of Level-1B Image data on 6 February 2021. The data is available in the folder ../../eodata/1_satellite/viirs/level_1b/2021/02/06/
. Let us load the data. First, we specify the file path and create a variable with the name file_name
.
file_names = glob.glob('../../../eodata/dust/part1/1_satellite/viirs/level_1b/2021/02/06/*A2021037*')
file_names
['../../../eodata/dust/part1/1_satellite/viirs/level_1b/2021/02/06/VNP02IMG.A2021037.1306.002.2021128093451.nc', '../../../eodata/dust/part1/1_satellite/viirs/level_1b/2021/02/06/VNP02DNB.A2021037.1306.002.2021128093451.nc', '../../../eodata/dust/part1/1_satellite/viirs/level_1b/2021/02/06/VNP02MOD.A2021037.1306.002.2021128093451.nc', '../../../eodata/dust/part1/1_satellite/viirs/level_1b/2021/02/06/VNP03IMG.A2021037.1306.002.2021126191945.nc', '../../../eodata/dust/part1/1_satellite/viirs/level_1b/2021/02/06/VNP03MOD.A2021037.1306.002.2021126191945.nc', '../../../eodata/dust/part1/1_satellite/viirs/level_1b/2021/02/06/VNP03DNB.A2021037.1306.002.2021126191945.nc']
In a next step, we use the Scene
constructor from the satpy library. Once loaded, a Scene
object represents a single geographic region of data, typically at a single continuous time range.
You have to specify the two keyword arguments reader
and filenames
in order to successfully load a scene. As mentioned above, for VIIRS Level-1B data, you can use the viirs_l1b
reader.
scn =Scene(filenames=file_names,reader='viirs_l1b')
scn
<satpy.scene.Scene at 0x7f1594590fa0>
A Scene
object is a collection of different bands, with the function available_dataset_names()
, you can see the available bands of the scene. To learn more about the bands of VIIRS, visit this website.
scn.available_dataset_names()
['DNB', 'I01', 'I02', 'I03', 'I04', 'I05', 'M01', 'M02', 'M03', 'M04', 'M05', 'M06', 'M07', 'M08', 'M09', 'M10', 'M11', 'M12', 'M13', 'M14', 'M15', 'M16', 'dnb_lat', 'dnb_lon', 'dnb_lunar_azimuth_angle', 'dnb_lunar_zenith_angle', 'dnb_moon_illumination_fraction', 'dnb_satellite_azimuth_angle', 'dnb_satellite_zenith_angle', 'dnb_solar_azimuth_angle', 'dnb_solar_zenith_angle', 'i_lat', 'i_lon', 'i_satellite_azimuth_angle', 'i_satellite_zenith_angle', 'i_solar_azimuth_angle', 'i_solar_zenith_angle', 'm_lat', 'm_lon', 'satellite_azimuth_angle', 'satellite_zenith_angle', 'solar_azimuth_angle', 'solar_zenith_angle']
The underlying container for data in satpy is the xarray.DataArray
. With the function load()
, you can specify an individual band by name, e.g. I01
and load the data. If you then select the loaded band, you see that the band object is a xarray.DataArray
. Band I01 has a bandwidth or wavelength of 0.6 - 0.68μm.
scn.load(['I01'])
scn['I01']
<xarray.DataArray 'I01' (y: 6496, x: 6400)> dask.array<add, shape=(6496, 6400), dtype=float64, chunksize=(4096, 4096), chunktype=numpy.ndarray> Coordinates: crs object +proj=latlong +datum=WGS84 +ellps=WGS84 +type=crs Dimensions without coordinates: y, x Attributes: start_orbit: 48084 end_orbit: 48084 long_name: I-band 01 observations at pixel locations _FillValue: 65535 valid_min: 0 valid_max: 65527 scale_factor: 1.9991758e-05 add_offset: 0.0 radiance_scale_factor: 0.010676029 radiance_add_offset: 0.0 radiance_units: Watts/meter^2/steradian/micrometer name: I01 wavelength: 0.64 µm (0.6-0.68 µm) resolution: 371 calibration: reflectance coordinates: ('i_lon', 'i_lat') file_type: vl1bi standard_name: toa_bidirectional_reflectance units: % modifiers: () rows_per_scan: 32 shape: (6496, 6400) file_units: 1 platform_name: Suomi-NPP sensor: VIIRS start_time: 2021-02-06 13:06:00 end_time: 2021-02-06 13:12:00 reader: viirs_l1b area: Shape: (6496, 6400)\nLons: <xarray.DataArray 'lon... _satpy_id: DataID(name='I01', wavelength=WavelengthRange(min... ancillary_variables: []
|
array(<Geographic 2D CRS: +proj=latlong +datum=WGS84 +ellps=WGS84 +type=crs> Name: unknown Axis Info [ellipsoidal]: - lon[east]: Longitude (degree) - lat[north]: Latitude (degree) Area of Use: - undefined Datum: World Geodetic System 1984 - Ellipsoid: WGS 84 - Prime Meridian: Greenwich , dtype=object)
With an xarray data structure, you can handle the object as a xarray.DataArray
. For example, you can print a list of available attributes with the function attrs.keys()
.
scn['I01'].attrs.keys()
dict_keys(['start_orbit', 'end_orbit', 'long_name', '_FillValue', 'valid_min', 'valid_max', 'scale_factor', 'add_offset', 'radiance_scale_factor', 'radiance_add_offset', 'radiance_units', 'name', 'wavelength', 'resolution', 'calibration', 'coordinates', 'file_type', 'standard_name', 'units', 'modifiers', 'rows_per_scan', 'shape', 'file_units', 'platform_name', 'sensor', 'start_time', 'end_time', 'reader', 'area', '_satpy_id', 'ancillary_variables'])
With the attrs()
function, you can also access individual metadata information, e.g. start_time
and end_time
.
scn['I01'].start_time, scn['I01'].end_time
(datetime.datetime(2021, 2, 6, 13, 6), datetime.datetime(2021, 2, 6, 13, 12))
RGB composites combine different channels of satellite data in order to get e.g. a true-color image of the scene. Depending on which channel combination is used, different features can be highlighted in the composite, e.g. dust. SatPy offers several predefined RGB composite options. The function available_composite_ids()
returns a list of available composite IDs.
scn.available_composite_ids()
[DataID(name='adaptive_dnb'), DataID(name='ash'), DataID(name='cloudtop'), DataID(name='cloudtop_daytime'), DataID(name='day_microphysics'), DataID(name='dust'), DataID(name='dynamic_dnb'), DataID(name='false_color'), DataID(name='fire_temperature'), DataID(name='fire_temperature_39refl'), DataID(name='fire_temperature_awips'), DataID(name='fire_temperature_eumetsat'), DataID(name='fog'), DataID(name='green_snow'), DataID(name='histogram_dnb'), DataID(name='hncc_dnb'), DataID(name='hr_cloudtop_daytime'), DataID(name='hr_overview'), DataID(name='ir108_3d'), DataID(name='ir_cloud_day'), DataID(name='natural_color'), DataID(name='natural_color_raw'), DataID(name='natural_color_sun_lowres'), DataID(name='natural_enh'), DataID(name='natural_with_night_fog'), DataID(name='night_fog'), DataID(name='night_microphysics'), DataID(name='night_overview'), DataID(name='ocean_color'), DataID(name='overview'), DataID(name='overview_raw'), DataID(name='snow_age'), DataID(name='snow_lowres'), DataID(name='ssec_fog'), DataID(name='true_color'), DataID(name='true_color_crefl'), DataID(name='true_color_lowres'), DataID(name='true_color_lowres_crefl'), DataID(name='true_color_lowres_land'), DataID(name='true_color_lowres_marine_tropical'), DataID(name='true_color_raw')]
Let us define a list with the composite ID natural_color
and dust
. This list (composite_id
) can then be passed to the function load()
. Per default, scenes are loaded with the north pole facing downwards. You can specify the keyword argument upper_right_corner='NE'
in order to turn the image around and have the north pole facing upwards.
composite_ids = ['natural_color','dust']
scn.load(composite_ids, upper_right_corner='NE')
Inconsistent sensor/satellite input - sensor set to viirs
A print of the Scene object scn
shows you that three bands / composites are available: natural_color
, dust
and band I01
.
print(scn)
<xarray.DataArray 'I01' (y: 6496, x: 6400)> dask.array<add, shape=(6496, 6400), dtype=float64, chunksize=(4096, 4096), chunktype=numpy.ndarray> Coordinates: crs object +proj=latlong +datum=WGS84 +ellps=WGS84 +type=crs Dimensions without coordinates: y, x Attributes: start_orbit: 48084 end_orbit: 48084 long_name: I-band 01 observations at pixel locations _FillValue: 65535 valid_min: 0 valid_max: 65527 scale_factor: 1.9991758e-05 add_offset: 0.0 radiance_scale_factor: 0.010676029 radiance_add_offset: 0.0 radiance_units: Watts/meter^2/steradian/micrometer name: I01 wavelength: 0.64 µm (0.6-0.68 µm) resolution: 371 calibration: reflectance coordinates: ('i_lon', 'i_lat') file_type: vl1bi standard_name: toa_bidirectional_reflectance units: % modifiers: () rows_per_scan: 32 shape: (6496, 6400) file_units: 1 platform_name: Suomi-NPP sensor: VIIRS start_time: 2021-02-06 13:06:00 end_time: 2021-02-06 13:12:00 reader: viirs_l1b area: Shape: (6496, 6400)\nLons: <xarray.DataArray 'lon... _satpy_id: DataID(name='I01', wavelength=WavelengthRange(min... ancillary_variables: [] <xarray.DataArray 'concatenate-46f724bb7a6b05252ce1225392d80d11' (bands: 3, y: 3248, x: 3200)> dask.array<concatenate, shape=(3, 3248, 3200), dtype=float64, chunksize=(1, 3248, 3200), chunktype=numpy.ndarray> Coordinates: crs object +proj=latlong +datum=WGS84 +ellps=WGS84 +type=crs * bands (bands) <U1 'R' 'G' 'B' Dimensions without coordinates: y, x Attributes: start_orbit: 48084 add_offset: 0.0 start_time: 2021-02-06 13:06:00 sensor: VIIRS end_time: 2021-02-06 13:12:00 scale_factor: 1.9991758e-05 standard_name: natural_color resolution: 742 end_orbit: 48084 valid_min: 0 _FillValue: 65535 rows_per_scan: 16 file_units: 1 platform_name: Suomi-NPP valid_max: 65527 reader: viirs_l1b file_type: vl1bm area: Shape: (3248, 3200)\nLons: <xarray.DataArray 'lo... radiance_units: Watts/meter^2/steradian/micrometer radiance_add_offset: 0.0 shape: (3248, 3200) coordinates: ('m_lon', 'm_lat') ancillary_variables: [] wavelength: None calibration: reflectance units: % _satpy_id: DataID(name='natural_color', resolution=742) name: natural_color prerequisites: [DataQuery(name='M10', modifiers=('sunz_correcte... optional_prerequisites: [DataQuery(name='I01', modifiers=('sunz_correcte... mode: RGB <xarray.DataArray 'where-7e75f1ee81973adfbe1fa6dd2bf51ac5' (bands: 3, y: 3248, x: 3200)> dask.array<where, shape=(3, 3248, 3200), dtype=float32, chunksize=(1, 3248, 3200), chunktype=numpy.ndarray> Coordinates: crs object +proj=latlong +datum=WGS84 +ellps=WGS84 +type=crs * bands (bands) <U1 'R' 'G' 'B' Dimensions without coordinates: y, x Attributes: start_orbit: 48084 start_time: 2021-02-06 13:06:00 sensor: VIIRS end_time: 2021-02-06 13:12:00 standard_name: dust resolution: 742 end_orbit: 48084 valid_min: 0 _FillValue: 65535 rows_per_scan: 16 file_units: Watts/meter^2/steradian/micrometer platform_name: Suomi-NPP valid_max: 65527 reader: viirs_l1b file_type: vl1bm area: Shape: (3248, 3200)\nLons: <xarray.DataArray 'lo... shape: (3248, 3200) coordinates: ('m_lon', 'm_lat') ancillary_variables: [] wavelength: None optional_datasets: [] name: dust _satpy_id: DataID(name='dust', resolution=742) prerequisites: [DataQuery(name='_dust_dep_0'), DataQuery(name='... optional_prerequisites: [] mode: RGB
Often, you might want to highlight a specific geographical region. Let us generate a geographical subset around Sardinia, Italy where a dust cloud is visible. You can do this with the function stored in the coord2area_def.py
script, which converts human coordinates (longitude and latitude) to an area definition.
We need to define the following arguments:
name
:the name of the area definition, set this to sardinia_area_1km
proj
: the projection, set this to laea
which stands for the Lambert azimuthal equal-area projectionmin_lat
: the minimum latitude value, set this to 35
max_lat
: the maximum latitude value, set this to 45
min_lon
: the minimum longitude value, set this to 0
max_lon
: the maximum longitude value, set this to 15
resolution(km)
: the resolution in kilometres, set this to 1
Afterwards, you can visualize the resampled image with the function show()
.
%run coord2area_def.py sardinia_area_1km laea 35 45 0 15 1
### +proj=laea +lat_0=40.0 +lon_0=7.5 +ellps=WGS84 sardinia_area_1km: description: sardinia_area_1km projection: proj: laea ellps: WGS84 lat_0: 40.0 lon_0: 7.5 shape: height: 1107 width: 1368 area_extent: lower_left_xy: [-684216.045232, -526823.964278] upper_right_xy: [684216.045232, 580660.441353]
From the values generated by coord2area_def.py
, we copy and paste several into the template below.
We need to define the following arguments in the code block template below:
area_id
(string): the name of the area definition, set this to sardinia_area_1km
x_size
(integer): the number of values for the width, set this to the value of the shape width
, which is 1368
y_size
(integer): the number of values for the height, set this to the value of the shape height
, which is 1107
area_extent
(set of coordinates in brackets): the extent of the map is defined by 2 sets of coordinates, within a set of brackets ()
paste in the values of the lower_left_xy
from the area_extent above, followed by the upper_right_xy
values. You should end up with (-684216.045232, -526823.964278, 684216.045232, 580660.441353)
.projection
(string): the projection, paste in the first line after ###
starting with +proj
description
(string): Give this a generic name for the region,proj_id
(string): A recommended format is the projection short code followed by lat_0 and lon_0, e.g. laea_40.0_7.5
You should end up with the following code block.
from pyresample import get_area_def
area_id = 'sardinia_area_1km'
x_size = 1368
y_size = 1107
area_extent = (-684216.045232, -526823.964278, 684216.045232, 580660.441353)
projection = '+proj=laea +lat_0=40.0 +lon_0=7.5 +ellps=WGS84'
description = "Southern Europe"
proj_id = 'laea_40.0_7.5'
areadef = get_area_def(area_id, description, proj_id, projection,x_size, y_size, area_extent)
Next, you can use the area definition above in order to resample the loaded Scene object. You can use the function resample()
to do so.
scn_resample_nc = scn.resample(areadef)
Afterwards, you can visualize the resampled natural_color
RGB with the function show()
. You see see dust intrusions over the Mediterranean sea.
scn_resample_nc.show('natural_color')
SatPy's built-in visualization function is nice, but often you want to make use of additonal features, such as country borders. The library Cartopy offers powerful functions that enable the visualization of geospatial data in different projections and to add additional features to a plot. Below, we will show you how you can visualize the natural_color
composite with the two Python packages matplotlib and Cartopy.
As a first step, we have to convert the Scene
object into a numpy array. The numpy array additionally needs to be transposed to a shape that can be interpreted by matplotlib's function imshow()
: (M,N,3). You can convert a Scene object into a numpy.array
object with the function np.asarray()
. We have to transpose the array and add index=0 on index position 3.
image = np.asarray(scn_resample_nc['natural_color']).transpose(1,2,0)
The next step is then to replace all nan values with 0. You can do this with the numpy function nan_to_num()
. In a subsequent step, we then scale the values to the range between 0 and 1, clipping the lower and upper percentiles so that a potential contrast decrease caused by outliers is eliminated.
image = np.nan_to_num(image)
image = np.interp(image, (np.percentile(image,1), np.percentile(image,99)), (0, 1))
Let us now also define a variable for the coordinate reference system
. We take the area
attribute from she scn_resample_nc
Scene and convert it with the function to_cartopy_crs()
into a format Cartopy can read. We will use the crs
information for plotting.
crs = scn_resample_nc['natural_color'].attrs['area'].to_cartopy_crs()
Now, we can visualize the natural_color
composite. The plotting code can be divided in four main parts:
imshow()
# ... and use it to generate an axes in our figure with the same CRS
fig = plt.figure(figsize=(16, 10))
ax = fig.add_subplot(1, 1, 1, projection=crs)
# Now we can add some coastlines...
ax.coastlines(resolution="10m", color="white")
# ... and a lat/lon grid:
gl = ax.gridlines(draw_labels=True, linestyle='--', xlocs=range(-10,11,5), ylocs=range(30,50,5))
gl.top_labels=False
gl.right_labels=False
gl.xformatter=LONGITUDE_FORMATTER
gl.yformatter=LATITUDE_FORMATTER
gl.xlabel_style={'size':14}
gl.ylabel_style={'size':14}
# In the end, we can plot our image data...
ax.imshow(image, transform=crs, extent=crs.bounds, origin="upper")
# and add a title to our plot
plt.title("Natural color composite around Sardinia, Italy, recorded by VIIRS at " + scn_resample_nc.start_time.strftime("%Y-%m-%d %H:%M"), fontsize=20, pad=20.0)
Text(0.5, 1.0, 'Natural color composite around Sardinia, Italy, recorded by VIIRS at 2021-02-06 13:06')
In a final step, we would like to have a closer look at the Dust RGB composite, whose primary aim is to support the detection of dust in the atmosphere. The dust composite makes use of three window channels:
IR12.0-IR10.8
IR10.8-IR.8.7
IR10.8
Since we already loaded the dust composite to our Scene object scn_resample_nc
, we can visualize the dust composite with Satpy's built-in visualization function show()
. The colours of the dust RGB can be interpreted as follows:
Magenta
: Dust or ash cloudsBlack
: Cirrus cloudsDark red
: Thick, high and cold ice cloudsYellow
: Thick mid-level cloudsDarkblue
: Humid air in lower levelsLilac
: Dry air in lower levelsGet more information about the Dust RGB in the SEVIRI Dust RGB Quick Guide.
The dust composite of 6 February at 13:06 UTC clearly show a dust cloud below the thicker cold ice cloud over northern Italy.
scn_resample_nc.show('dust')
VNP02IMG - VIIRS/NPP Imagery Resolution 6-Min L1B Swath 375m distributed from NASA LAADS DAAC. Available on-line https://ladsweb.modaps.eosdis.nasa.gov/search/order/2/VNP02IMG--5200. doi:10.5067/VIIRS/VNP02IMG.002
VNP02DNB - VIIRS/NPP Day/Night Band 6-Min L1B Swath 750m distributed from NASA LAADS DAAC. Available on-line https://ladsweb.modaps.eosdis.nasa.gov/search/order/2/VNP02DNB--5200. doi:10.5067/VIIRS/VNP02IMG.002
VNP02MOD - VIIRS/NPP Moderate Resolution 6-Min L1B Swath 750m distributed from NASA LAADS DAAC. Available on-line https://ladsweb.modaps.eosdis.nasa.gov/search/order/2/VNP02MOD--5200. doi:10.5067/VIIRS/VNP02MOD.002
VNP03IMG - VIIRS/NPP Imagery Resolution Terrain-Corrected Geolocation 6-Min L1 Swath 375m distributed from NASA LAADS DAAC. Available on-line https://ladsweb.modaps.eosdis.nasa.gov/search/order/2/VNP03IMG--5200. doi:10.5067/VIIRS/VNP03IMG.002
VNP03DNB - VIIRS/NPP Day/Night Band Moderate Resolution Terrain-Corrected Geolocation 6-Min L1 Swath 750m distributed from NASA LAADS DAAC. Available on-line https://ladsweb.modaps.eosdis.nasa.gov/search/order/2/VNP02DNB--5200. doi:10.5067/VIIRS/VNP02DNB.002
VNP03MOD - VIIRS/NPP Moderate Resolution Terrain-Corrected Geolocation 6-Min L1 Swath 750m distributed from NASA LAADS DAAC. Available on-line https://ladsweb.modaps.eosdis.nasa.gov/missions-and-measurements/products/VNP03MOD. doi:10.5067/VIIRS/VNP03MOD.002
Some code in this notebook was adapted from the following sources:
This project is licensed under GNU General Public License v3.0 only and is developed under a Copernicus contract.