Simple Basemap Example

You can do a lot with Python basemaps. Here is an example of a few things you can do...

Plotting maps with a loop: Best Practices

Best practices for plotting data on a map with a loop...

If you want to plot many data sets with different times on a base map, only plot the base map once!

Some psudo code:

fig = plt.figure(1)
ax  = fig.add_subplot(111)

## Draw Background basemap
m = Basemap(resolution='i',projection='cyl',\
    llcrnrlon=bot_left_lon,llcrnrlat=bot_left_lat,\
    urcrnrlon=top_right_lon,urcrnrlat=top_right_lat,)
m.drawcoastlines()
m.drawcountries()
m.drawcounties()
m.drawstates()

# Brackground image from GIS
m.arcgisimage(service='NatGeo_World_Map', xpixels = 1000, dpi=100)

# Add a shape file for Fire perimiter
p4 = m.readshapefile('/uufs/chpc.utah.edu/common/home/u0553130/fire_shape/perim_'+perimiter,'perim',drawbounds=False)

# Plot station locations
m.scatter(lons,lats)
plt.text(lons[0],lats[0],stnid[0])
plt.text(lons[1],lats[1],stnid[1])

for analysis_h in range(1,5):

    # code to get lat/lon and data

    # plot pcolormesh, and colorbar
    PC = m.pcolormesh(lon,lat,data)
    cbar = plt.colorbar(orientation='horizontal',pad=0.05,shrink=.8,ticks=range(0,80,5))

    # plot contour
    CC = m.contour(lon,lat,data)

    # plot HRRR 10-m winds
    WB = m.barbs(lon,lat,u,v,)
    

    ### Before you plot the next loop, remove the data plots
    cbar.remove()                  # remove colorbar
    PC.remove()                    # remove pcolomesh
    WB[0].remove()              #remove wind barbs
    WB[1].remove()              #   (takes two steps)
    for cc in CC.collections:  #remove each line in the contour collection
        cc.remove()

Plotting the fire perimeter from a shape file

There is very useful documentation related to plotting shape files on a basemap here: http://basemaptutorial.readthedocs.io/en/latest/shapefile.html

I needed to plot fire perimeter, plotted on this map of Idaho and shaded in red, on a basemap.

This is how I did it...

from matplotlib.patches import Polygon
from matplotlib.collections import PatchCollection
## My file 
# Plot Fire Perimeter                               
perimiter = '160806'
p4 = m.readshapefile('[PATH]/perim_'+perimiter,'perim',drawbounds=False)
                           
patches   = []
for info, shape in zip(m.perim_info, m.perim):
      if info['FIRENAME'] == 'PIONEER' and info['SHAPENUM']==1772:
             x, y = zip(*shape) 
              #print info
              #m.plot(x, y, marker=None,color='maroon',linewidth=2)
              patches.append(Polygon(np.array(shape), True) )
              ax.add_collection(PatchCollection(patches, facecolor= 'maroon', alpha=.65, edgecolor='k', linewidths=1.5, zorder=1))

Plotting radar data with MetPy, pyproj, Basemap

MetPy radar plots

The MetPy python package is really helpful for atmospheric scientists which allows you to plot radar data (specifically Level 3 data) downloaded from the THREDDS server whose files are in the .nids format. How to use MetPy to plot these data is explained on the MetPy docs, and you'll create a plots of reflectivity and radial velocity like this...

There are two things the example doesn't tell you...

1) The data value units are arbitrary. Notice that the values of dbz color scale is the same as the values in the radial velocity scale on the right. What the MetPy example doesn't tell you is that you need to use the map_data instance to convert the data to appropriate units.

f = Level3File(FILE)

# Pull the data out of the file object (this stuff is in 8-bit)
datadict = f.sym_block[0][0]

# Turn into an array, then mask
data = f.map_data(datadict['data']) # I'm confident this maps the data to floats with units of m/s
data = ma.array(data)
data[np.isnan(data)] = ma.masked

2) You don't get a sense of where the radar and its observations are located in relation to places, just a distance from the center. You can convert the azimuth and range into a latitude and longitude by using the pyproj module. Then you can use the lon and lat variables to plot the data on a basemap. Given the center latitude, center longitude, azimuth (degrees from north), and the range from center lat/lon, you can calculate the latitude and longitude of a distant point like this...

from pyproj import Geod

# Grab azimuths and calculate a range based on number of gates
az = np.array(datadict['start_az'] + [datadict['end_az'][-1]])
rng = np.linspace(0, f.max_range, data.shape[-1] + 1)

