Polygon data access and visualization

Submitted by pwadmin on Tue, 07/11/2017 - 10:57

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.

In [13]:
# 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).

In [14]:
Multi = MultiPolygon([shape(pol['geometry']) for pol in fiona.open('/home/jenn/aerdData/CCAMLR_Statistical_Areas/statistical_areasPolygon.shp')])

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.

In [15]:
i = 0
for polygon in Multi:
    if i == 0:
        project = partial( 
            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
In [16]:
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:

    # load the data into python as a netcdf_dataset
    data4326 = netcdf_dataset(file)
    print('Data not returned. See examples on generating data queries for help.')
    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.

In [17]:
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.

In [18]:
thisproj = ccrs.SouthPolarStereo() #define projection
fig = plt.figure()    
ax1 = plt.axes(projection = thisproj, ) # set projection

# 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
<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.