print "import modules...", import sys sys.stdout.flush() import pickle from mpl_toolkits.basemap import Basemap, shiftgrid, cm import matplotlib import matplotlib.pyplot as plt import numpy as np from netCDF4 import Dataset print "ok" # Lovecraft: 47:9'S 126:43'W lovecraft_lat = -47.9 lovecraft_lon = -126.43 # August Derleth: 49:51'S 128:34'W derleth_lat = -49.51 derleth_lon = -128.34 # Nemo point: 48:52.6'S 123:23.6'W nemo_lat = -48.526 nemo_lon = -123.236 # The Bloop: # Appears to be way too far from the Nemo point to be interesting in a R'lyeh context # bransfield_strait_lat=-63 # bransfield_strait_lon=-59 # ross_sea_lat = -75 # ross_sea_lon = -175 # cape_adare_lat = -71.17 # cape_adare_lon = -170.14 # Not necessary, because the default projection is ortho, # but can be useful if you want another one. def equi(m, centerlon, centerlat, radius, *args, **kwargs): """ Drawing circles of a given radius around any point on earth, given the current projection. http://www.geophysique.be/2011/02/20/matplotlib-basemap-tutorial-09-drawing-circles/ """ glon1 = centerlon glat1 = centerlat X = [] Y = [] for azimuth in range(0, 360): glon2, glat2, baz = shoot(glon1, glat1, azimuth, radius) X.append(glon2) Y.append(glat2) X.append(X[0]) Y.append(Y[0]) #m.plot(X,Y,**kwargs) #Should work, but doesn't... X,Y = m(X,Y) plt.plot(X,Y,**kwargs) def shoot(lon, lat, azimuth, maxdist=None): """Shooter Function Plotting great circles with Basemap, but knowing only the longitude, latitude, the azimuth and a distance. Only the origin point is known. Original javascript on http://williams.best.vwh.net/gccalc.htm Translated to python by Thomas Lecocq : http://www.geophysique.be/2011/02/19/matplotlib-basemap-tutorial-08-shooting-great-circles/ """ glat1 = lat * np.pi / 180. glon1 = lon * np.pi / 180. s = maxdist / 1.852 faz = azimuth * np.pi / 180. EPS= 0.00000000005 if ((np.abs(np.cos(glat1))<EPS) and not (np.abs(np.sin(faz))<EPS)): alert("Only N-S courses are meaningful, starting at a pole!") a=6378.13/1.852 f=1/298.257223563 r = 1 - f tu = r * np.tan(glat1) sf = np.sin(faz) cf = np.cos(faz) if (cf==0): b=0. else: b=2. * np.arctan2 (tu, cf) cu = 1. / np.sqrt(1 + tu * tu) su = tu * cu sa = cu * sf c2a = 1 - sa * sa x = 1. + np.sqrt(1. + c2a * (1. / (r * r) - 1.)) x = (x - 2.) / x c = 1. - x c = (x * x / 4. + 1.) / c d = (0.375 * x * x - 1.) * x tu = s / (r * a * c) y = tu c = y + 1 while (np.abs (y - c) > EPS): sy = np.sin(y) cy = np.cos(y) cz = np.cos(b + y) e = 2. * cz * cz - 1. c = y x = e * cy y = e + e - 1. y = (((sy * sy * 4. - 3.) * y * cz * d / 6. + x) * d / 4. - cz) * sy * d + tu b = cu * cy * cf - su * sy c = r * np.sqrt(sa * sa + b * b) d = su * cy + cu * sy * cf glat2 = (np.arctan2(d, c) + np.pi) % (2*np.pi) - np.pi c = cu * cy - su * sy * cf x = np.arctan2(sy * sf, c) c = ((-3. * c2a + 4.) * f + 4.) * c2a * f / 16. d = ((e * cy * c + cz) * sy * c + y) * sa glon2 = ((glon1 + x - (1. - c) * d * f + np.pi) % (2*np.pi)) - np.pi baz = (np.arctan2(sa, b) + np.pi) % (2 * np.pi) glon2 *= 180./np.pi glat2 *= 180./np.pi baz *= 180./np.pi return (glon2, glat2, baz) print "read in etopo5 topography/bathymetry" url = 'http://ferret.pmel.noaa.gov/thredds/dodsC/data/PMEL/etopo5.nc' etopodata = Dataset(url) print "get data" def topopickle(etopodata,name): import sys print "\t"+name+"...", sys.stdout.flush() filename = "rlyeh_"+name+".pickle" try: with open(filename,"r") as fd: print "load...", var = pickle.load(fd) except IOError: print "copy...", var = etopodata.variables[name][:] with open(filename,"w") as fd: print "dump...", pickle.dump(var,fd) print "ok" return var topoin = topopickle(etopodata,"ROSE") lons = topopickle(etopodata,"ETOPO05_X") lats = topopickle(etopodata,"ETOPO05_Y") print "shift data so lons go from -180 to 180 instead of 20 to 380...", sys.stdout.flush() topoin,lons = shiftgrid(180.,topoin,lons,start=False) print "ok" # create the figure and axes instances. fig = plt.figure() ax = fig.add_axes([0.1,0.1,0.8,0.8]) print "setup basemap" # set up orthographic m projection with # perspective of satellite looking down at 50N, 100W. # use low resolution coastlines. m = Basemap(projection='ortho',lat_0=nemo_lat,lon_0=nemo_lon,resolution='l') m.bluemarble() # Generic Mapping Tools colormaps: # GMT_drywet, GMT_gebco, GMT_globe, GMT_haxby GMT_no_green, GMT_ocean, GMT_polar, # GMT_red2green, GMT_relief, GMT_split, GMT_wysiwyg print "transform to nx x ny regularly spaced native projection grid" # step=5000. step=10000. nx = int((m.xmax-m.xmin)/step)+1; ny = int((m.ymax-m.ymin)/step)+1 topodat = m.transform_scalar(topoin,lons,lats,nx,ny) print "plot topography/bathymetry as shadows" from matplotlib.colors import LightSource ls = LightSource(azdeg = 45, altdeg = 220, hsv_min_val=0.0, hsv_max_val=1.0, hsv_min_sat=0.0, hsv_max_sat=1.0) # convert data to rgb array including shading from light source. # (must specify color m) rgb = ls.shade(topodat, cm.GMT_ocean) im = m.imshow(rgb, alpha=0.15) print "draw coastlines, country boundaries, fill continents" m.drawcoastlines(linewidth=0.25) # draw the edge of the map projection region m.drawmapboundary(fill_color='white') # draw lat/lon grid lines every 30 degrees. m.drawmeridians(np.arange( 0,360,30), color="black" ) m.drawparallels(np.arange(-90,90 ,30), color="black" ) print "draw points" psize=5 font = {'family' : 'serif', 'weight' : 'normal', 'size' : 18} matplotlib.rc('font', **font) x,y = m( lovecraft_lon, lovecraft_lat ) m.scatter(x,y,psize,marker='o', color='white') plt.text(x+50000,y+50000+50000, "Lovecraft", color='white') x,y = m( derleth_lon, derleth_lat ) m.scatter(x,y,psize,marker='o',color='white') plt.text(x+50000-120000,y+50000, "Derleth", color='white', horizontalalignment="right") x,y = m( nemo_lon, nemo_lat ) m.scatter(x,y,psize*3,marker='+',color='#555555') plt.text(x+50000+50000,y+50000-80000, "Nemo", color="#555555", verticalalignment="top") equi(m, nemo_lon, nemo_lat, radius=2688, color="#555555" ) # x,y = m( bransfield_strait_lon, bransfield_strait_lat ) # m.scatter(x,y,psize*3,marker='+',color='#555555') # plt.text(x+50000+20000,y+50000-80000, "bransfield_strait", color="#555555", verticalalignment="baseline") # x,y = m( ross_sea_lon, ross_sea_lat ) # m.scatter(x,y,psize*3,marker='+',color='#555555') # plt.text(x+50000+20000,y+50000-80000, "ross_sea", color="#555555", verticalalignment="baseline") # x,y = m( cape_adare_lon, cape_adare_lat ) # m.scatter(x,y,psize*3,marker='+',color='#555555') # plt.text(x+50000+20000,y+50000-80000, "cape_adare", color="#555555", verticalalignment="baseline") plt.savefig("R'lyeh_locations.png", dpi=600, bbox_inches='tight') # plt.show()