How to save output raster to my workspace

[slesinskip]

Hi I'm testing S2 cloudless mosaic script. Every thing works fine but I can't save my output to my workspace as geotiff. After trying to save outputs I got an error:

AttributeError Traceback (most recent call last)
Cell In [23], line 1
----> 1 image.odc.write_cog(
2 fname,
3 overwrite=True,
4 blocksize=2048,
5 ovr_blocksize=1024,
6 overview_resampling="average",
7 intermediate_compression={"compress": "zstd", "zstd_level": 1},
8 use_windowed_writes=True,
9 compress="zstd",
10 zstd_level=6,
11 BIGTIFF=True,
12 SPARSE_OK=True,
13 NUM_THREADS=4,
14 )

File /srv/conda/envs/notebook/lib/python3.10/site-packages/xarray/core/common.py:256, in AttrAccessMixin.getattr(self, name)
254 with suppress(KeyError):
255 return source[name]
--> 256 raise AttributeError(
257 f"{type(self).name!r} object has no attribute {name!r}"
258 )

AttributeError: 'DataArray' object has no attribute 'odc'

Below is my code:

import numpy as np
import xarray as xr
import rasterio.features
import stackstac
import pystac_client
import planetary_computer
import xrspatial.multispectral as ms
from dask_gateway import GatewayCluster

cluster = GatewayCluster() # Creates the Dask Scheduler. Might take a minute.
client = cluster.get_client()
cluster.adapt(minimum=4, maximum=24)
print(cluster.dashboard_link)

area_of_interest = {
"type": "Polygon",
"coordinates": [
[
[14.1226, 49.0022],
[15.1455, 54.8399],
]
],
}
bbox = rasterio.features.bounds(area_of_interest)

stac = pystac_client.Client.open(
"https://planetarycomputer.microsoft.com/api/stac/v1",
modifier=planetary_computer.sign_inplace,
)
search = stac.search(
bbox=bbox,
datetime="2022-03-15/2022-04-25",
collections=["sentinel-2-l2a"],
query={"eo:cloud_cover": {"lt": 25}},
)
items = search.item_collection()
print(len(items))

data = (
stackstac.stack(
items,
assets=["B04", "B03", "B02"], # red, green, blue
chunksize=4096,
resolution=100,
epsg=2180,
)
.where(lambda x: x > 0, other=np.nan) # sentinel-2 uses 0 as nodata
.assign_coords(band=lambda x: x.common_name.rename("bands")) # use common names
)
data

data = data.persist()

median = data.median(dim="time").compute()

image = ms.true_color(*median) # expects red, green, blue DataArrays

import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(8, 8))
ax.set_axis_off()
image.plot.imshow(ax=ax);

fname = f"PL.tif"
print(f"Saving to: '{fname}'")

image.odc.write_cog(
fname,
overwrite=True,
blocksize=2048,
ovr_blocksize=1024,
overview_resampling="average",
intermediate_compression={"compress": "zstd", "zstd_level": 1},
use_windowed_writes=True,
compress="zstd",
zstd_level=6,
BIGTIFF=True,
SPARSE_OK=True,
NUM_THREADS=4,
)

[TomAugspurger]

The line

AttributeError: 'DataArray' object has no attribute 'odc'

indicates that you need to import the odc.geo package to register its accessor (https://odc-geo.readthedocs.io/en/latest/). Adding an import odc.geo should fix things.

By the way, we don't recommend storing outputs (at least large outputs) in your home directory on the Hub. See https://planetarycomputer.microsoft.com/docs/overview/environment/#understanding-the-file-system for more.