import urllib2

import cartopy.crs as ccrs
from cartopy.feature import NaturalEarthFeature
import matplotlib.pyplot as plt
from netCDF4 import Dataset

def fetch_gini(dir_string=('http://motherlode.ucar.edu:8080/'
                           'thredds/dodsC/satellite/IR/EAST-CONUS_4km/current/'),
               pattern_match = 'EAST-CONUS_4km_IR_'):
    dirlisting = urllib2.urlopen(dir_string)
    all_lines = dirlisting.read().splitlines()
    has_gini = []
    for line in all_lines:
        if '.gini' in 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]

url = fetch_gini()
print 'retreiving data from: {}'.format(url)
ds = Dataset(url)
ir = ds.variables['IR'][0].astype(np.uint8)
x = ds.variables['x'][:] * 1000
y = ds.variables['y'][:] * 1000

grid = ds.variables['LambertConformal']
lon0 = grid.longitude_of_central_meridian
lat0 = grid.latitude_of_projection_origin

data_crs = ccrs.LambertConformal(central_latitude=lat0, central_longitude=lon0,
                                 secant_latitudes=(lat0, lat0))

fig = plt.figure(figsize=(15, 11))

# Miller projection map over Lake Michigan:
ax = plt.axes(projection=ccrs.Miller())
ax.set_extent([-90, -85, 40, 45], crs=ccrs.Geodetic())

# Add state boundaries etc.:
us_states = NaturalEarthFeature('cultural',
                                'admin_1_states_provinces_lakes_shp',
                                '50m',
                                facecolor='none')
ax.add_feature(us_states)

# Draw some grid lines:
ax.gridlines()

# Plot the IR data, making sure to specify the data CRS:
pc = ax.pcolormesh(x, y, ir, transform=data_crs,
                   cmap=plt.get_cmap('gray'), vmin=105)

# Add a colorbar:
cb = plt.colorbar(pc, orientation='vertical')

fig = plt.figure(figsize=(15, 11))

# Orthographic projection:
ax = plt.axes(projection=ccrs.Orthographic(central_longitude=-95,
                                           central_latitude=25))
ax.set_global()

# Draw coastlines on the map:
ax.coastlines()

# Draw some grid lines:
ax.gridlines()

# Plot the IR data, making sure to specify the data CRS:
pc = ax.pcolormesh(x, y, ir, transform=data_crs,
                   cmap=plt.get_cmap('gray'))

# Add a colorbar:
cb = plt.colorbar(pc, orientation='vertical')

# Plot corner points accoriding to http://www.nws.noaa.gov/noaaport/html/icdtb48e.html:
for point in [(-113.133, 16.369), (-65.091, 14.335),
              (-49.385, 57.289), (-123.044, 59.844)]:
    ax.plot(point[0], point[1], transform=ccrs.Geodetic(),
            markersize=10, marker='o', markerfacecolor='y')