# Convert az, range to a lat/lon
g = Geod(ellps='clrk66') # This is the type of ellipse the earth is projected on. 
                         # There are other types of ellipses you can use,
                         # but the differences will be small
center_lat = np.ones([len(az),len(rng)])*f.lat    
center_lon = np.ones([len(az),len(rng)])*f.lon
az2D = np.ones_like(center_lat)*az[:,None]
rng2D = np.ones_like(center_lat)*np.transpose(rng[:,None])*1000
lon,lat,back=g.fwd(center_lon,center_lat,az2D,rng2D)
#...(add code to create map)
m.pcolormesh(lon,lat,data,cmap="BrBG_r",vmax=10,vmin=-10)

You can replace the xlocs and ylocs in the MetPy example with the lon and lat variables, and use those lons/lats to plot on a basemap.
(Note: back_az is the azimuth from the point back to the center location)

TDWR plot on Basemap

I like to look at the TDWR base radial velocity scans to look for lake breezes in the Salt Lake Valley. Below is plotted the radial velocity on June 19, 2015 at 00:24 UTC using MetPy without converting to correct units using map_data and has not been plotted on a basemap (two color schemes used).

both plots are the same, but the color bar is different.

It is so much nicer to view radar data on a basemap plot, as shown below, so you get a sense of there the observations are in space. You can see a lake breeze boundary in the center of the valley (below center of the image where color changes from dark green to light green). I have also used the map_data instance to covert the data to units of ms-1. Each plot is the same with different color schemes. Left is the Brown/BlueGreen that I like, and right is the NWS Red/Green scheme built into MetPy.

 

The full code to do this is on github:
Plotting TDWR Radar data:
https://github.com/blaylockbk/pyBKB_v2/blob/master/TDWR/plot_TDWR_reflectivity.py

Plotting NEXRAD Level 2 data:

https://github.com/blaylockbk/pyPlots_v2/blob/master/NEXRAD_II/plot_NEXRAD_II_data.py

Extra things you can put on the map are a shape file of the roads, and weather observations from the MesoWest API. Let me know if you find this useful :)

Python Basemap Background Image from ArcGIS MapServer

Click here for Jupyter Notebook

https://github.com/blaylockbk/pyBKB_v2/blob/master/demos/add_arcgis_image_to_basemap.ipynb
Python's Basemap toolkit has some default backgrounds you can load. To draw the background on your map you use m.bluemarble(), and m.shadedrelief(). These are really low resolution, so they aren't very useful if you zoom in on a specific area.

Alternatively, you can download higher resolution maps from the ArcGIS maps server. Here's how...
m.arcgisimage(service=maps, xpixels = 3500, dpi=500, verbose= True)

For the Salt Lake Valley, here's what each of the map backgrounds look like. You can plot anything on top of this :) 

 

To plot all of the availble maps I've made this script:
import matplotlib.pyplot as plt
import numpy as np
from functions_domains_models import *

from mpl_toolkits.basemap import Basemap


# Define the map boundaries lat/lon
domain = get_domain('salt_lake_valley')
top_right_lat = domain['top_right_lat']+.1
top_right_lon = domain['top_right_lon']-.1
bot_left_lat = domain['bot_left_lat']
bot_left_lon = domain['bot_left_lon']

## Map in cylindrical projection (data points may apear skewed)
m = Basemap(resolution='i',projection='cyl',\
    llcrnrlon=bot_left_lon,llcrnrlat=bot_left_lat,\
    urcrnrlon=top_right_lon,urcrnrlat=top_right_lat,)


map_list = [
'ESRI_Imagery_World_2D',    # 0
'ESRI_StreetMap_World_2D',  # 1
'NatGeo_World_Map',         # 2
'NGS_Topo_US_2D',           # 3
#'Ocean_Basemap',            # 4
'USA_Topo_Maps',            # 5
'World_Imagery',            # 6
'World_Physical_Map',       # 7     Still blurry
'World_Shaded_Relief',      # 8
'World_Street_Map',         # 9
'World_Terrain_Base',       # 10
'World_Topo_Map'            # 11
]

for maps in map_list: 
    plt.figure(figsize=[10,20])    
    ## Instead of using WRF terrain fields you can get a high resolution image from ESRI
    m.arcgisimage(service=maps, xpixels = 3500, dpi=500, verbose= True)
    m.drawstates()
    plt.title(maps)
    
    plt.savefig('00'+maps,bbox_inches="tight")