<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Transforming points in one crs to another, e.g. OSGB (Ordinance Survey Great Britain to Geodetic
Step1:  You can project lines with multiple vertices also
Step2:  Look at this line on a map
Step3:  Using cartopy's io shapefile interface to read Hunza basin shapefile
<ASSISTANT_TASK:>
Python Code:
osgb = ccrs.OSGB()
geod = ccrs.Geodetic()
# Convert from Ordnance Survey GB to lon/lat:
easting = 291813.424
northing = 92098.387
lon, lat = geod.transform_point(
    x=easting, y=northing, src_crs=osgb)
print(lon, lat)
# check with mapx, this is UL corner of EASE-Grid 2.0 N!
e2n = ccrs.LambertAzimuthalEqualArea(central_latitude=90.0)
lon, lat = geod.transform_point(
    x = -9000000.,
    y = 9000000.,
    src_crs=e2n)
print(lon, lat)
Explanation: Transforming points in one crs to another, e.g. OSGB (Ordinance Survey Great Britain to Geodetic
End of explanation
import shapely.geometry as sgeom
new_york = [-74.0060, 40.7128] # lon, lat
honolulu = [-157.8583, 21.3069] # lon, lat
line = sgeom.LineString([new_york, honolulu])
pc = ccrs.PlateCarree()
lines = pc.project_geometry(line, pc)
lines
help(lines)
Explanation: You can project lines with multiple vertices also:
End of explanation
%matplotlib notebook
import matplotlib.pyplot as plt
plt.figure()
ax = plt.axes(projection=pc)
ax.add_geometries(
    [lines], pc,
    edgecolor='blue', facecolor='none', lw=2)
ax.coastlines()
plt.show()
# Will this work on E2N? No:  I think it's because original line coords are in lat/lon?
#e2n_lines = e2n.project_geometry(line, e2n)
#plt.figure()
#ax = plt.axes(projection=e2n)
#ax.add_geometries(
#    [e2n_lines], e2n,
#    edgecolor='blue', facecolor='none', lw=2)
#ax.coastlines()
#plt.show()
Explanation: Look at this line on a map:
End of explanation
import cartopy.crs as ccrs
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)
e2n = ccrs.LambertAzimuthalEqualArea(central_latitude=90.)
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)
ax.set_extent([-9000000., 9000000., -9000000., 9000000.], crs=e2n)
ax.gridlines(color='gray', linestyle='--')
ax.coastlines()
# ax.imshow(tb, transform=e2n)
#ax.add_geometries(
#    [e2n_basin], e2n,
#    edgecolors='red', facecolor='none', lw=2)
plt.tight_layout()
plt.show()
#ax.imshow(tb, extent=extent, transform=rob, origin='upper')
Explanation: Using cartopy's io shapefile interface to read Hunza basin shapefile
End of explanation
<END_TASK>