# 1. How to work with satellite data in Python 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 of sea surface temperature (SST) around the main Hawaiian islands. This tutorial is also available as a [Jupyter notebook](https://github.com/melhawaii/python-satellite-course/blob/master/OW_Tutorial1.ipynb). You will need to install the `urllib`, `xarray` and `NetCDF4` packages if you don't already have them. ## Downloading data in Python Because ERDDAP includes RESTful services, you can download data listed on any ERDDAP platform from R using the URL structure. For example, the following page allows you to subset monthly SST data:\ \ Select your region and date range of interest, then select the '.nc' (NetCDF) file type and click on "Just Generate the URL". ![](/files/-LzF9qqbKZxTfVMPubXN) In this specific example, the URL we generated is : You can also edit this URL manually. In Python, run the following to download the data using the generated URL : `import urllib.request` \ \ `url="https://oceanwatch.pifsc.noaa.gov/erddap/griddap/CRW_sst_v1_0_monthly.nc?analysed_sst[(2018-01-01T12:00:00Z):1:(2018-12-01T12:00:00Z)][(17):1:(30)][(195):1:(210)]"` \ \ `urllib.request.urlretrieve(url, "sst.nc")` ## Importing the downloaded data in R 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:\ `import xarray as xr` \ `import netCDF4 as nc` * Open the file and load it as an `xarray` dataset: `ds = xr.open_dataset('sst.nc',decode_cf=False)` * examine the data structure: `ds` `Dimensions: (latitude: 261, longitude: 301, time: 12)`\ `Coordinates:`\ `*time (time) float64 1.515e+09 1.517e+09 ... 1.541e+09 1.544e+09`\ `*latitude (latitude) float32 17.025 17.075 17.125 ... 29.975 30.025`\ `*longitude (longitude) float32 195.025 195.075 ... 209.975 210.025`\ `Data variables:` \ `analysed_sst (time, latitude, longitude) float64 ...`\ `Attributes:`\ `acknowledgement: NOAA Coral Reef Watch Program`\ `cdm_data_type: Grid`\ `comment: This product is designed to improve on and re...`\ `contributor_name: NOAA Coral Reef Watch program`\ `contributor_role: Collecting source data and deriving products;...`\ `Conventions: CF-1.6, ACDD-1.3, COARDS`\ `creator_email: coralreefwatch@noaa.gov`\ `creator_institution: NOAA/NESDIS/STAR Coral Reef Watch program`\ `creator_name: NOAA Coral Reef Watch program`\ `creator_type: group`\ `creator_url: https://coralreefwatch.noaa.gov/` \ `data_source: NOAA Daily Global 5km Geo-Polar Blended Night...`\ `date_created: 2018-01-01T00:00:00Z`\ `date_issued: 2018-12-02T15:20:07Z`\ `date_metadata_modified: 2018-09-01T00:00:00Z`\ `date_modified: 2018-01-01T00:00:00Z`\ `Easternmost_Easting: 210.025`\ `geospatial_bounds: "POLYGON((-90.0 360.0, 90.0 360.0, 90.0 0.0, ...`\ `geospatial_bounds_crs: EPSG:32663`\ `geospatial_lat_max: 30.025`\ `geospatial_lat_min: 17.025`\ `geospatial_lat_resolution: 0.049999999999999996`\ `geospatial_lat_units: degrees_north`\ `geospatial_lon_max: 210.025`\ `geospatial_lon_min: 195.025`\ `geospatial_lon_resolution: 0.05000000000000001`\ `geospatial_lon_units: degrees_east`\ `history: Mon Mar 2 06:00:23 2020: ncatted -O -a geosp...`\ `id: CoralTemp-v1.0`\ `infoUrl: https://coralreefwatch.noaa.gov/satellite/ble...`\ `institution: NOAA/NESDIS/STAR Coral Reef Watch program`\ `instrument: ATSR-1, ATSR-2, AATSR, AVHRR, AVHRR-2, AVHRR-...`\ `instrument_vocabulary: NOAA NODC Ocean Archive System Instruments`\ `keywords: 5km, analysed, analysed_sst, analysis, blende...`\ `keywords_vocabulary: GCMD Science Keywords`\ `license: OSTIA Usage Statement (1985-2002): IMPORTANT ...`\ `metadata_link: https://coralreefwatch.noaa.gov/satellite/ble...`\ `naming_authority: gov.noaa.coralreefwatch`\ `NCO: 4.3.7`\ `nco_openmp_thread_number: 1`\ `Northernmost_Northing: 30.025`\ `platform: Ships, drifting buoys, moored buoys, TOGA-TAO...`\ `platform_vocabulary: NOAA NODC Ocean Archive System Platforms`\ `processing_level: L4`\ `product_version: 1.0`\ `program: NOAA Coral Reef Watch program`\ `project: NOAA Coral Reef Watch program`\ `publisher_email: coralreefwatch@noaa.gov`\ `publisher_institution: NOAA/NESDIS/STAR Coral Reef Watch program`\ `publisher_name: NOAA Coral Reef Watch program`\ `publisher_type: group`\ `publisher_url: https://coralreefwatch.noaa.gov/` \ `references: Donlon, et al., 2011. The Operational Sea Sur...`\ `source: OSTIA Sea Surface Temperature Reanalysis (nig...`\ `sourceUrl: (local files)`\ `Southernmost_Northing: 17.025`\ `standard_name_vocabulary: CF Standard Name Table v27`\ `summary: CoralTemp 5km gap-free analysed blended sea s...`\ `time_coverage_duration: P1D`\ `time_coverage_end: 2018-12-01T12:00:00Z`\ `time_coverage_resolution: P1D`\ `time_coverage_start: 2018-01-01T12:00:00Z`\ `title: Sea Surface Temperature, Coral Reef Watch, Co...`\ `Westernmost_Easting: 195.025` * examine which coordinates and variables are included in the dataset: `ds.coords` `Coordinates:`\ `*time (time) float64 1.515e+09 1.517e+09 ... 1.541e+09 1.544e+09`\ `*latitude (latitude) float32 17.025 17.075 17.125 ... 29.925 29.975 30.025`\ `*longitude (longitude) float32 195.025 195.075 195.125 ... 209.975 210.025` `ds.data_vars` `Data variables:` \ `analysed_sst (time, latitude, longitude) float64 ...` * examine the structure of `analysed_sst:` `ds.analysed_sst.shape` `(12, 261, 301)` Our dataset is a 3-D array with 261 rows corresponding to latitudes and 301 columns corresponding to longitudes, for each of the 12 time steps. * get the dates for each time step: `dates=nc.num2date(ds.time,ds.time.units)` \ `dates` `array([datetime.datetime(2018, 1, 1, 12, 0),`\ `datetime.datetime(2018, 2, 1, 12, 0),`\ `datetime.datetime(2018, 3, 1, 12, 0),`\ `datetime.datetime(2018, 4, 1, 12, 0),`\ `datetime.datetime(2018, 5, 1, 12, 0),`\ `datetime.datetime(2018, 6, 1, 12, 0),`\ `datetime.datetime(2018, 7, 1, 12, 0),`\ `datetime.datetime(2018, 8, 1, 12, 0),`\ `datetime.datetime(2018, 9, 1, 12, 0),`\ `datetime.datetime(2018, 10, 1, 12, 0),`\ `datetime.datetime(2018, 11, 1, 12, 0),`\ `datetime.datetime(2018, 12, 1, 12, 0)], dtype=object)` ## Working with the extracted data ### Creating a map for one time step Let's create a map of SST for January 2018 (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')` * set some color breaks `np.nanmin(ds.analysed_sst)` `17.922142857142862` `np.nanmax(ds.analysed_sst)` `28.390645161290323` `levs = np.arange(17.5, 28.5, 0.05)` * define a color palette `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))` * plot the SST map `plt.contourf(ds.longitude, ds.latitude, ds.analysed_sst[0,:,:], levs,cmap=cm)` * plot color scale `plt.colorbar()` * example of how to add points to the map `plt.scatter(range(202,206),np.repeat(26,4),c='black')` * example of how to add a contour line `plt.contour(ds.longitude, ds.latitude, ds.analysed_sst[0,:,:],levels=20,linewidths=1)` * plot title `plt.title("Monthly Sea Surface Temperature " + dates[0].strftime('%b %Y'))`\ `plt.show()` ![](/files/-M6VyJsTB4DbCXQRbiV3) ### Plotting a time series Let's pick the following box : 18-23N, 200-206E. We are going to generate a time series of mean SST within that box. * first, let subset our data: `lat_bnds, lon_bnds = [18, 23], [200, 206]`\ `da=ds.sel(latitude=slice(lat_bnds), longitude=slice(lon_bnds))` * let's plot the subset: ![](/files/-M6Vy7JJO7GF0QIYJTFk) * let's compute the monthly mean over the bounding region: `res=np.mean(da.analysed_sst,axis=(1,2))` * let's plot the time-series: `plt.figure(figsize=(8,4))` \ `plt.scatter(dates,res)` \ `plt.ylabel('SST (ºC)')` ![](/files/-M6Vyil4flw-T8qBxERX) ### Creating a map of average SST over a year * let's compute the yearly mean for the region: `mean_sst=np.mean(ds.analysed_sst,axis=0)` * let's plot the map of the 2018 average SST in the region: `plt.contourf(ds.longitude, ds.latitude, mean_sst, levs,cmap=cm)` \ `plt.colorbar()` \ `plt.title("Mean SST " + dates[0].strftime('%Y-%m')+' - '+dates[11].strftime('%Y-%m'))` \ `plt.show()` ![](/files/-M6VzPOV2zHTqCK7zFTA) --- # Agent Instructions: Querying This Documentation If you need additional information that is not directly available in this page, you can query the documentation dynamically by asking a question. Perform an HTTP GET request on the current page URL with the `ask` query parameter: ``` GET https://coastwatch.gitbook.io/satellite-course/tutorials/python-tutorial/1.-how-to-work-with-satellite-data-in-python.md?ask= ``` The question should be specific, self-contained, and written in natural language. The response will contain a direct answer to the question and relevant excerpts and sources from the documentation. Use this mechanism when the answer is not explicitly present in the current page, you need clarification or additional context, or you want to retrieve related documentation sections.