lmes_64_wbackground_SST_avg_2025-2050

How to plot CMIP5 SST data within polygons of large marine ecosystems

I have started using the next IPCC AR5 model simulations that have been made available online. As I described in previous post I have managed to download and organise all of these data so that I now have one file containing the ensemble average SST climatology of 11 ESM models for the period 1961-1990 as well as one file containing the ensemble average projections for 2006-2025 and 2025-2050 for the RCP8.5 simulations. My problem was that I needed to calculate the average SST within each polygon defining a large marine ecosystems (LME) . The polygons for the LMEs were available as GIS Shapefiles which I had never worked with before. I struggled for a long time to get this to work, however I got some good help through Stackoverflow.com .

I used the OGR module or Python to read the Shapefile and extract each polygon in each feature layer contained in the Shapefile. In most cases, the LMEs consisted of several polygons (multi polygons). The basic method I used to read the points defining each polygon was:

if (geom.GetGeometryType() == ogr.wkbMultiPolygon):
for i in range(geom.GetGeometryCount()):
r = geom.GetGeometryRef(i)
for part in r:
x = [part.GetX(j) for j in range(part.GetPointCount())]
y = [part.GetY(j) for j in range(part.GetPointCount())]

codes += [mpath.Path.MOVETO] + (len(x)-1)*[mpath.Path.LINETO]
all_x += x
all_y += y

Since I wanted to plot the SST within the polygons on a map, I used Basemap to keep track of the map coordinates.

mymap = Basemap(resolution='l',projection='robin',lon_0=0)

This map object can then be used to convert the longitude-latitude position of each polygon to map coordinates, which makes it much easier to work with maps on a global scale.

pathX,pathY=mymap(all_x,all_y)

I stored my SST anomalies in numpy arrays that I flattened to make it easy to extract only the SST grid points found within each polygon so that I could calculate the average SST anomaly. To find which points are inside a polygon I used the brilliant method inside_poly from matplotlib

insideBoolean = plt.mlab.inside_poly(np.c_[mymapgrid[:,0],mymapgrid[:,1]],np.c_[X,Y])

Each polygon is stored as a path object in map coordinates and can be defined and stored as a patch in a list. I also store a list of SST values (one for each polygon) to be used as facecolor for the polygon:

mypath = mpath.Path(np.column_stack((all_x,all_y)), codes)
path_patch = patches.PathPatch(mypath, lw=1)
mypatches.append(path_patch)
myAVG=np.ma.mean(myaveragesanom)
mycolors.append(myAVG)

Finally, all of the polygons are plotted simultaneously using PatchCollection:

p = PatchCollection(mypatches, cmap=plt.get_cmap('RdYlBu_r'), alpha=1.0)
p.set_array(array(mycolors))

The result is seen above and the full script used to create the map with polygons for LMEs is available for download.

import numpy as np
import math
from pylab import *
import matplotlib.patches as patches
from matplotlib.collections import PatchCollection
import string, os, sys
import datetime, types
from netCDF4 import Dataset
import matplotlib.nxutils as nx
from mpl_toolkits.basemap import Basemap, shiftgrid
import ogr
import matplotlib.path as mpath
import matplotlib.patches as patches
import datetime as datetime

__author__ = 'Trond Kristiansen'
__email__ = 'me (at) trondkristiansen.com'
__created__ = datetime.datetime(2013, 5, 10)
__modified__ = datetime.datetime(2013, 5, 10)
__version__ = "1.1.5"
__status__ = "Production"

"""This script reads a Shapefile containing polygons that edfine the boundaries
of large marine ecosystems (LMEs). Each polygon is read using the OGR module
and extracted point by point and then re-assembled as a polygon as described here:

http://stackoverflow.com/questions/16326068/create-closed-polygon-from-boundary-points

For each polygon, the number of points from a global sea surface temperature
netcdf file is estimated .The grid of the global file is assumed to be a
rectangular grid of size 360x180, but this could be arbitraruy as long as the
longitude latitude data is available. Also note that this program shifts the grid from 0-360 to
-180 to 180 using the shiftgrid function of basemap. If the input grid is already in this
format then remove the calls to shiftgrid.

The output is a CSV file containing the average SST values for each polygon. The output
file contains the names and region codes of LMEs or FAO regions depending on what sort
of input Shapefile you are using.
"""


