Download Satellite Images from NASA Worldview, Add a watermark with the date

Partly taken form some old "watermark" code I found several years ago, this script allows you do download MODIS images from NASA's Worldview image viewer. I then compiled the images into an animated GIF.

# Brian Blaylock
# University of Utah

# Download True Color Images from NASA WorldView
# and add the time stamp to the image

## Images downloaded from NASA's Worldview
## https://earthdata.nasa.gov/labs/worldview/
# 1) zoom in to desired location
# 2) click "take a snapshot" icon
# 3) Draw the region you wish to take a photo
# 4) Modify the URL below, particularly the lat and lon
# 5) Edit the start and end dates

import urllib
from datetime import datetime, timedelta
import numpy as np

from PIL import Image, ImageDraw, ImageFont, ImageEnhance

FONT = 'Arial.ttf'
watermark_label = "THE DATE GOES HERE"


def add_watermark(in_file, text, out_file='watermark.jpg', angle=0, opacity=0.8):
  img = Image.open(in_file).convert('RGB')
  watermark = Image.new('RGBA', img.size, (0,0,0,0))
  size = 2
  n_font = ImageFont.truetype(FONT, size)
  n_width, n_height = n_font.getsize(text)
  text_size_scale = 1.5
  while n_width+n_height < watermark.size[0]/text_size_scale:
    size += 2
    n_font = ImageFont.truetype(FONT, size)
    n_width, n_height = n_font.getsize(text)
  draw = ImageDraw.Draw(watermark, 'RGBA')
  draw.text(((watermark.size[0] - n_width) / 10,
            (watermark.size[1] - n_height) / 100),
  text, font=n_font)
  watermark = watermark.rotate(angle,Image.BICUBIC)
  alpha = watermark.split()[3 ]
  alpha = ImageEnhance.Brightness(alpha).enhance(opacity)
  watermark.putalpha(alpha)
  Image.composite(watermark, img, watermark,).save(out_file, 'JPEG')

##--------------------------------------

outdir = 'satellite_img/'

start_date = datetime(2015,8,15)
end_date = datetime(2015,8,25)


# specify the dates you want to retrieve

date = start_date

while end_date >= date:
  for sat in np.array(["Terra","Aqua"]):
    year = str(date.year)
    dayofyear = str(date.timetuple().tm_yday)
    stringdate = datetime.strftime(date,"%Y-%m-%d")

    URL = "http://map2.vis.earthdata.nasa.gov/image-downloadTIME="+year+dayofyear+"&extent=-116.6972484336448,35.054396923451016,-106.9589671836448,43.896193798451016&epsg=4326&layers=MODIS_"+sat+"_CorrectedReflectance_TrueColor,Coastlines,Reference_Features,MODIS_Fires_All&opacities=1,1,1,1&worldfile=false&format=image/jpeg&width=4433&height=4024"

    #Since the Terra satellite passes over before Aqua
    # we need to save it before Aqua (alphebetical order puts it in wrong order)
    if sat == "Terra":
      sat_order="1Terra"
    if sat == "Aqua":
      sat_order="2Aqua"
    image_name = 'MODIS_TrueColor_'+stringdate+'_'+sat_order+'.jpg'
    urllib.urlretrieve(URL, outdir+image_name)
    add_watermark(outdir+image_name,'MODIS '+sat+' True Color: '+stringdate,outdir+image_name)
    print "Saved", image_name
  date = date+timedelta(days=1)

Quick WRF Plot

Plotting WRF Data in python without NetCDF4:

Scrapped this code from places on line. Reads in netcdf data with scipy. Plots specific humidity with the lake mask outlined on top. See blog post here for details on the plot.

# Brian Blaylock
#
# Plotting WRF netCDF output
# Parts and Pieces of this code is from Luke Madaus (University of Washington)


import sys,getopt
#from netCDF4 import Dataset # use scipy instead
from scipy.io import netcdf
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import os
from datetime import datetime, timedelta
import coltbls as coltbls
from mpl_toolkits.basemap import Basemap
from matplotlib.colors import LinearSegmentedColormap
from mpl_toolkits.axes_grid import make_axes_locatable
import matplotlib.axes as maxes


# Running Script example:
# $ python my_WRF_map.py wrfout_d0x_YYYY-MM-DD_HH:MM:SS

# Open file in a netCDF reader
directory = ''
wrf_file_name = 'M11KS3~E.NC' #NC files are named differently in a windows system
nc = netcdf.netcdf_file(wrf_file_name,'r')

# Grab these variables for now
landmask = nc.variables['LANDMASK'] # HRRR metgrid landmask




#-------------------------------------------
# BEGIN ACTUAL PROCESSING HERE
#-------------------------------------------

# x_dim and y_dim are the x and y dimensions of the model
# domain in gridpoints
x_dim = nc.dimensions['west_east']
y_dim = nc.dimensions['south_north']

# Get the grid spacing
dx = float(nc.DX)
dy = float(nc.DY)

width_meters = dx * (x_dim - 1) #Domain Width
height_meters = dy * (y_dim - 1) #Domain Height

