This notebook needs the dev version of rasterio
conda install -c conda-forge/label/dev rasterio
import xarray as xr
The NLCD dataset was downloaded and transformed to cloud-optimized geotiff thusly:
wget --no-clobber --no-verbose \
'http://www.landfire.gov/bulk/downloadfile.php?TYPE=nlcd2011&FNAME=nlcd_2011_landcover_2011_edition_2014_10_10.zip' \
--output-document=NLCD2011.zip
unzip NLCD2011.zip
cd nlcd_2011_landcover_2011_edition_2014_10_10
conda install gdal
gdal_translate nlcd_2011_landcover_2011_edition_2014_10_10.img \
nlcd_2011_landcover_2011_edition_2014_10_10.tif -co TILED=YES -co COMPRESS=LZW
rclone copy nlcd_2011_landcover_2011_edition_2014_10_10.tif s3:rsignell/NLCD/
For private access buckets, need to provide AWS credentials (either in ~/.config/aws or here)
%env AWS_ACCESS_KEY_ID= AWS_SECRET_ACCESS_KEY=
env: AWS_ACCESS_KEY_ID=AWS_SECRET_ACCESS_KEY=
fname = "s3://rsignell/NLCD/nlcd_2011_landcover_2011_edition_2014_10_10.tif"
For public access buckets, we need to currently use the HTTPS URL as per this github issue: https://github.com/mapbox/rasterio/issues/1362#issuecomment-393653563
fname = "https://rsignell.s3.amazonaws.com/NLCD/nlcd_2011_landcover_2011_edition_2014_10_10.tif"
ds = xr.open_rasterio(fname, chunks={'x':256, 'y':256})
ds
<xarray.DataArray (band: 1, y: 104424, x: 161190)> dask.array<shape=(1, 104424, 161190), dtype=uint8, chunksize=(1, 256, 256)> Coordinates: * band (band) int64 1 * y (y) float64 3.31e+06 3.31e+06 3.31e+06 3.31e+06 3.31e+06 ... * x (x) float64 -2.493e+06 -2.493e+06 -2.493e+06 -2.493e+06 ... Attributes: transform: (-2493045.0, 30.0, 0.0, 3310005.0, 0.0, -30.0) crs: +ellps=GRS80 +lat_0=23 +lat_1=29.5 +lat_2=45.5 +lon_0=-96 +n... res: (30.0, 30.0) is_tiled: 1 nodatavals: (nan,)
ds.crs
'+ellps=GRS80 +lat_0=23 +lat_1=29.5 +lat_2=45.5 +lon_0=-96 +no_defs +proj=aea +towgs84=0,0,0,0,0,0,0 +units=m +x_0=0 +y_0=0'
ds[0,50000:51000,80000:81000].plot.imshow();
# top left corner of region of interest
e0 = -93030.; n0 = 1809990
# size of region in meters
xsize, ysize = 30000., 30000.
e1 = e0 + xsize
n1 = n0 - ysize
ds[0,:,:].sel(x=slice(e0,e1), y=slice(n0,n1)).plot.imshow();