def getLMEpolygon(coordinatefile,mymap,index,first,typename,filetype):

ds = ogr.Open(coordinatefile)
lyr = ds.GetLayer(0)
numberOfPolygons=lyr.GetFeatureCount()

if first is False:
ft = lyr.GetFeature(index)
geom = ft.GetGeometryRef()
polygonName=ft.items()[typename]

if filetype=="LME":
polygonNumber=ft.items()["LME_NUMBER"]
polygonCode=ft.items()["LME_code"]

if filetype=="FAO":
polygonNumber=ft.items()["F_AREA"]
polygonCode=ft.items()["SUBOCEAN"]

"""Create lists to store values for polygons"""
codes = []; all_x = []; all_y = []; all_XY= []

"""Check if polygon:"""
if (geom.GetGeometryType() == ogr.wkbPolygon):
#print "Found polygon:", polygonName
for i in xrange(geom.GetGeometryCount()):

r = geom.GetGeometryRef(i)
x = [r.GetX(j) for j in xrange(r.GetPointCount())]
y = [r.GetY(j) for j in xrange(r.GetPointCount())]

codes += [mpath.Path.MOVETO] + (len(x)-1)*[mpath.Path.LINETO]
all_x += x
all_y += y
all_XY +=mymap(x,y)

"""Check if multipolygon:"""
if (geom.GetGeometryType() == ogr.wkbMultiPolygon):
#print "Geometric features", geom.GetGeometryCount()
for i in range(geom.GetGeometryCount()):
r = geom.GetGeometryRef(i)
for part in r:
x = [part.GetX(j) for j in range(part.GetPointCount())]
y = [part.GetY(j) for j in range(part.GetPointCount())]

codes += [mpath.Path.MOVETO] + (len(x)-1)*[mpath.Path.LINETO]
all_x += x
all_y += y

if len(all_x)==0:
all_XY=None
mypoly=None
else:
pathX,pathY=mymap(all_x,all_y)
mypath = mpath.Path(np.column_stack((all_x,all_y)), codes)
mymappath = mpath.Path(np.column_stack((pathX,pathY)), codes)
mypoly = mpath.Path.to_polygons(mypath)
mymappoly = mpath.Path.to_polygons(mymappath)
mypolysX=[];mypolysY=[]

for poly in xrange(geom.GetGeometryCount()):
myarray=np.squeeze(mypoly[:][poly])
myX,myY=mymap(myarray[:,0],myarray[:,1])

mypolysX.append(np.squeeze(myX))
mypolysY.append(np.squeeze(myY))

if first is True or polygonName in ["West Bering Sea", "Red Sea"]:
print "Will extract data for %s polygons"%(numberOfPolygons)
mypolysX=None
mypolysY=None
mypoly=None
mypath=None
mymappath=None
mymappoly=None
polygonName=None
polygonCode=None
polygonNumber=None
first=False

return mypolysX, mypolysY, mypoly, mymappath, mymappoly, first, numberOfPolygons, polygonName, polygonNumber, polygonCode

"""Function that opens a CMIP5 file and reads the contents. The innput
is assumed to be on grid 0-360 so all values are shifted to new grid
on format -180 to 180 using the shiftgrid function of basemap."""

def openCMIP5file(CMIP5name,myvar,mymap):
if os.path.exists(CMIP5name):
myfile=Dataset(CMIP5name)
print "Opened CMIP5 file: %s"%(CMIP5name)
else:
print "Could not find CMIP5 input file %s : abort"%(CMIP5name)
sys.exit()
myTEMP=np.squeeze(myfile.variables[myvar][0,:,:])
myTEMP=myTEMP-273.15

lonCMIP5=np.squeeze(myfile.variables["lon"][:])
latCMIP5=np.squeeze(myfile.variables["lat"][:])

myTEMPin, lonsin = shiftgrid(180, myTEMP, lonCMIP5, start=False)
lons,lats=np.meshgrid(lonsin,latCMIP5)

lons=lons.flatten()
lats=lats.flatten()
mygrid=np.empty((len(lats),2))
mymapgrid=np.empty((len(lats),2))