cen_lat = float(nc.CEN_LAT)
cen_lon = float(nc.CEN_LON)
truelat1 = float(nc.TRUELAT1)
truelat2 = float(nc.TRUELAT2)
standlon = float(nc.STAND_LON)

# Draw the base map behind it with the lats and
# lons calculated earlier
m = Basemap(resolution='i',projection='lcc',\
width=width_meters,height=height_meters,\
lat_0=cen_lat,lon_0=cen_lon,lat_1=truelat1,\
lat_2=truelat2)

# This sets the standard grid point structure at full resolution
# Converts WRF lat and long to the maps x an y coordinate
x,y = m(nc.variables['XLONG_M'][0],nc.variables['XLAT_M'][0])


# This sets a thinn-ed out grid point structure for plotting
# wind barbs at the interval specified in "thin"
thin = 5
x_th,y_th = m(nc.variables['XLONG_M'][0,::thin,::thin],\
nc.variables['XLAT_M'][0,::thin,::thin])

# Set universal figure margins
width = 10
height = 8

plt.figure(figsize=(width,height))
plt.rc("figure.subplot", left = .001) #gets rid of white space in plot
plt.rc("figure.subplot", right = .999)
plt.rc("figure.subplot", bottom = .001)
plt.rc("figure.subplot", top = .999)

F = plt.gcf() # Gets the current figure

m.drawstates(color='k', linewidth=1.25)
m.drawcoastlines(color='k')
m.drawcountries(color='k', linewidth=1.25)

height = nc.variables['HGT_M'][0,:,:]
landmask = nc.variables['LANDMASK'][0,:,:]
spechumd = nc.variables['SPECHUMD'][0][0,:,:]

plt.figure(1)
plt.title("HRRR Land Mask and Specific Humidity\nAugust 10, 2015 23:00z")
plt.pcolormesh(x,y,landmask)
plt.pcolormesh(x,y,spechumd*1000, cmap = plt.cm.cool)
cbar_loc = plt.colorbar()
cbar_loc.ax.set_ylabel('Specific Humidity g/kg')
plt.contour(x,y,landmask, [0,1], linewidths=1, colors="k")
#plt.contour(x,y,height)

xs,ys = m(-111.97,40.78) #buffr sounding coordinates
plt.scatter(xs,ys, s=70, c='w')







plt.show()

Plotting Sounding Data from University of Wyoming's Website

Plotting sounding data can be difficult, but it is made easy with this python module: https://pypi.python.org/pypi/SkewT

It requires a .txt file with data in the format of the University of Wyoming sounding data archive (example) or you can create your own dictionary.

Here I grab the sounding data from the University of Wyoming's website and process it with the following steps:

1. Use urllib2 to open the url and read the data
2. Parse out the html tags using BeautifulSoup
3. Separate the text by a new line "\n"
4. Write each line of text to a new .txt file reinserting the new line with +"\n" and skipping the first three lines of the file (this puts it in the same format as the University of Wyoming website)
5. Draw the Skew-T Plot (NOTE: requires the SkewT package here)

Here is the download function on github. 
As always, please share other solutions you may come up with :)

# Brian Blaylock

# July 6, 2015

# Download, process, and Plot Sounding Data from University of Wyoming

import urllib2
from bs4 import BeautifulSoup
from skewt import SkewT

stn = '72572' #72572 is ID for SLC
year= '2015'
month = '06'
day = '12'
hour = '12' #either 12 or 00

# 1)
# Wyoming URL to download Sounding from
url = 'http://weather.uwyo.edu/cgi-bin/sounding?region=naconf&TYPE=TEXT%3ALIST&YEAR='+year+'&MONTH='+month+'&FROM='+day+hour+'&TO='+day+hour+'&STNM='+stn
content = urllib2.urlopen(url).read()

# 2)
# Remove the html tags
soup = BeautifulSoup(content)
data_text = soup.get_text()

# 3)
# Split the content by new line.
splitted = data_text.split("\n",data_text.count("\n"))

# 4)
# Write this splitted text to a .txt document
Sounding_filename = str(stn)+'.'+str(year)+str(month)+str(day)+str(hour)+'.txt'
f = open(Sounding_filename,'w')
for line in splitted[4:]:
    f.write(line+'\n')
f.close()

# 5)
sounding = SkewT.Sounding(filename=Sounding_filename)
sounding.plot_skewt()


And if you wish to plot more than one sounding on the same plot do this...
S = SkewT.Sounding(filename=Sounding_filename)
T = SkewT.Sounding(filename="72572.2015061212.txt")
S.make_skewt_axes(tmax=55)
S.add_profile(color='r',linewidth=5,bloc=0)
S.soundingdata=T.soundingdata
S. add_profile(color="b",bloc=1)

Just made my first python website with cgi

http://meso2.chpc.utah.edu/gslso3s/cgi-bin/ozone_rose.cgi

Check out the neat webpage I made using cgi and python. This page allows a user to generate ozone pollution roses based on a range of dates. It uses the windrose.py functions to make the plots. Give it a try!!

Python Legend, put legend outside of plot display

Often when you legend is too large, it is convenient to put it off to the side. Here's how...

from matplotlib import pyplot as plt

plt.plot(a,b, label="legend label")
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5),prop={'size':10})

An example here: