import rasterio
import xarray as xr
import hvplot.xarray
from dask.distributed import Client, progress
import cartopy.crs as ccrs
print('Xarray version: ', xr.__version__)
print('Rasterio version: ', rasterio.__version__)
print('Hvplot version: ', hvplot.__version__)
Xarray version: 0.10.8 Rasterio version: 1.0.3 Hvplot version: 0.2.1
client = Client()
client
/opt/conda/lib/python3.6/site-packages/distributed/bokeh/core.py:55: UserWarning: Port 8787 is already in use. Perhaps you already have a cluster running? Hosting the diagnostics dashboard on a random port instead. warnings.warn('\n' + msg)
Client
|
Cluster
|
image_url = 'https://storage.googleapis.com/gcp-public-data-landsat/LC08/01/047/027/LC08_L1TP_047027_20130421_20170310_01_T1/LC08_L1TP_047027_20130421_20170310_01_T1_B4.TIF'
with rasterio.open(image_url) as src:
print(src.profile)
{'driver': 'GTiff', 'dtype': 'uint16', 'nodata': None, 'width': 7751, 'height': 7531, 'count': 1, 'crs': CRS({'init': 'epsg:32610'}), 'transform': Affine(30.0, 0.0, 356685.0, 0.0, -30.0, 5367615.0), 'blockxsize': 256, 'blockysize': 256, 'tiled': True, 'compress': 'lzw', 'interleave': 'band'}
# Note that the blocksize of the image is 256 by 256, so we want xarray to use some multiple of that
xchunk = 2048
ychunk = 2048
da = xr.open_rasterio(image_url, chunks={'band': 1, 'x': xchunk, 'y': ychunk})
da
<xarray.DataArray (band: 1, y: 7531, x: 7751)> dask.array<shape=(1, 7531, 7751), dtype=uint16, chunksize=(1, 2048, 2048)> Coordinates: * band (band) int64 1 * y (y) float64 5.368e+06 5.368e+06 5.368e+06 5.368e+06 5.367e+06 ... * x (x) float64 3.567e+05 3.567e+05 3.568e+05 3.568e+05 3.568e+05 ... Attributes: transform: (30.0, 0.0, 356685.0, 0.0, -30.0, 5367615.0, 0.0, 0.0, 1.0) crs: +init=epsg:32610 res: (30.0, 30.0) is_tiled: 1 nodatavals: (nan,)
band1 = da.sel(band=1).persist()
progress(band1)
VBox()
%%time
display(band1.hvplot(rasterize=True, width=600, height=400, cmap='viridis'))
CPU times: user 4.01 s, sys: 628 ms, total: 4.64 s Wall time: 3.66 s
%%time
crs = ccrs.UTM(zone='10n')
display(band1.hvplot(crs=crs, rasterize=True, width=600, height=400, cmap='viridis'))
CPU times: user 3min 37s, sys: 7.82 s, total: 3min 45s Wall time: 3min 34s