Python Tutorial - How to work with CoastWatch data in Python

This tutorial is based on the OceanWatch tutorial meterial edited with Great Lakes data. This tutorial will show the steps to grab data in ERDDAP from Python, how to work with NetCDF files in Python and how to make some maps and time-series water surface temperature (sst) in Lake Erie.

1. Downlading data from Python

Because ERDDAP includes RESTful services, you can download data listed on any ERDDAP platform from Python using the URL structure. For example, the following page allows you to subset daily water surface temperature data from the dataset GLSEA_ACSPO_GCS

In this specific example, the URL we generated is :

https://apps.glerl.noaa.gov/erddap/griddap/GLSEA_ACSPO_GCS.nc?sst%5B(2023-06-01T12:00:00Z):1:(2023-06-30T12:00:00Z)%5D%5B(41):1:(43)%5D%5B(-83.5):1:(-78.5)%5D

In Python, run the following to download the data using the generated URL. Note: replace coastwatch.glerl.noaa.gov with apps.glerl.noaa.gov) :

import urllib.request

url="https://apps.glerl.noaa.gov/erddap/griddap/GLSEA_ACSPO_GCS.nc?sst%5B(2023-06-01T12:00:00Z):1:(2023-06-30T12:00:00Z)%5D%5B(41):1:(43)%5D%5B(-83.5):1:(-78.5)%5D"
urllib.request.urlretrieve(url, "e_sst.nc")

('e_sst.nc', <http.client.HTTPMessage at 0x1b7d5152610>)

2. Importing NetCDF4 data in Python

Now that we’ve downloaded the data locally, we can import it and extract our variables of interest.

The xarray package makes it very convenient to work with NetCDF files. Documentation is available here: http://xarray.pydata.org/en/stable/why-xarray.html

import xarray as xr
import netCDF4 as nc

- Open the file and load it as an xarray dataset:

ds = xr.open_dataset('e_sst.nc',decode_cf=False)
#ds = xr.open_dataset('e_sst.nc')

- Examine the data structure:

ds

<xarray.Dataset>
Dimensions:    (time: 30, latitude: 143, longitude: 358)
Coordinates:
  * time       (time) float64 1.686e+09 1.686e+09 ... 1.688e+09 1.688e+09
  * latitude   (latitude) float64 41.01 41.02 41.03 41.05 ... 42.97 42.98 43.0
  * longitude  (longitude) float64 -83.51 -83.49 -83.48 ... -78.53 -78.52 -78.5
Data variables:
    sst        (time, latitude, longitude) float32 ...
Attributes: (12/37)
    cdm_data_type:                             Grid
    Conventions:                               CF-1.6, COARDS, ACDD-1.3
    Easternmost_Easting:                       -78.5025420282523
    GDAL:                                      GDAL 3.4.3, released 2022/04/22
    geospatial_lat_max:                        42.9955492889468
    geospatial_lat_min:                        41.0053458284746
    ...                                        ...
    summary:                                   Sea Surface Temperature (SST) ...
    testOutOfDate:                             now-3days
    time_coverage_end:                         2023-06-30T12:00:00Z
    time_coverage_start:                       2023-06-01T12:00:00Z
    title:                                     Sea Surface Temperature (SST) ...
    Westernmost_Easting:                       -83.5060817140874

xarray.Dataset

• Dimensions:
  • time: 30
  • latitude: 143
  • longitude: 358