for i in xrange(len(lats)):
mygrid[i,0]=lons[i]
mygrid[i,1]=lats[i]
X,Y=mymap(lons[i],lats[i])
mymapgrid[i,0]=X
mymapgrid[i,1]=Y

return myTEMPin, mygrid, mymapgrid,lonsin,latCMIP5

"""Map related functions:"""
def drawMap(NUM_COLORS):

plt.figure(figsize=(12,10), frameon=False)
ax = plt.subplot(111)
cm = plt.get_cmap('RdBu')
ax.set_color_cycle([cm(1.*j/NUM_COLORS) for j in range(NUM_COLORS)])

mymap = Basemap(resolution='l',projection='robin',lon_0=0)

return ax, mymap, cm

def finishMap(mymap):
mymap.drawcountries()
mymap.drawcoastlines()
mymap.fillcontinents(color='lightgrey',lake_color='white')
mymap.drawmapboundary(fill_color='white')

""" --------------------"""

"""Edit the correct names below:"""

LMEcoordinatefiles=[]
LMEcoordinatefiles.append('ShapefileBoundaries/lmes_64.shp')
LMEcoordinatefiles.append('ShapefileBoundaries/LME_300.shp')
LMEcoordinatefiles.append('ShapefileBoundaries/LME_0_300.shp')
LMEcoordinatefiles.append('ShapefileBoundaries/FAO_total_clip.shp')

start_year=2006
end_year=2025

workDir="/Users/trond/Projects/ESM-PP/SST/RCP85/MODELS/REGRID/YRLYAVG/"
workDirHist="/Users/trond/Projects/ESM-PP/SST/Historical/MODELS/REGRID/YRLYAVG/"

CMIP5file=workDir+"ensemble_average_"+str(start_year)+"_"+str(end_year)+".nc"
CMIP5HISTfile=workDirHist+"ensemble_average_1961_1990.nc"

mydebug=False
doPoints=False
doSSTbackground=True

for LMEcoordinatefile in LMEcoordinatefiles:
first=True
print "Running: %s"%(LMEcoordinatefile)

"""Save filenames and paths for creating output filenames"""
filetype=(os.path.basename(LMEcoordinatefile)[0:3]).upper()
if filetype=="FAO": typename="OCEAN"
else: typename="LME_NAME"
mainfilename=os.path.basename(LMEcoordinatefile)[0:-4]

"""Open output file where stats are written"""
outputfile="results/"+mainfilename+"_anomaliesSST_year-"+str(start_year)+"-"+str(end_year)+".csv"
if os.path.exists(outputfile): os.remove(outputfile)
myoutfile=open(outputfile,'a')
print "Results will be stored into file: %s"%(outputfile)

"""initialize the map:"""
mymap=None; mypatches=[]; mycolors=[]
mypolysX, mypolysY, mypoly, mypath, mymappoly,first, numberOfPolygons, polygonName, polygonNumber, polygonCode = getLMEpolygon(LMEcoordinatefile, mymap, 0,first,typename,filetype)
NUM_COLORS=numberOfPolygons
ax, mymap, cm = drawMap(NUM_COLORS)

"""Get the CMIP5 data together with the grid"""
SST,mygrid, mymapgrid,lonCMIP5,latCMIP5 = openCMIP5file(CMIP5file,"tos",mymap)
SSThist,mygrid, mymapgrid,lonCMIP5,latCMIP5 = openCMIP5file(CMIP5HISTfile,"tos",mymap)

"""Calculate the SST anomalies"""
SSTanom=np.subtract(SST,SSThist)

print "Anomalies mean: %s range: %s to %s "%(SSTanom.mean(), SSTanom.min(), SSTanom.max())
print "Projections SST mean: %s range: %s to %s "%(SST.mean(), SST.min(), SST.max())
print "Historical SST mean: %s range: %s to %s "%(SSThist.mean(), SSThist.min(), SSThist.max())

"""For each LME of interest create a polygon of coordinates defining the boundaries"""
for counter in xrange(numberOfPolygons-1):

mypolysX, mypolysY, mypoly, mypath, mymappoly,first,numberOfPolygons, polygonName, polygonNumber, polygonCode = getLMEpolygon(LMEcoordinatefile, mymap,counter,first,typename,filetype)

