import cartopy.io.shapereader as shpreader
import shapely.geometry as sgeom

bfile = '/Users/brodzik/Desktop/GIS_data/basins/IN_Hunza_at_DainyorBridge.shp'

reader = shpreader.Reader(bfile)
record = next(reader.records())

record

record.attributes

record.bounds

record.geometry

help(record)

"""
Explanation: Using cartopy's io shapefile interface to read Hunza basin shapefile
End of explanation
"""


e2n_basin = e2n.project_geometry(record.geometry)

plt.figure()

ax = plt.axes(projection=e2n)
ax.set_extent([0., 9000000., -9000000., 0.], crs=e2n)
ax.gridlines(color='gray', linestyle='--')
ax.coastlines()

ax.add_geometries(
    [e2n_basin], e2n,
    edgecolors='red', facecolor='none', lw=2)


plt.tight_layout()
plt.show()


from netCDF4 import Dataset

cetbfile = '/Users/brodzik/cetb_data/v1.3/AQUA_AMSRE/2004164/NSIDC-0630-EASE2_N3.125km-AQUA_AMSRE-2004164-36H-M-SIR-RSS-v1.3.nc'
f = Dataset(cetbfile, 'r', 'netCDF4')
f

import numpy as np
tb = np.squeeze(f.variables['TB'][:])
tb.shape

plt.figure()

ax = plt.axes(projection=e2n)
extent = [-9000000., 9000000., -9000000., 9000000.] # [xmin, xmax, ymin, ymax] of UL corner
ax.set_extent(extent, crs=e2n)
ax.gridlines(color='gray', linestyle='--')
ax.coastlines()

ax.imshow(tb, extent=extent, transform=e2n, origin='upper')

ax.add_geometries(
    [e2n_basin], e2n,
    edgecolors='red', facecolor='none', lw=2)

plt.tight_layout()
plt.show()


"""
Explanation: What does it mean to get a projected geometry?
Converts the shapefile geod  to E2N?
End of explanation
"""