• Coordinates: (3)
  • time
    (time)
    float64
    1.686e+09 1.686e+09 ... 1.688e+09

    _CoordinateAxisType :

    Time

    actual_range :

    [1.6856208e+09 1.6881264e+09]

    axis :

    T

    calendar :

    Gregorian

    ioos_category :

    Time

    long_name :

    Easten Time

    standard_name :

    time

    time_origin :

    01-JAN-1970 00:00:00

    units :

    seconds since 1970-01-01T00:00:00Z

    array([1.685621e+09, 1.685707e+09, 1.685794e+09, 1.685880e+09, 1.685966e+09,
           1.686053e+09, 1.686139e+09, 1.686226e+09, 1.686312e+09, 1.686398e+09,
           1.686485e+09, 1.686571e+09, 1.686658e+09, 1.686744e+09, 1.686830e+09,
           1.686917e+09, 1.687003e+09, 1.687090e+09, 1.687176e+09, 1.687262e+09,
           1.687349e+09, 1.687435e+09, 1.687522e+09, 1.687608e+09, 1.687694e+09,
           1.687781e+09, 1.687867e+09, 1.687954e+09, 1.688040e+09, 1.688126e+09])

  • latitude
    (latitude)
    float64
    41.01 41.02 41.03 ... 42.98 43.0

    _CoordinateAxisType :

    Lat

    actual_range :

    [41.00534583 42.99554929]

    axis :

    Y

    ioos_category :

    Location

    long_name :

    Latitude

    standard_name :

    latitude

    units :

    degrees_north

    array([41.005346, 41.019361, 41.033377, 41.047392, 41.061408, 41.075423,
           41.089439, 41.103454, 41.11747 , 41.131485, 41.145501, 41.159517,
           41.173532, 41.187548, 41.201563, 41.215579, 41.229594, 41.24361 ,
           41.257625, 41.271641, 41.285656, 41.299672, 41.313687, 41.327703,
           41.341718, 41.355734, 41.369749, 41.383765, 41.39778 , 41.411796,
           41.425811, 41.439827, 41.453842, 41.467858, 41.481873, 41.495889,
           41.509904, 41.52392 , 41.537935, 41.551951, 41.565967, 41.579982,
           41.593998, 41.608013, 41.622029, 41.636044, 41.65006 , 41.664075,
           41.678091, 41.692106, 41.706122, 41.720137, 41.734153, 41.748168,
           41.762184, 41.776199, 41.790215, 41.80423 , 41.818246, 41.832261,
           41.846277, 41.860292, 41.874308, 41.888323, 41.902339, 41.916354,
           41.93037 , 41.944385, 41.958401, 41.972417, 41.986432, 42.000448,
           42.014463, 42.028479, 42.042494, 42.05651 , 42.070525, 42.084541,
           42.098556, 42.112572, 42.126587, 42.140603, 42.154618, 42.168634,
           42.182649, 42.196665, 42.21068 , 42.224696, 42.238711, 42.252727,
           42.266742, 42.280758, 42.294773, 42.308789, 42.322804, 42.33682 ,
           42.350835, 42.364851, 42.378867, 42.392882, 42.406898, 42.420913,
           42.434929, 42.448944, 42.46296 , 42.476975, 42.490991, 42.505006,
           42.519022, 42.533037, 42.547053, 42.561068, 42.575084, 42.589099,
           42.603115, 42.61713 , 42.631146, 42.645161, 42.659177, 42.673192,
           42.687208, 42.701223, 42.715239, 42.729254, 42.74327 , 42.757285,
           42.771301, 42.785317, 42.799332, 42.813348, 42.827363, 42.841379,
           42.855394, 42.86941 , 42.883425, 42.897441, 42.911456, 42.925472,
           42.939487, 42.953503, 42.967518, 42.981534, 42.995549])

  • longitude
    (longitude)
    float64
    -83.51 -83.49 ... -78.52 -78.5

    _CoordinateAxisType :

    Lon

    actual_range :

    [-83.50608171 -78.50254203]

    axis :

    X

    ioos_category :

    Location

    long_name :

    Longitude

    standard_name :

    longitude

    units :

    degrees_east

    array([-83.506082, -83.492066, -83.478051, ..., -78.530573, -78.516558,
           -78.502542])

• Data variables: (1)
  • sst
    (time, latitude, longitude)
    float32
    ...

    _FillValue :

    -99999.0

    colorBarMaximum :

    32.0

    colorBarMinimum :

    0.0

    ioos_category :

    Temperature

    long_name :

    Temperature

    standard_name :

    sea_water_temperature

    units :

    degree_C

    [1535820 values with dtype=float32]

• Indexes: (3)
  • time
    PandasIndex

    PandasIndex(Float64Index([1685620800.0, 1685707200.0, 1685793600.0, 1685880000.0,
                  1685966400.0, 1686052800.0, 1686139200.0, 1686225600.0,
                  1686312000.0, 1686398400.0, 1686484800.0, 1686571200.0,
                  1686657600.0, 1686744000.0, 1686830400.0, 1686916800.0,
                  1687003200.0, 1687089600.0, 1687176000.0, 1687262400.0,
                  1687348800.0, 1687435200.0, 1687521600.0, 1687608000.0,
                  1687694400.0, 1687780800.0, 1687867200.0, 1687953600.0,
                  1688040000.0, 1688126400.0],
                 dtype='float64', name='time'))

  • latitude
    PandasIndex

    PandasIndex(Float64Index([41.0053458284746, 41.0193613458019, 41.0333768631292,
                  41.0473923804564, 41.0614078977837,  41.075423415111,
                  41.0894389324382, 41.1034544497655, 41.1174699670928,
                    41.13148548442,
                  ...
                  42.8694096330014, 42.8834251503287, 42.8974406676559,
                  42.9114561849832, 42.9254717023105, 42.9394872196377,
                   42.953502736965, 42.9675182542923, 42.9815337716195,
                  42.9955492889468],
                 dtype='float64', name='latitude', length=143))

  • longitude
    PandasIndex

    PandasIndex(Float64Index([-83.5060817140874, -83.4920661967601, -83.4780506794328,
                  -83.4640351621056, -83.4500196447783,  -83.436004127451,
                  -83.4219886101238, -83.4079730927965, -83.3939575754692,
                  -83.3799420581419,
                  ...
                  -78.6286816841978, -78.6146661668705, -78.6006506495432,
                   -78.586635132216, -78.5726196148887, -78.5586040975614,
                  -78.5445885802342, -78.5305730629069, -78.5165575455796,
                  -78.5025420282523],
                 dtype='float64', name='longitude', length=358))

