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()