Ok so here is the issue, GOES IR images from OpenDAP have some serious hot pixesls.. Lets fix that
#first load modules.. We are going to use Py-ART as a holder for the data 'cause Py-ART is Cool ok..
import pyart
import sys, os, urllib, urllib2
from matplotlib import pyplot as plt
from matplotlib import colors
import numpy as np
from mpl_toolkits.basemap import Basemap, addcyclic, pyproj
import netCDF4
import copy
from scipy import ndimage
%matplotlib inline
Ok.. here are some functions.. One checks the opendap directory from ucar.. and works out the latest file.. These are not ideal.. the data is stored in unsigned integers.. and are counts.. -If you know a better opendap source let me know!-
Next function fills out a Py-ART Grid object with Lat Lon data and IR data in counts.. Again, I REALLY want to find out the scaling from counts to irradience! Oh and for those that now know how to use PyProj to convert GOES coords to Lat Lon: You are Welcome!
Why a (Py-ART object? Because we have done the hard yards in IO and this means I can now just use PyART.IO to write the file out..
Final function use Matplotlib's basemap to plot to a nice custom map..
def fetch_gini(dir_string = 'http://motherlode.ucar.edu:8080/thredds/dodsC/satellite/IR/EAST-CONUS_4km/current/',
pattern_match = 'EAST-CONUS_4km_IR_'):
#dir_string = 'http://motherlode.ucar.edu:8080/thredds/dodsC/satellite/3.9/EAST-CONUS_4km/current/'
#pattern_match = 'EAST-CONUS_4km_3.9_'
dirlisting = urllib2.urlopen(dir_string)
all_lines = dirlisting.read().splitlines() #EAST-CONUS_4km_3.9_20140123_2245.gini
has_gini = []
for line in all_lines:
if '.gini' in line:
#print line
happy_place = line.find(pattern_match)
end_place = line.find('.gini')
my_upper_name = line[happy_place:end_place]+'.gini'
url = dir_string + my_upper_name
has_gini.append(url)
return has_gini[0]
def gini_grid(open_dap_dataset):
xg, yg = np.meshgrid(open_dap_dataset.variables['x'][:]*1000.0, open_dap_dataset.variables['y'][:]*1000.0)
vals = np.uint8(open_dap_dataset.variables['IR'][:])
pnyc = pyproj.Proj(proj = 'lcc',
lat_1 = open_dap_dataset.variables['LambertConformal'].latitude_of_projection_origin,
lat_2 = open_dap_dataset.variables['LambertConformal'].latitude_of_projection_origin,
lat_0 = open_dap_dataset.variables['LambertConformal'].latitude_of_projection_origin,
lon_0 = open_dap_dataset.variables['LambertConformal'].longitude_of_central_meridian )
lon, lat = pnyc(xg, yg, inverse = True)
mng = pyart.testing.make_empty_grid([1, xg.shape[0], xg.shape[1]],
( (0,0), (yg.min(),yg.max()),(xg.min(), xg.max()) ))
mng.axes['lat']['data'] = open_dap_dataset.variables['LambertConformal'].latitude_of_projection_origin
mng.axes['lon']['data'] = open_dap_dataset.variables['LambertConformal'].longitude_of_central_meridian
mng.axes['alt']['data'] = 0.0
ir_fld = { 'data' : vals,
'units' : 'Counts',
'standard_name' : 'Sensor Counts', #non CF
'long_name' : 'Number_of_counts_in_channel',
'valid_max' : 256,
'valid_min' : 0,
'_FillValue' : 9999}
lat_fld = {'data' : lat,
'units' : 'degree_north',
'standard_name' : 'latitude',
'long_name' : 'latitude_of_grid_point',
'valid_max' : 90,
'valid_min' : -90}
lon_fld = {'data' : lon,
'units' : 'degree_east',
'standard_name' : 'longitude',
'long_name' : 'longitude_of_grid_point',
'valid_max' : 180,
'valid_min' : -180}
mng.fields.update({'IR' : ir_fld, 'lat': lat_fld, 'lon': lon_fld})
return mng
def plot_gini_grid(gini_grid, box = [-110, -70, 20, 52], resolution = 'l',
parallels = np.linspace(10,50, 9),
meridians = np.linspace(-110, -80,7),
vmin = None, vmax = None,
fld = 'IR'):
m = Basemap(llcrnrlon = box[0] ,llcrnrlat = box[2] , urcrnrlon = box[1],
urcrnrlat = box[3] , projection = 'mill', area_thresh =1000 ,
resolution='l')
x, y = m(gini_grid.fields['lon']['data'], gini_grid.fields['lat']['data'])
# create figure.
m.drawparallels(np.linspace(10,50, 9) ,labels=[1,1,0,0])
m.drawmeridians(np.linspace(-110, -80,7),labels=[0,0,0,1])
pc = m.pcolormesh(x, y , gini_grid.fields[fld]['data'][0,:], cmap=plt.get_cmap('gray'),
vmin = vmin, vmax = vmax)
m.drawcoastlines(linewidth=1.25)
m.drawstates()
plt.colorbar(mappable=pc)
Ok lets read, dump tp grid and plot!
open_dap_11microns = fetch_gini()
data_11=netCDF4.Dataset(open_dap_11microns)
grid_11 = gini_grid(data_11)
f = plt.figure(figsize = [15,11])
plot_gini_grid(grid_11, box = [-90, -85, 40, 45], vmin = 105, resolution = 'h', fld = 'IR')
Crud! That looks problematic.. Do some googling and you find that the 11 micron imager has some issues.. So how do we heal this. Here is my take.. all those hot pixels are over 195 ish counts.. to lets selectively apply a median filter..
new_ir = copy.deepcopy(grid_11.fields['IR'])
img = new_ir['data'][0,:]
level = 80
img[np.where(img > img.max()-level)] = ndimage.median_filter(img, 2)[np.where(img > img.max()-level)]
new_ir['data'][0, :] = img
grid_11.fields.update({'IR_filt' : new_ir})
Lets PLOT!
f = plt.figure(figsize = [15,11])
plot_gini_grid(grid_11, box = [-90, -85, 40, 45], vmin = 105, resolution = 'h', fld = 'IR_filt')