• Attributes: (37)

  cdm_data_type :

  Grid

  Conventions :

  CF-1.6, COARDS, ACDD-1.3

  Easternmost_Easting :

  -78.5025420282523

  GDAL :

  GDAL 3.4.3, released 2022/04/22

  geospatial_lat_max :

  42.9955492889468

  geospatial_lat_min :

  41.0053458284746

  geospatial_lat_resolution :

  0.014015517327269056

  geospatial_lat_units :

  degrees_north

  geospatial_lon_max :

  -78.5025420282523

  geospatial_lon_min :

  -83.5060817140874

  geospatial_lon_resolution :

  0.01401551732726889

  geospatial_lon_units :

  degrees_east

  grid_mapping_GeoTransform :

  -92.42695849289404 0.01401551732726498 0 50.61298295630733 0 -0.01401551732726498

  grid_mapping_inverse_flattening :

  298.257223563

  grid_mapping_long_name :

  CRS definition

  grid_mapping_longitude_of_prime_meridian :

  0.0

  grid_mapping_name :

  latitude_longitude

  grid_mapping_semi_major_axis :

  6378137.0

  grid_mapping_spatial_ref :

  GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]

  history :

  Sun Apr 28 03:50:43 2024: ncrename -v Band1,sst /home/cw/opr/glsea3_b/mk_nc/output/2024_108_glsea_sst.nc Sun Apr 28 03:50:43 2024: GDAL Create( /home/cw/opr/glsea3_b/mk_nc/output/2024_108_glsea_sst.nc, ... ) 2024-04-29T15:13:03Z (local files) 2024-04-29T15:13:03Z https://apps.glerl.noaa.gov/erddap/griddap/GLSEA_ACSPO_GCS.nc?sst%5B(2023-06-01T12:00:00Z):1:(2023-06-30T12:00:00Z)%5D%5B(41):1:(43)%5D%5B(-83.5):1:(-78.5)%5D

  infoUrl :

  https://coastwatch.glerl.noaa.gov/glsea/glsea.html

  institution :

  CoastWatch Great Lakes Node

  keywords :

  ACSPO, analysis, data, earth, Earth Science > Oceans > Ocean Temperature > Water Temperature, environmental, glsea, great, great lakes, lakes, ocean, oceans, science, sea, sea_water_temperature, seawater, sst, surface, temperature, time, water

  keywords_vocabulary :

  GCMD Science Keywords

  license :

  The data may be used and redistributed for free but is not intended for legal use, since it may contain inaccuracies. Neither the data Contributor, ERD, NOAA, nor the United States Government, nor any of their employees or contractors, makes any warranty, express or implied, including warranties of merchantability and fitness for a particular purpose, or assumes any legal liability for the accuracy, completeness, or usefulness, of this information.

  NCO :

  netCDF Operators version 4.7.5 (Homepage = http://nco.sf.net, Code = https://github.com/nco/nco)

  Northernmost_Northing :

  42.9955492889468

  source :

  2024_108_glsea.asc

  sourceUrl :

  (local files)

  Southernmost_Northing :

  41.0053458284746

  standard_name_vocabulary :

  CF Standard Name Table v70

  summary :

  Sea Surface Temperature (SST) from Great Lakes Surface Environmental Analysis (GLSEA - ACSPO)

  testOutOfDate :

  now-3days

  time_coverage_end :

  2023-06-30T12:00:00Z

  time_coverage_start :

  2023-06-01T12:00:00Z

  title :

  Sea Surface Temperature (SST) from Great Lakes Surface Environmental Analysis (ACSPO GLSEA), The ACSPO GLSEA is the developed product that uses the Advanced Clear-sky Processor for Oceans (ACSPO) L3S-LEO SST as the input data source. The L3S-LEO SST is the level 3 Super Collated data from Low Earth Orbit satellites sensor including NPP, N20, and MetOp A/B/C. Geodetic coordinate system (LAT, LON), 2006-present

  Westernmost_Easting :

  -83.5060817140874

print(ds)

<xarray.Dataset>
Dimensions:    (time: 30, latitude: 143, longitude: 358)
Coordinates:
  * time       (time) float64 1.686e+09 1.686e+09 ... 1.688e+09 1.688e+09
  * latitude   (latitude) float64 41.01 41.02 41.03 41.05 ... 42.97 42.98 43.0
  * longitude  (longitude) float64 -83.51 -83.49 -83.48 ... -78.53 -78.52 -78.5
Data variables:
    sst        (time, latitude, longitude) float32 ...
Attributes: (12/37)
    cdm_data_type:                             Grid
    Conventions:                               CF-1.6, COARDS, ACDD-1.3
    Easternmost_Easting:                       -78.5025420282523
    GDAL:                                      GDAL 3.4.3, released 2022/04/22
    geospatial_lat_max:                        42.9955492889468
    geospatial_lat_min:                        41.0053458284746
    ...                                        ...
    summary:                                   Sea Surface Temperature (SST) ...
    testOutOfDate:                             now-3days
    time_coverage_end:                         2023-06-30T12:00:00Z
    time_coverage_start:                       2023-06-01T12:00:00Z
    title:                                     Sea Surface Temperature (SST) ...
    Westernmost_Easting:                       -83.5060817140874

- Examine which coordinates and variables are included in the dataset:

#ds.dims

ds.coords

Coordinates:
  * time       (time) float64 1.686e+09 1.686e+09 ... 1.688e+09 1.688e+09
  * latitude   (latitude) float64 41.01 41.02 41.03 41.05 ... 42.97 42.98 43.0
  * longitude  (longitude) float64 -83.51 -83.49 -83.48 ... -78.53 -78.52 -78.5

ds.data_vars

Data variables:
    sst      (time, latitude, longitude) float32 ...

ds.attrs

{'cdm_data_type': 'Grid',
 'Conventions': 'CF-1.6, COARDS, ACDD-1.3',
 'Easternmost_Easting': -78.5025420282523,
 'GDAL': 'GDAL 3.4.3, released 2022/04/22',
 'geospatial_lat_max': 42.9955492889468,
 'geospatial_lat_min': 41.0053458284746,
 'geospatial_lat_resolution': 0.014015517327269056,
 'geospatial_lat_units': 'degrees_north',
 'geospatial_lon_max': -78.5025420282523,
 'geospatial_lon_min': -83.5060817140874,
 'geospatial_lon_resolution': 0.01401551732726889,
 'geospatial_lon_units': 'degrees_east',
 'grid_mapping_GeoTransform': '-92.42695849289404 0.01401551732726498 0 50.61298295630733 0 -0.01401551732726498',
 'grid_mapping_inverse_flattening': 298.257223563,
 'grid_mapping_long_name': 'CRS definition',
 'grid_mapping_longitude_of_prime_meridian': 0.0,
 'grid_mapping_name': 'latitude_longitude',
 'grid_mapping_semi_major_axis': 6378137.0,
 'grid_mapping_spatial_ref': 'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]',
 'history': 'Sun Apr 28 03:50:43 2024: ncrename -v Band1,sst /home/cw/opr/glsea3_b/mk_nc/output/2024_108_glsea_sst.nc\nSun Apr 28 03:50:43 2024: GDAL Create( /home/cw/opr/glsea3_b/mk_nc/output/2024_108_glsea_sst.nc, ... )\n2024-04-29T15:13:03Z (local files)\n2024-04-29T15:13:03Z https://apps.glerl.noaa.gov/erddap/griddap/GLSEA_ACSPO_GCS.nc?sst%5B(2023-06-01T12:00:00Z):1:(2023-06-30T12:00:00Z)%5D%5B(41):1:(43)%5D%5B(-83.5):1:(-78.5)%5D',
 'infoUrl': 'https://coastwatch.glerl.noaa.gov/glsea/glsea.html',
 'institution': 'CoastWatch Great Lakes Node',
 'keywords': 'ACSPO, analysis, data, earth, Earth Science > Oceans > Ocean Temperature > Water Temperature, environmental, glsea, great, great lakes, lakes, ocean, oceans, science, sea, sea_water_temperature, seawater, sst, surface, temperature, time, water',
 'keywords_vocabulary': 'GCMD Science Keywords',
 'license': 'The data may be used and redistributed for free but is not intended\nfor legal use, since it may contain inaccuracies. Neither the data\nContributor, ERD, NOAA, nor the United States Government, nor any\nof their employees or contractors, makes any warranty, express or\nimplied, including warranties of merchantability and fitness for a\nparticular purpose, or assumes any legal liability for the accuracy,\ncompleteness, or usefulness, of this information.',
 'NCO': 'netCDF Operators version 4.7.5 (Homepage = http://nco.sf.net, Code = https://github.com/nco/nco)',
 'Northernmost_Northing': 42.9955492889468,
 'source': '2024_108_glsea.asc',
 'sourceUrl': '(local files)',
 'Southernmost_Northing': 41.0053458284746,
 'standard_name_vocabulary': 'CF Standard Name Table v70',
 'summary': 'Sea Surface Temperature (SST) from Great Lakes Surface Environmental Analysis (GLSEA - ACSPO)',
 'testOutOfDate': 'now-3days',
 'time_coverage_end': '2023-06-30T12:00:00Z',
 'time_coverage_start': '2023-06-01T12:00:00Z',
 'title': 'Sea Surface Temperature (SST) from Great Lakes Surface Environmental Analysis (ACSPO GLSEA), The ACSPO GLSEA is the developed product that uses the Advanced Clear-sky Processor for Oceans (ACSPO) L3S-LEO SST as the input data source. The L3S-LEO SST is the level 3 Super Collated data from Low Earth Orbit satellites sensor including NPP, N20, and MetOp A/B/C. Geodetic coordinate system (LAT, LON), 2006-present',
 'Westernmost_Easting': -83.5060817140874}

- Examine the structure of sst:

ds.sst.shape

(30, 143, 358)

Our dataset is a 3-D array with 143 rows corresponding to latitudes and 358 columns corresponding to longitudes, for each of the 30 time steps.

- Get the dates for each time step:

ds.time

<xarray.DataArray 'time' (time: 30)>
array([1.685621e+09, 1.685707e+09, 1.685794e+09, 1.685880e+09, 1.685966e+09,
       1.686053e+09, 1.686139e+09, 1.686226e+09, 1.686312e+09, 1.686398e+09,
       1.686485e+09, 1.686571e+09, 1.686658e+09, 1.686744e+09, 1.686830e+09,
       1.686917e+09, 1.687003e+09, 1.687090e+09, 1.687176e+09, 1.687262e+09,
       1.687349e+09, 1.687435e+09, 1.687522e+09, 1.687608e+09, 1.687694e+09,
       1.687781e+09, 1.687867e+09, 1.687954e+09, 1.688040e+09, 1.688126e+09])
Coordinates:
  * time     (time) float64 1.686e+09 1.686e+09 ... 1.688e+09 1.688e+09
Attributes:
    _CoordinateAxisType:  Time
    actual_range:         [1.6856208e+09 1.6881264e+09]
    axis:                 T
    calendar:             Gregorian
    ioos_category:        Time
    long_name:            Easten Time
    standard_name:        time
    time_origin:          01-JAN-1970 00:00:00
    units:                seconds since 1970-01-01T00:00:00Z

xarray.DataArray

'time'

• time: 30

• 1.686e+09 1.686e+09 1.686e+09 ... 1.688e+09 1.688e+09 1.688e+09

  array([1.685621e+09, 1.685707e+09, 1.685794e+09, 1.685880e+09, 1.685966e+09,
         1.686053e+09, 1.686139e+09, 1.686226e+09, 1.686312e+09, 1.686398e+09,
         1.686485e+09, 1.686571e+09, 1.686658e+09, 1.686744e+09, 1.686830e+09,
         1.686917e+09, 1.687003e+09, 1.687090e+09, 1.687176e+09, 1.687262e+09,
         1.687349e+09, 1.687435e+09, 1.687522e+09, 1.687608e+09, 1.687694e+09,
         1.687781e+09, 1.687867e+09, 1.687954e+09, 1.688040e+09, 1.688126e+09])

• Coordinates: (1)
  • time
    (time)
    float64
    1.686e+09 1.686e+09 ... 1.688e+09

    _CoordinateAxisType :

    Time

    actual_range :

    [1.6856208e+09 1.6881264e+09]

    axis :

    T

    calendar :

    Gregorian

    ioos_category :

    Time

    long_name :

    Easten Time

    standard_name :

    time

    time_origin :

    01-JAN-1970 00:00:00

    units :

    seconds since 1970-01-01T00:00:00Z

    array([1.685621e+09, 1.685707e+09, 1.685794e+09, 1.685880e+09, 1.685966e+09,
           1.686053e+09, 1.686139e+09, 1.686226e+09, 1.686312e+09, 1.686398e+09,
           1.686485e+09, 1.686571e+09, 1.686658e+09, 1.686744e+09, 1.686830e+09,
           1.686917e+09, 1.687003e+09, 1.687090e+09, 1.687176e+09, 1.687262e+09,
           1.687349e+09, 1.687435e+09, 1.687522e+09, 1.687608e+09, 1.687694e+09,
           1.687781e+09, 1.687867e+09, 1.687954e+09, 1.688040e+09, 1.688126e+09])

• Indexes: (1)
  • time
    PandasIndex

    PandasIndex(Float64Index([1685620800.0, 1685707200.0, 1685793600.0, 1685880000.0,
                  1685966400.0, 1686052800.0, 1686139200.0, 1686225600.0,
                  1686312000.0, 1686398400.0, 1686484800.0, 1686571200.0,
                  1686657600.0, 1686744000.0, 1686830400.0, 1686916800.0,
                  1687003200.0, 1687089600.0, 1687176000.0, 1687262400.0,
                  1687348800.0, 1687435200.0, 1687521600.0, 1687608000.0,
                  1687694400.0, 1687780800.0, 1687867200.0, 1687953600.0,
                  1688040000.0, 1688126400.0],
                 dtype='float64', name='time'))

• Attributes: (9)

  _CoordinateAxisType :

  Time

  actual_range :

  [1.6856208e+09 1.6881264e+09]

  axis :

  T

  calendar :

  Gregorian

  ioos_category :

  Time

  long_name :

  Easten Time

  standard_name :

  time

  time_origin :

  01-JAN-1970 00:00:00

  units :

  seconds since 1970-01-01T00:00:00Z

ds.time.attrs

{'_CoordinateAxisType': 'Time',
 'actual_range': array([1.6856208e+09, 1.6881264e+09]),
 'axis': 'T',
 'calendar': 'Gregorian',
 'ioos_category': 'Time',
 'long_name': 'Easten Time',
 'standard_name': 'time',
 'time_origin': '01-JAN-1970 00:00:00',
 'units': 'seconds since 1970-01-01T00:00:00Z'}

the time units is seconds, we need to convert the seconds to dates.

dates=nc.num2date(ds.time,ds.time.units,only_use_cftime_datetimes=False, 
                        only_use_python_datetimes=True )
dates

array([real_datetime(2023, 6, 1, 12, 0), real_datetime(2023, 6, 2, 12, 0),
       real_datetime(2023, 6, 3, 12, 0), real_datetime(2023, 6, 4, 12, 0),
       real_datetime(2023, 6, 5, 12, 0), real_datetime(2023, 6, 6, 12, 0),
       real_datetime(2023, 6, 7, 12, 0), real_datetime(2023, 6, 8, 12, 0),
       real_datetime(2023, 6, 9, 12, 0),
       real_datetime(2023, 6, 10, 12, 0),
       real_datetime(2023, 6, 11, 12, 0),
       real_datetime(2023, 6, 12, 12, 0),
       real_datetime(2023, 6, 13, 12, 0),
       real_datetime(2023, 6, 14, 12, 0),
       real_datetime(2023, 6, 15, 12, 0),
       real_datetime(2023, 6, 16, 12, 0),
       real_datetime(2023, 6, 17, 12, 0),
       real_datetime(2023, 6, 18, 12, 0),
       real_datetime(2023, 6, 19, 12, 0),
       real_datetime(2023, 6, 20, 12, 0),
       real_datetime(2023, 6, 21, 12, 0),
       real_datetime(2023, 6, 22, 12, 0),
       real_datetime(2023, 6, 23, 12, 0),
       real_datetime(2023, 6, 24, 12, 0),
       real_datetime(2023, 6, 25, 12, 0),
       real_datetime(2023, 6, 26, 12, 0),
       real_datetime(2023, 6, 27, 12, 0),
       real_datetime(2023, 6, 28, 12, 0),
       real_datetime(2023, 6, 29, 12, 0),
       real_datetime(2023, 6, 30, 12, 0)], dtype=object)

The datetime object includes year, month, hour, minutes, eg. 2021, 6, 12, 0.

Working with the extracted data

Creating a map for one time step

Let’s create a map of SST for June 1, 2021 (our first time step).

import pandas as pd
import numpy as np
from matplotlib import pyplot as plt
from matplotlib.colors import LinearSegmentedColormap
#np.warnings.filterwarnings('ignore')

- Examine the values of sst:

ds.sst.values

array([[[-99999., -99999., -99999., ..., -99999., -99999., -99999.],
        [-99999., -99999., -99999., ..., -99999., -99999., -99999.],
        [-99999., -99999., -99999., ..., -99999., -99999., -99999.],
        ...,
        [-99999., -99999., -99999., ..., -99999., -99999., -99999.],
        [-99999., -99999., -99999., ..., -99999., -99999., -99999.],
        [-99999., -99999., -99999., ..., -99999., -99999., -99999.]],

       [[-99999., -99999., -99999., ..., -99999., -99999., -99999.],
        [-99999., -99999., -99999., ..., -99999., -99999., -99999.],
        [-99999., -99999., -99999., ..., -99999., -99999., -99999.],
        ...,
        [-99999., -99999., -99999., ..., -99999., -99999., -99999.],
        [-99999., -99999., -99999., ..., -99999., -99999., -99999.],
        [-99999., -99999., -99999., ..., -99999., -99999., -99999.]],

       [[-99999., -99999., -99999., ..., -99999., -99999., -99999.],
        [-99999., -99999., -99999., ..., -99999., -99999., -99999.],
        [-99999., -99999., -99999., ..., -99999., -99999., -99999.],
        ...,
        [-99999., -99999., -99999., ..., -99999., -99999., -99999.],
        [-99999., -99999., -99999., ..., -99999., -99999., -99999.],
        [-99999., -99999., -99999., ..., -99999., -99999., -99999.]],

       ...,

       [[-99999., -99999., -99999., ..., -99999., -99999., -99999.],
        [-99999., -99999., -99999., ..., -99999., -99999., -99999.],
        [-99999., -99999., -99999., ..., -99999., -99999., -99999.],
        ...,
        [-99999., -99999., -99999., ..., -99999., -99999., -99999.],
        [-99999., -99999., -99999., ..., -99999., -99999., -99999.],
        [-99999., -99999., -99999., ..., -99999., -99999., -99999.]],

       [[-99999., -99999., -99999., ..., -99999., -99999., -99999.],
        [-99999., -99999., -99999., ..., -99999., -99999., -99999.],
        [-99999., -99999., -99999., ..., -99999., -99999., -99999.],
        ...,
        [-99999., -99999., -99999., ..., -99999., -99999., -99999.],
        [-99999., -99999., -99999., ..., -99999., -99999., -99999.],
        [-99999., -99999., -99999., ..., -99999., -99999., -99999.]],

       [[-99999., -99999., -99999., ..., -99999., -99999., -99999.],
        [-99999., -99999., -99999., ..., -99999., -99999., -99999.],
        [-99999., -99999., -99999., ..., -99999., -99999., -99999.],
        ...,
        [-99999., -99999., -99999., ..., -99999., -99999., -99999.],
        [-99999., -99999., -99999., ..., -99999., -99999., -99999.],
        [-99999., -99999., -99999., ..., -99999., -99999., -99999.]]],
      dtype=float32)

ds.sst.attrs

{'_FillValue': -99999.0,
 'colorBarMaximum': 32.0,
 'colorBarMinimum': 0.0,
 'ioos_category': 'Temperature',
 'long_name': 'Temperature',
 'standard_name': 'sea_water_temperature',
 'units': 'degree_C'}

ds.sst.attrs['_FillValue']

-99999.0

#ds.sst.dims
#ds.sst.coords

- Make a new sst DataArray and replace _fillValue with NaN

#nan_sst = ds.sst.where(ds.sst.values != -99999.0)
nan_sst = ds.sst.where(ds.sst.values != ds.sst.attrs['_FillValue'])

# nan_sst[time][latitude][longitude]
#print(nan_sst[10][100][200])

print(nan_sst)

<xarray.DataArray 'sst' (time: 30, latitude: 143, longitude: 358)>
array([[[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        ...,
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan]],

       [[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        ...,
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan]],

       [[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        ...,
...
        ...,
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan]],

       [[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        ...,
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan]],

       [[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        ...,
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan]]], dtype=float32)
Coordinates:
  * time       (time) float64 1.686e+09 1.686e+09 ... 1.688e+09 1.688e+09
  * latitude   (latitude) float64 41.01 41.02 41.03 41.05 ... 42.97 42.98 43.0
  * longitude  (longitude) float64 -83.51 -83.49 -83.48 ... -78.53 -78.52 -78.5
Attributes:
    _FillValue:       -99999.0
    colorBarMaximum:  32.0
    colorBarMinimum:  0.0
    ioos_category:    Temperature
    long_name:        Temperature
    standard_name:    sea_water_temperature
    units:            degree_C

- Set some color breaks

np.nanmin(ds.sst)

-99999.0

# find min value in man_sst
np.nanmin(nan_sst)

13.25

np.nanmax(nan_sst)

23.35

levs = np.arange(13.25, 23.35, 0.05)
len(levs)

203

- Define a color palette

# init a color list
jet=["blue", "#007FFF", "cyan","#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000"]

- Set color scale using the jet palette

cm = LinearSegmentedColormap.from_list('my_jet', jet, N=len(levs))

#https://www.youtube.com/watch?v=qk0n-YaKIkY

- plot the SST map

np.linspace(-82.5,-80,num=4)

array([-82.5       , -81.66666667, -80.83333333, -80.        ])

plt.subplots(figsize=(10, 5))

#plot first sst image: nan_sst[0,:,:]
plt.contourf(nan_sst.longitude, nan_sst.latitude, nan_sst[0,:,:], levs,cmap=cm)

#plot the color scale
plt.colorbar()

#example of how to add points to the map
plt.scatter(np.linspace(-82,-80.5,num=4),np.repeat(42,4),c='black')

#example of how to add a contour line
step = np.arange(9,26, 1)

plt.contour(ds.longitude, ds.latitude, ds.sst[0,:,:],levels=step,linewidths=1)

#plot title
plt.title("Lake Erie Water Surface Temperature - " + dates[0].strftime('%b %d, %Y'))
plt.show()

Plotting a time series

Let’s pick the following box : 41.75-42.0N, 83.0-83.5W. We are going to generate a time series of mean SST within that box.

- first, let’s subset our data:

lat_bnds, lon_bnds = [41.75, 42.0], [-83.5, -83.0]
a_sst=nan_sst.sel(latitude=slice(*lat_bnds), longitude=slice(*lon_bnds))
print(a_sst)

<xarray.DataArray 'sst' (time: 30, latitude: 17, longitude: 36)>
array([[[  nan,   nan,   nan, ..., 18.86, 18.84, 18.84],
        [  nan,   nan,   nan, ..., 18.82, 18.8 , 18.8 ],
        [  nan,   nan,   nan, ..., 18.78, 18.76, 18.76],
        ...,
        [  nan,   nan,   nan, ..., 18.5 , 18.46, 18.46],
        [  nan,   nan,   nan, ..., 18.5 , 18.46, 18.46],
        [  nan,   nan,   nan, ..., 18.53, 18.49, 18.49]],

       [[  nan,   nan,   nan, ..., 19.45, 19.44, 19.44],
        [  nan,   nan,   nan, ..., 19.42, 19.41, 19.41],
        [  nan,   nan,   nan, ..., 19.39, 19.38, 19.38],
        ...,
        [  nan,   nan,   nan, ..., 19.13, 19.1 , 19.1 ],
        [  nan,   nan,   nan, ..., 19.07, 19.06, 19.06],
        [  nan,   nan,   nan, ..., 19.06, 19.04, 19.04]],

       [[  nan,   nan,   nan, ..., 19.63, 19.61, 19.61],
        [  nan,   nan,   nan, ..., 19.6 , 19.58, 19.58],
        [  nan,   nan,   nan, ..., 19.56, 19.54, 19.54],
        ...,
...
        ...,
        [  nan,   nan,   nan, ..., 21.24, 21.31, 21.31],
        [  nan,   nan,   nan, ..., 21.24, 21.34, 21.34],
        [  nan,   nan,   nan, ..., 21.29, 21.42, 21.42]],

       [[  nan,   nan,   nan, ..., 21.8 , 21.75, 21.75],
        [  nan,   nan,   nan, ..., 21.7 , 21.67, 21.67],
        [  nan,   nan,   nan, ..., 21.6 , 21.59, 21.59],
        ...,
        [  nan,   nan,   nan, ..., 21.36, 21.41, 21.41],
        [  nan,   nan,   nan, ..., 21.36, 21.44, 21.44],
        [  nan,   nan,   nan, ..., 21.44, 21.53, 21.53]],

       [[  nan,   nan,   nan, ..., 21.93, 21.82, 21.82],
        [  nan,   nan,   nan, ..., 21.86, 21.79, 21.79],
        [  nan,   nan,   nan, ..., 21.78, 21.75, 21.75],
        ...,
        [  nan,   nan,   nan, ..., 21.49, 21.52, 21.52],
        [  nan,   nan,   nan, ..., 21.51, 21.55, 21.55],
        [  nan,   nan,   nan, ..., 21.63, 21.66, 21.66]]], dtype=float32)
Coordinates:
  * time       (time) float64 1.686e+09 1.686e+09 ... 1.688e+09 1.688e+09
  * latitude   (latitude) float64 41.76 41.78 41.79 41.8 ... 41.96 41.97 41.99
  * longitude  (longitude) float64 -83.49 -83.48 -83.46 ... -83.03 -83.02 -83.0
Attributes:
    _FillValue:       -99999.0
    colorBarMaximum:  32.0
    colorBarMinimum:  0.0
    ioos_category:    Temperature
    long_name:        Temperature
    standard_name:    sea_water_temperature
    units:            degree_C

- let’s plot the subset:

#plot first image of the a_sst array
plt.contourf(a_sst.longitude, a_sst.latitude, a_sst[0,:,:], levs,cmap=cm)
plt.colorbar()
plt.title("Subset of Lake Erie Water Surface Temperature " + dates[0].strftime('%b %d, %Y'))
plt.show()

- let’s compute the daily mean over the bounding region:

res=np.nanmean(a_sst,axis=(1,2))
res

array([19.252283, 19.84715 , 20.109388, 20.06981 , 20.01454 , 20.00379 ,
       19.958447, 19.921082, 19.963646, 19.956564, 19.885717, 19.525248,
       19.253717, 19.179811, 19.149954, 19.190142, 19.41266 , 19.634283,
       19.90713 , 20.220236, 20.525293, 20.747812, 21.015364, 21.460989,
       21.71833 , 21.710024, 21.643694, 21.593224, 21.749012, 22.0068  ],
      dtype=float32)

- let’s plot the time-series:

plt.figure(figsize=(8,4))

plt.scatter(dates,res)

degree_sign = u"\N{DEGREE SIGN}"
plt.ylabel('SST (' + degree_sign + 'C)')

plt.xlim(dates[0], dates[-1])

plt.xticks(dates,rotation=70, fontsize=10 )
plt.show()

Creating a map of average SST over a month

- let’s compute the monthly mean for the region:

import warnings

warnings.filterwarnings('ignore')
mean_sst=np.nanmean(nan_sst,axis=0)

mean_sst.shape

(143, 358)

- let’s plot the map of the average SST in the region for 2021 June:

plt.subplots(figsize=(10, 5))

plt.contourf(ds.longitude, ds.latitude, mean_sst, levs,cmap=cm)

cbar = plt.colorbar()
cbar.set_label('SST')

plt.title("Mean SST " + dates[0].strftime('%Y-%m-%d')+' - '+dates[-1].strftime('%Y-%m-%d'))
plt.show()

!jupyter nbconvert --to html GL_python_tutorial1.ipynb

[NbConvertApp] Converting notebook GL_python_tutorial1.ipynb to html
[NbConvertApp] Writing 974120 bytes to GL_python_tutorial1.html