This example is a python notebook, you can run and modify this example directly on your computer, the code is available for download here.
In this example we will download and plot Antarctic satellite temperature data for a polygon of interest (a CCAMLR statistical area). We will load a shapefile of multiple polygons, pull out one polygon, extract corresponding satellite data from an ERDDAP server and plot the returned data in polar projection.
This example runs in python 3 with the following geospatial packages: fiona which opens and writes files, shapely to work with polygons, pyproj to convert between projections and cartopy for map based plotting.
Additional examples in this series will build on this one and will include how to genarate the dataset request url, how to add an ice mask, how to calculate statistics on the temperature data and how to output data in geojson format for leaflet.
# Modules for importing and working with shape files import fiona from shapely.geometry import shape from shapely.geometry import MultiPolygon from functools import partial import pyproj from shapely.ops import transform # Modules for fetching and working with the satellite data import urllib3 from netCDF4 import Dataset as netcdf_dataset import numpy as np # modules for masking and plotting the data from matplotlib import path import matplotlib.pyplot as plt import cartopy import cartopy.crs as ccrs
Now we will use fiona to open a Shapefile that contains multipolygons and read the coverage area geometries. This shapefile is CCAMLR statistical areas and was downloaded from the CCAMLR GIS webpage (https://gis.ccamlr.org/home).
Multi = MultiPolygon([shape(pol['geometry']) for pol in fiona.open('/home/jenn/aerdData/CCAMLR_Statistical_Areas/statistical_areasPolygon.shp')]) #Multi.wkt
For this example we will be using the first polygon from the multipolygon as an example. We will uses pyproj and shapely to reproject the shapefile polygon into EPSG:4326 and turn it into a shapely polygon. The projections of the shapefile data, satellite data and our plot are all different, this step transforms the shapefile data into the 4326 projection of the satellite data. Note that requests to the satellite data server need to be in EPSG:4326.
i = 0 for polygon in Multi: if i == 0: project = partial( pyproj.transform, pyproj.Proj(init='esri:102020'), # source coordinate system pyproj.Proj(init='epsg:4326')) # destination coordinate system p4326 = transform(project, polygon) # new shapely polygon in new projection i = i + 1
dataRequestUrl = 'http://coastwatch.pfeg.noaa.gov/erddap/griddap/jplMURSST41mday.nc?sst[(2017-05-16T00:00:00Z):40:(2017-05-16T00:00:00Z)][(-75.428900):40:(-60.000000)][(-105.000000):40:(-68.123990)]' file = 'dataset.nc' # the filename, local location for the satellite data download http = urllib3.PoolManager() r = http.request('GET', dataRequestUrl, preload_content=False) # Open the file and write the data to it with open(file, 'wb') as out: while True: data = r.read(1024*1024) if not data: break out.write(data) r.release_conn() try: # load the data into python as a netcdf_dataset data4326 = netcdf_dataset(file) except: print('Data not returned. See examples on generating data queries for help.') else: print ('\nSatellite Data File Retrieved') # Pull out temperature data parameter = 'sst' if len(np.shape(data4326.variables[parameter])) == 3: data = data4326.variables[parameter][0, :, :] elif len(np.shape(data4326.variables[parameter])) == 4: data = data4326.variables[parameter][0, 0,:, :] lats = data4326.variables['latitude'][:] lons = data4326.variables['longitude'][:] print('\nQuery was for '+ str(len(lats)*len(lons)) + ' data points')
Satellite Data File Retrieved Query was for 3627 data points
Next we will create a mask using the CCAMLR polygon. This will tell us which data is inside or outside the polygon and we will use it for plotting the data in the next step. We get the perimiter of the CCAMLR polygon, use it as a path for the mask, and reformat the satellite temperature data into a masked array.
polyListx, polyListy = p4326.exterior.xy # perimeter of polygon polyList = list(zip(list(polyListx),list(polyListy))) # formatted perimeter p = path.Path(polyList) # path for mask X, Y = np.meshgrid(lons, lats) # create the grid points = np.array((X.flatten(), Y.flatten())).T # break it down mask = p.contains_points(points).reshape(X.shape) # calc and grid a mask based on polygon datamasked = np.ma.masked_where(mask != True, data) # create masked data array
Now we will make a plot of the temperature data in polar stereographic projection. We will use a colormesh plot and cartopy to create a simple plot of the data the is only within the CCAMLR polygon.
thisproj = ccrs.SouthPolarStereo() #define projection fig = plt.figure() ax1 = plt.axes(projection = thisproj, ) # set projection ax1.set_global() # set type of plot as colormesh, indicate it's in 4326 (Plate Carree) and needs to be projected dataplt = ax1.pcolormesh(X,Y,datamasked,transform=ccrs.PlateCarree()) ax1.set_extent([-4000000, 2500000, -2500000, 2500000], thisproj) # reduce extent/bounds of plot ax1.add_feature( cartopy.feature.LAND, zorder=1, edgecolor='none', facecolor='#fae5c9') # add Land ax1.coastlines(color='#f0c7ab') # add coastline ax1.gridlines(alpha='0.3') # add lat lon rings fig.set_dpi(200) plt.show() plt.clf()
<matplotlib.figure.Figure at 0x7fadd0474b70>
Stay tuned for more examples that will show how to generate the dataset request url, show how to add an ice mask, how to calculate statistics on the temperature data and how to output data in geojson format for leaflet.