Here, we create a map without radar data to concentrate on the other
layers.

In [1]:

importwradlibaswrlimportmatplotlib.pyplotasplimportwarningswarnings.filterwarnings('ignore')try:get_ipython().magic("matplotlib inline")except:pl.ion()importnumpyasnp# Some more matplotlib tools we will need...importmatplotlib.tickerastickerfrommatplotlib.colorsimportLogNormfrommpl_toolkits.axes_grid1importmake_axes_locatable

Note: we organise the code in functions which we can re-use in this
notebook.

In [2]:

defplot_dem(ax):filename=wrl.util.get_wradlib_data_file('geo/bangladesh.tif')ds=wrl.io.open_raster(filename)# pixel_spacing is in output units (lonlat)ds=wrl.georef.reproject_raster_dataset(ds,spacing=0.005)rastervalues,rastercoords,proj=wrl.georef.extract_raster_dataset(ds)# specify kwargs for plotting, using terrain colormap and LogNormdem=ax.pcolormesh(rastercoords[...,0],rastercoords[...,1],rastervalues,cmap=pl.cm.terrain,norm=LogNorm(),vmin=1,vmax=3000)# make some space on the right for colorbar axisdiv1=make_axes_locatable(ax)cax1=div1.append_axes("right",size="5%",pad=0.1)# add colorbar and title# we use LogLocator for colorbarcb=pl.gcf().colorbar(dem,cax=cax1,ticks=ticker.LogLocator(subs=range(10)))cb.set_label('terrain height [m]')

defplot_borders(ax):# country listcountries=['India','Nepal','Bhutan','Myanmar']# open the input data source and get the layerfilename=wrl.util.get_wradlib_data_file('geo/ne_10m_admin_0_boundary_''lines_land.shp')dataset,inLayer=wrl.io.open_vector(filename)# iterate over countries, filter accordingly, get coordinates and plotforitemincountries:# SQL-like selection syntaxfattr="(adm0_left = '"+item+"' or adm0_right = '"+item+"')"inLayer.SetAttributeFilter(fattr)# get borders and namesborders,keys=wrl.georef.get_vector_coordinates(inLayer,key='name')wrl.vis.add_lines(ax,borders,color='black',lw=2,zorder=4)ax.autoscale()

defplot_rivers(ax):# plot rivers from esri vector shape, filter spatially# http://www.fao.org/geonetwork/srv/en/metadata.show?id=37331# open the input data source and get the layerfilename=wrl.util.get_wradlib_data_file('geo/rivers_asia_37331.shp')dataset,inLayer=wrl.io.open_vector(filename)# do spatial filtering to get only geometries inside bounding boxinLayer.SetSpatialFilterRect(88,20,93,27)rivers,keys=wrl.georef.get_vector_coordinates(inLayer,key='MAJ_NAME')# plot on ax1, and ax4wrl.vis.add_lines(ax,rivers,color=pl.cm.terrain(0.),lw=0.5,zorder=3)

The 5 biggest cities of bangladesh are added using simple matplotlib
functions.

In [10]:

defplot_cities(ax):# plot city dots with annotation, finalize plot# lat/lon coordinates of five cities in Bangladeshlats=[23.73,22.32,22.83,24.37,24.90]lons=[90.40,91.82,89.55,88.60,91.87]cities=['Dhaka','Chittagong','Khulna','Rajshahi','Sylhet']forlon,lat,cityinzip(lons,lats,cities):ax.plot(lon,lat,'ro',zorder=5)ax.text(lon+0.01,lat+0.01,city,fontsize='large')

defplot_wgs84(ax):fromosgeoimportosrwgs84=osr.SpatialReference()wgs84.ImportFromEPSG(4326)# some testing on additional axes# add Bangladesh to countriescountries=['India','Nepal','Bhutan','Myanmar','Bangladesh']# create colors for country-patchescm=pl.cm.jetcolors=[]foriinrange(len(countries)):colors.append(cm(1.*i/len(countries)))# open the input data source and get the layerfilename=wrl.util.get_wradlib_data_file('geo/ne_10m_admin_0_''countries.shp')dataset,layer=wrl.io.open_vector(filename)# filter spatially and plot as PatchCollection on ax3layer.SetSpatialFilterRect(88,20,93,27)patches,keys=wrl.georef.get_vector_coordinates(layer,dest_srs=wgs84,key='name')i=0forname,patchinzip(keys,patches):# why comes the US in here?ifnameincountries:wrl.vis.add_patches(ax,patch,facecolor=colors[i],cmap=pl.cm.viridis,alpha=0.4)i+=1ax.autoscale(True)ax.set_aspect('equal')ax.set_xlabel('Longitude')ax.set_ylabel('Latitude')ax.set_title('South Asia - WGS 84')