GoogleEarth
GoogleEarth layer of risk maps
Kmz
Open with GoogleEarth: KMZ
Code¶
import sys
import matplotlib
matplotlib.use('Agg') # Must be before importing matplotlib.pyplot or pylab!
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap,maskoceans
from matplotlib.path import Path
from matplotlib.patches import PathPatch
#from osgeo import gdal
import numpy
import shapefile
#from netCDF4 import Dataset,num2date
from scipy.io import netcdf
import numpy as N
import scipy.ndimage
from scipy.ndimage.filters import gaussian_filter
import matplotlib.mlab as ml
from scipy.interpolate import griddata
from scipy import interpolate
import math
import datetime
from simplekml import (Kml, OverlayXY, ScreenXY, Units, RotationXY,
AltitudeMode, Camera)
#import mpl_toolkits.basemap.pyproj as pyproj
import numpy as np
params = { 'legend.fontsize': 10,\
'figure.subplot.left': 0.05,\
'figure.subplot.right': 0.95,\
'figure.subplot.bottom': 0.05,\
'figure.subplot.top': 0.95,\
'font.family': 'serif'
}
plt.rcParams.update(params)
def deg2num(lat_deg, lon_deg, zoom):
lat_rad = math.radians(lat_deg)
n = 2.0 ** zoom
xtile = int((lon_deg + 180.0) / 360.0 * n)
ytile = int((1.0 - math.log(math.tan(lat_rad) + (1 / math.cos(lat_rad))) / math.pi) / 2.0 * n)
return (xtile, ytile)
def make_kml(llcrnrlon, llcrnrlat, urcrnrlon, urcrnrlat,
figs, timestamps,kmzname,para,colorbar=None, **kw):
"""TODO: LatLon bbox, list of figs, optional colorbar figure,
and several simplekml kw..."""
kml = Kml()
altitude = kw.pop('altitude', 2e7)
roll = kw.pop('roll', 0)
tilt = kw.pop('tilt', 0)
altitudemode = kw.pop('altitudemode', AltitudeMode.relativetoground)
camera = Camera(latitude=np.mean([urcrnrlat, llcrnrlat]),
longitude=np.mean([urcrnrlon, llcrnrlon]),
altitude=altitude, roll=roll, tilt=tilt,
altitudemode=altitudemode)
kml.document.camera = camera
draworder = 0
for fig in figs: # NOTE: Overlays are limited to the same bbox.
draworder += 1
ground = kml.newgroundoverlay(name=timestamps[draworder-1])
ground.draworder = draworder
ground.visibility = kw.pop('visibility', 2)
ground.name = kw.pop('name',str(timestamps[draworder-1]))
ground.color = kw.pop('color', '9effffff')
ground.atomauthor = kw.pop('author', 'P. Hoffmann')
ground.latlonbox.rotation = kw.pop('rotation', 0)
ground.description = kw.pop('description', 'Matplotlib figure')
ground.gxaltitudemode = kw.pop('gxaltitudemode','clampToSeaFloor')
ground.icon.href = fig
ground.latlonbox.east = llcrnrlon
ground.latlonbox.south = llcrnrlat
ground.latlonbox.north = urcrnrlat
ground.latlonbox.west = urcrnrlon
ta = str('%4i-%02i-%02i'%(timestamps[draworder-1],1,1))
te = str('%4i-%02i-%02i'%(timestamps[draworder-1]+9,12,31))
ground.timespan.begin = ta
ground.timespan.end = te
# pnt = kml.newpoint(name="Oberwiesenthal", coords=[(12.9833,50.4167)])
# pnt.description = '<img src="%s_Oberwiesenthal.png" alt="picture" width="600" height="400" align="left"/>'%para
# pnt.style.iconstyle.icon.href = 'http://maps.google.com/mapfiles/kml/shapes/mountains.png'
# pnt = kml.newpoint(name="Oberhof", coords=[(10.7333,50.7167)])
# pnt.description = '<img src="%s_Oberhof.png" alt="picture" width="600" height="400" align="left"/>'%para
# pnt.style.iconstyle.icon.href = 'http://maps.google.com/mapfiles/kml/shapes/mountains.png'
# pnt = kml.newpoint(name="Brocken", coords=[(10.615559,51.799087)])
# pnt.description = '<img src="%s_Brocken.png" alt="picture" width="600" height="400" align="left"/>'%para
# pnt.style.iconstyle.icon.href = 'http://maps.google.com/mapfiles/kml/shapes/mountains.png'
pnt = kml.newpoint(name="Potsdam", coords=[(13.0,52.5)])
pnt.description = '<img src="timeseries_Potsdam.png" alt="picture" width="600" height="400" align="left"/>'
pnt.style.iconstyle.icon.href = 'http://maps.google.com/mapfiles/kml/pal4/icon48.png'
if colorbar: # Options for colorbar are hard-coded (to avoid a big mess).
screen = kml.newscreenoverlay(name='ScreenOverlay')
screen.icon.href = colorbar
screen.overlayxy = OverlayXY(x=0, y=0,
xunits=Units.fraction,
yunits=Units.fraction)
screen.screenxy = ScreenXY(x=0.015, y=0.075,
xunits=Units.fraction,
yunits=Units.fraction)
screen.rotationXY = RotationXY(x=0.5, y=0.5,
xunits=Units.fraction,
yunits=Units.fraction)
screen.size.x = 0
screen.size.y = 0
screen.size.xunits = Units.fraction
screen.size.yunits = Units.fraction
screen.visibility = 1
kmzfile = kw.pop('kmzfile',kmzname)
kml.savekmz(kmzfile)
def gearth_fig(llcrnrlon, llcrnrlat, urcrnrlon, urcrnrlat, pixels=1024):
"""Return a Matplotlib `fig` and `ax` handles for a Google-Earth Image."""
aspect = np.cos(np.mean([llcrnrlat, urcrnrlat]) * np.pi/180.0)
xsize = np.ptp([urcrnrlon, llcrnrlon]) * aspect
ysize = np.ptp([urcrnrlat, llcrnrlat])
aspect = ysize / xsize
if aspect > 1.0:
figsize = (10.0 / aspect, 10.0)
else:
figsize = (10.0, 10.0 * aspect)
if False:
plt.ioff() # Make `True` to prevent the KML components from poping-up.
fig = plt.figure(figsize=figsize,
frameon=False,
dpi=pixels//10)
# KML friendly image. If using basemap try: `fix_aspect=False`.
ax = fig.add_axes([0, 0, 1, 1])
ax.set_xlim(llcrnrlon, urcrnrlon)
ax.set_ylim(llcrnrlat, urcrnrlat)
return fig, ax
################################
para = 'prc_ext_ja'
figs = []
timestamps = []
for p in ['1961-1990','1990-2019']:
nc = netcdf.netcdf_file('../csv/gen_%s/dw/prc_ext_ja.cdf'%p)
fac = nc.variables['pr'].scale_factor
add = nc.variables['pr'].add_offset
print (add,fac)
lon = N.array(nc.variables['longitude'][:])
lat = N.array(nc.variables['latitude'][:])
dat = N.array(nc.variables['pr'][:])*fac+add
dat = dat*24.*3600.
nc.close()
print (dat.shape,dat.max())
tmp = datetime.date(int(p[0:4]),1,1)#datetime.datetime.strptime('01/01/%s'%jj[d], "%d/%m/%Y").timestamp()
timestamps.append(int(p[0:4]))
fig, ax = gearth_fig(lon.min(),lat.min(),lon.max(),lat.max(), pixels=1024)
m = Basemap(projection='merc',llcrnrlat=lat.min(),urcrnrlat=lat.max(),llcrnrlon=lon.min(),urcrnrlon=lon.max(),resolution='h',epsg=4326)
lons,lats = N.meshgrid(lon,lat)
dats = maskoceans(lons,lats,dat[0,:,:])
xx,yy = m(lons,lats)
cs = m.contourf(xx,yy,dats,levels=[50,100,150,200,300],cmap=plt.get_cmap('Purples'),extend='max',alpha=1)
m.contour(xx,yy,dats,levels=[50,100,150,200,300],colors='k')
x,y = m(-10,45)
plt.text(x,y,'%s'%p,fontsize=16,weight='bold',color='w')
plt.axis('off')
fname = '%s.png'%(p)
figs.append(fname)
plt.savefig(fname,dpi=240,transparent=True,bbox_inches='tight',pad_inches=0.0)
plt.close()
nc.close()
fig = plt.figure(figsize=(1.0, 4.0), facecolor=None, frameon=False)
ax = fig.add_axes([0.0, 0.05, 0.2, 0.9])
cb = fig.colorbar(cs, cax=ax)
#cb.set_label('colorbar label', color=fg_color)
cb.ax.tick_params(labelsize=14)
cb.ax.yaxis.set_tick_params(color='w')
plt.setp(plt.getp(cb.ax.axes, 'yticklabels'), color='w',weight='bold')
#cb.ax.set_xticklabels([5,10,15,20,25,30],colors='w')
cb.ax.xaxis.set_tick_params(color='w')
cb.set_label('5-Days Extreme Precipitation (mm)', rotation=90,color='w',weight='bold',fontsize=12,labelpad=0.5)
fig.savefig('legend.png', transparent=False, format='png')
make_kml(lon.min(),lat.min(),lon.max(),lat.max(),figs=figs,timestamps=timestamps,kmzname='earth.kmz',para=para,colorbar='legend.png')