if mypolysX != None:
"""Find the indices inside the grid that are within the polygon"""
SSTanomflat=SSTanom.flatten()
SSTflat=SST.flatten()

myaveragesanom=[]; myaverages=[]; CMIP5points=[]
for myc in xrange(len(mypolysX)):
X= np.asarray(mypolysX[:][myc])
Y= np.asarray(mypolysY[:][myc])

"""Find CMIP5 grid points inside the polygon in question"""
insideBoolean = plt.mlab.inside_poly(np.c_[mymapgrid[:,0],mymapgrid[:,1]],np.c_[X,Y])

"""Calculate the average SST value of points inside the polygon"""
if len(insideBoolean) > 0:
CMIP5points.append(len(insideBoolean))
mymapgrid=np.c_[mymapgrid[:,0],mymapgrid[:,1]]

myaverageSSTanom=np.ma.mean(SSTanomflat[insideBoolean])
print "Polygon: %s average SST: %s"%(polygonName,myaverageSSTanom)
myaverageSST=np.ma.mean(SSTflat[insideBoolean])

myaverages.append(myaverageSST)
myaveragesanom.append(myaverageSSTanom)

if doPoints is True:

for point in xrange(len(insideBoolean)):
pointX=mymapgrid[insideBoolean[point],0]
pointY=mymapgrid[insideBoolean[point],1]
ax.scatter(pointX,pointY,2,color='red')
ax.hold(True)

"""Write the header for the file"""
if counter==0:
if filetype=="FAO":
myheader="Shapefile, polygon name, F_AREA, SUBOCEAN, CMIP5 points inside polygon, average SST, anomaly SST\n"
if filetype=="LME":
myheader="Shapefile, polygon name, LME_NUMBER, LME_CODE, CMIP5 points inside polygon, average SST, anomaly SST\n"

myoutfile.writelines(myheader)

mydata="%s, %40s, %3i, %3i, %3i, %3.3f, %3.3f\n"%(mainfilename,polygonName,int(polygonNumber),int(polygonCode),np.ma.sum(CMIP5points), np.ma.mean(myaverages), np.ma.mean(myaveragesanom))
myoutfile.writelines(mydata)

"""Add the path to a collection of paths to be drawn later. Also save the average value of SST
that will be used to determine the color of the polygon. See http://stackoverflow.com/questions/6028675/setting-color-range-in-matplotlib-patchcollection"""

myAVG=np.ma.mean(myaveragesanom)
mycolors.append(myAVG)
path_patch = patches.PathPatch(mypath, lw=1)
mypatches.append(path_patch)

"""Plot the background SST as grey shaded low alpha values"""
if doSSTbackground is True and filetype=="LME":
lons,lats=np.meshgrid(lonCMIP5,latCMIP5)
x,y=mymap(lons,lats)
levels = np.arange(0.2,2.5, 0.25)
mymap.contourf(x,y,SSTanom,levels,cmap=plt.get_cmap('Greys', len(levels)-1),alpha=0.5)

"""Plot the polygons"""
p = PatchCollection(mypatches, cmap=plt.get_cmap('RdYlBu_r'), alpha=1.0)
p.set_array(array(mycolors))

"""The polygon patches requires me to set the clim value unless all polygons are one color"""
p.set_clim([np.ma.min(mycolors),np.ma.max(mycolors)])
print "Range of colorbar is :%s to %s"%(np.ma.min(mycolors),np.ma.max(mycolors))
p.set_clim([0.2,2.5])

plt.colorbar(p,shrink=0.5)
if doPoints is False:
ax.add_collection(p)

"""Add some finishing touches to the map"""
finishMap(mymap)

plotname="figures/"+mainfilename+"_SST_avg_"+str(start_year)+"-"+str(end_year)+".png"
if doPoints is True:
plotname="figures/"+mainfilename+"_wpoints_SST_avg_"+str(start_year)+"-"+str(end_year)+".png"
if doSSTbackground is True and filetype=="LME":
plotname="figures/"+mainfilename+"_wbackground_SST_avg_"+str(start_year)+"-"+str(end_year)+".png"
plt.savefig(plotname,dpi=300)
print "Saved plot to file %s"%(plotname)
# plt.show()

clf()

Leave a Reply

Your email address will not be published. Required fields are marked *