diff --git a/docs/api/main.rst b/docs/api/main.rst index 12ade73..99eb9a5 100644 --- a/docs/api/main.rst +++ b/docs/api/main.rst @@ -46,4 +46,4 @@ Utilities to work with geospatial DataArrays produced by stackstac. reproject_array xyztile_of_array array_bounds - array_epsg + array_crs diff --git a/stackstac/__init__.py b/stackstac/__init__.py index b179390..f5a27cc 100644 --- a/stackstac/__init__.py +++ b/stackstac/__init__.py @@ -2,7 +2,7 @@ from .rio_reader import DEFAULT_GDAL_ENV, MULTITHREADED_DRIVER_ALLOWLIST from .stack import stack from .ops import mosaic -from .geom_utils import reproject_array, array_bounds, array_epsg, xyztile_of_array +from .geom_utils import reproject_array, array_bounds, array_crs, xyztile_of_array try: from . import show as _show # helpful for debugging @@ -42,7 +42,7 @@ def _missing_imports(*args, **kwargs): "mosaic", "reproject_array", "array_bounds", - "array_epsg", + "array_crs", "xyztile_of_array", "__version__", ] diff --git a/stackstac/geom_utils.py b/stackstac/geom_utils.py index b625a4b..40afc3d 100644 --- a/stackstac/geom_utils.py +++ b/stackstac/geom_utils.py @@ -16,7 +16,7 @@ def bounds_from_affine( - af: affine.Affine, ysize: int, xsize: int, from_epsg: int, to_epsg: int + af: affine.Affine, ysize: int, xsize: int, from_crs: Union[str, int], to_crs: Union[str, int] ) -> Bbox: ul_x, ul_y = af * (0, 0) ll_x, ll_y = af * (0, ysize) @@ -26,8 +26,8 @@ def bounds_from_affine( xs = [ul_x, ll_x, lr_x, ur_x] ys = [ul_y, ll_y, lr_y, ur_y] - if from_epsg != to_epsg: - transformer = cached_transformer(from_epsg, to_epsg, always_xy=True) + if from_crs != to_crs: + transformer = cached_transformer(from_crs, to_crs, always_xy=True) # TODO handle error xs_proj, ys_proj = transformer.transform(xs, ys, errcheck=True) else: @@ -39,8 +39,8 @@ def bounds_from_affine( # TODO: use `rio.warp.transform_bounds` instead, for densification? or add our own densification? # ours is likely faster due to transformer caching. -def reproject_bounds(bounds: Bbox, from_epsg: int, to_epsg: int) -> Bbox: - if from_epsg == to_epsg: +def reproject_bounds(bounds: Bbox, from_crs: Union[str, int], to_crs: Union[str, int]) -> Bbox: + if from_crs == to_crs: return bounds minx, miny, maxx, maxy = bounds @@ -48,7 +48,7 @@ def reproject_bounds(bounds: Bbox, from_epsg: int, to_epsg: int) -> Bbox: # read this in pairs, downward by column xs = [minx, minx, maxx, maxx] ys = [maxy, miny, miny, maxy] - transformer = cached_transformer(from_epsg, to_epsg, always_xy=True) + transformer = cached_transformer(from_crs, to_crs, always_xy=True) xs_proj, ys_proj = transformer.transform(xs, ys, errcheck=True) # TODO handle error return min(xs_proj), min(ys_proj), max(xs_proj), max(ys_proj) @@ -80,11 +80,11 @@ def snapped_bounds(bounds: Bbox, resolutions_xy: Resolutions) -> Bbox: return (minx, miny, maxx, maxy) -def array_epsg( +def array_crs( arr: xr.DataArray, default: Union[int, NO_DEFAULT_LITERAL] = NO_DEFAULT ) -> int: """ - Get the coordinate reference system of a `~xarray.DataArray` as an EPSG code. + Get the coordinate reference system of a `~xarray.DataArray` as an CRS code. Parameters ---------- @@ -97,33 +97,33 @@ def array_epsg( Returns ------- int: - EPSG code + CRS code Raises ------ ValueError: - If the `~xarray.DataArray` has no ``epsg`` coordinate, and ``default`` is not given. + If the `~xarray.DataArray` has no ``crs`` coordinate, and ``default`` is not given. """ # TODO look at `crs` in attrs; more compatibility with rioxarray data model try: - epsg = arr.epsg + crs = arr.crs except AttributeError: if default != NO_DEFAULT: return default else: - return epsg.item() + return crs.item() # NOTE: raise out here for shorter traceback raise ValueError( - f"DataArray {arr.name!r} does not have the `epsg` coordinate set, " + f"DataArray {arr.name!r} does not have the `crs` coordinate set, " "so we don't know what coordinate reference system it's in.\n" - "Please set it using `arr.assign_coords(epsg=)` " + "Please set it using `arr.assign_coords(crs=)` " "(note that this returns a new object)." ) -def array_bounds(arr: xr.DataArray, to_epsg: Optional[int] = None) -> Bbox: +def array_bounds(arr: xr.DataArray, to_crs: Optional[int] = None) -> Bbox: """ Get the bounds of a `~xarray.DataArray`. @@ -138,10 +138,10 @@ def array_bounds(arr: xr.DataArray, to_epsg: Optional[int] = None) -> Bbox: The ``x`` and ``y`` coordinates are assumed to indicate the top-left corner of each pixel, not the center. They're also assumed to be evenly spaced (constant resolution), and must be monotonic. - to_epsg: - CRS to reproject the bounds to, as an EPSG code. If None (default), + to_crs: + CRS to reproject the bounds to, as an CRS code. If None (default), the bounds are not reprojected. If not None, the `~xarray.DataArray` must have - an ``epsg`` coordinate. + an ``crs`` coordinate. Returns ------- @@ -157,7 +157,7 @@ def array_bounds(arr: xr.DataArray, to_epsg: Optional[int] = None) -> Bbox: ValueError: If ``x`` or ``y`` coordinates are not monotonic. ValueError: - If ``to_epsg`` is given, and the CRS of the ``~xarray.DataArray`` cannot be determined by `array_epsg`. + If ``to_crs`` is given, and the CRS of the ``~xarray.DataArray`` cannot be determined by `array_crs`. """ @@ -197,10 +197,10 @@ def array_bounds(arr: xr.DataArray, to_epsg: Optional[int] = None) -> Bbox: max(ys), ) - if to_epsg is None: + if to_crs is None: return bounds - return reproject_bounds(bounds, array_epsg(arr), to_epsg) + return reproject_bounds(bounds, array_crs(arr), to_crs) def reproject_array( @@ -227,7 +227,7 @@ def reproject_array( Parameters ---------- arr: - Array to reproject. It must have ``epsg``, ``x``, and ``y`` coordinates. + Array to reproject. It must have ``crs``, ``x``, and ``y`` coordinates. The ``x`` and ``y`` coordinates are assumed to indicate the top-left corner of each pixel, not the center. spec: @@ -247,9 +247,9 @@ def reproject_array( # so this both won't scale, and can be crazy slow in dask graph construction # (and the rechunk probably eliminates any hope of sending an HLG to the scheduler). - from_epsg = array_epsg(arr) + from_crs = array_crs(arr) if ( - from_epsg == spec.epsg + from_crs == spec.crs and array_bounds(arr) == spec.bounds and arr.shape[:-2] == spec.shape ): @@ -273,7 +273,7 @@ def reproject_array( x = np.linspace(minx, maxx, width, endpoint=False) y = np.linspace(maxy, miny, height, endpoint=False) - if from_epsg == spec.epsg: + if from_crs == spec.crs: # Simpler case: just interpolate within the same CRS result = arr.interp( x=x, y=y, method=interpolation, kwargs=dict(fill_value=fill_value) @@ -284,7 +284,7 @@ def reproject_array( # We do this by, for each point in the output grid, generating # the coordinates in the _input_ CRS that correspond to that point. - reverse_transformer = cached_transformer(spec.epsg, from_epsg, always_xy=True) + reverse_transformer = cached_transformer(spec.crs, from_crs, always_xy=True) xs, ys = np.meshgrid(x, y, copy=False) src_xs, src_ys = reverse_transformer.transform(xs, ys, errcheck=True) @@ -293,8 +293,8 @@ def reproject_array( ys_indexer = xr.DataArray(src_ys, dims=["y", "x"], coords=dict(y=y, x=x)) # TODO maybe just drop old dims instead? - old_xdim = f"x_{from_epsg}" - old_ydim = f"y_{from_epsg}" + old_xdim = f"x_{from_crs}" + old_ydim = f"y_{from_crs}" result = arr.rename(x=old_xdim, y=old_ydim).interp( {old_xdim: xs_indexer, old_ydim: ys_indexer}, @@ -329,7 +329,7 @@ def xyztile_of_array( Parameters ---------- arr: - DataArray to extract a map tile from. Must have ``x`` and ``y``, and the ``epsg`` coordinate set. + DataArray to extract a map tile from. Must have ``x`` and ``y``, and the ``crs`` coordinate set. x: X index of the tile y: @@ -356,7 +356,7 @@ def xyztile_of_array( ) bounds = mercantile.xy_bounds(mercantile.Tile(x, y, z)) - if not bounds_overlap(bounds, array_bounds(arr, to_epsg=3857)): + if not bounds_overlap(bounds, array_bounds(arr, to_crs=3857)): return None minx, miny, maxx, maxy = bounds @@ -367,7 +367,7 @@ def xyztile_of_array( tile = reproject_array( arr, RasterSpec( - epsg=3857, + crs=3857, bounds=bounds, resolutions_xy=( (maxx - minx) / tilesize, diff --git a/stackstac/prepare.py b/stackstac/prepare.py index 53a724a..48d91ce 100644 --- a/stackstac/prepare.py +++ b/stackstac/prepare.py @@ -61,7 +61,7 @@ def prepare_items( assets: Optional[Union[List[str], AbstractSet[str]]] = frozenset( ["image/tiff", "image/x.geotiff", "image/vnd.stac.geotiff", "image/jp2"] ), - epsg: Optional[int] = None, + crs: Optional[Union[str, int]] = None, resolution: Optional[Union[IntFloat, Resolutions]] = None, bounds: Optional[Bbox] = None, bounds_latlon: Optional[Bbox] = None, @@ -74,7 +74,7 @@ def prepare_items( f"Cannot give both `bounds` {bounds} and `bounds_latlon` {bounds_latlon}." ) - out_epsg = epsg + out_crs = crs out_bounds = bounds if resolution is not None and not isinstance(resolution, tuple): resolution = (resolution, resolution) @@ -131,7 +131,7 @@ def prepare_items( raise ValueError("Zero asset IDs requested") for item_i, item in enumerate(items): - item_epsg = item["properties"].get("proj:epsg") + item_crs = item["properties"].get("proj:code", item["properties"].get("proj:epsg")) item_bbox = item["properties"].get("proj:bbox") item_shape = item["properties"].get("proj:shape") item_transform = item["properties"].get("proj:transform") @@ -143,7 +143,7 @@ def prepare_items( except KeyError: continue - asset_epsg = asset.get("proj:epsg", item_epsg) + asset_crs = asset.get("proj:code", asset.get("proj:epsg", item_crs)) asset_bbox = asset.get("proj:bbox", item_bbox) asset_shape = asset.get("proj:shape", item_shape) asset_transform = asset.get("proj:transform", item_transform) @@ -185,29 +185,29 @@ def prepare_items( asset_affine = None # Auto-compute CRS - if epsg is None: - if asset_epsg is None: + if crs is None: + if asset_crs is None: raise ValueError( f"Cannot pick a common CRS, since asset {id!r} of item " f"{item_i} {item['id']!r} does not have one.\n\n" - "Please specify a CRS with the `epsg=` argument." + "Please specify a CRS with the `crs=` argument." ) - if out_epsg is None: - out_epsg = asset_epsg - elif out_epsg != asset_epsg: + if out_crs is None: + out_crs = asset_crs + elif out_crs != asset_crs: raise ValueError( f"Cannot pick a common CRS, since assets have multiple CRSs: asset {id!r} of item " - f"{item_i} {item['id']!r} is in EPSG:{asset_epsg}, " - f"but assets before it were in EPSG:{out_epsg}.\n\n" - "Please specify a CRS with the `epsg=` argument." + f"{item_i} {item['id']!r} is in {asset_crs}, " + f"but assets before it were in {out_crs}.\n\n" + "Please specify a CRS with the `crs=` argument." ) - assert isinstance(out_epsg, int), f"`out_epsg` not found. {out_epsg=}" + assert isinstance(out_crs, int), f"`out_crs` not found. {out_crs=}" # ^ because if it was None initially, and we didn't error out in the above check, it's now always set if bounds_latlon is not None and out_bounds is None: out_bounds = bounds = geom_utils.reproject_bounds( - bounds_latlon, 4326, out_epsg + bounds_latlon, 4326, out_crs ) # NOTE: we wait to reproject until now, so we can use the inferred CRS @@ -219,7 +219,7 @@ def prepare_items( # just reproject that and use it if ( asset_bbox is not None - and asset_epsg is not None + and asset_crs is not None and asset_transform == item_transform and asset_shape == item_shape # TODO this still misses the case where the asset overrides bbox, but not transform/shape. @@ -229,7 +229,7 @@ def prepare_items( # then transform from item, then latlon bbox from item ): asset_bbox_proj = geom_utils.reproject_bounds( - asset_bbox, asset_epsg, out_epsg + asset_bbox, asset_crs, out_crs ) # If there's no bbox (or asset-level metadata is more accurate), compute one from the shape and geotrans @@ -237,15 +237,15 @@ def prepare_items( if ( asset_transform is not None and asset_shape is not None - and asset_epsg is not None + and asset_crs is not None ): asset_affine = affine.Affine(*asset_transform[:6]) asset_bbox_proj = geom_utils.bounds_from_affine( asset_affine, asset_shape[0], asset_shape[1], - asset_epsg, - out_epsg, + asset_crs, + out_crs, ) # There's no bbox, nor shape and transform. The only info we have is `item.bbox` in lat-lon. @@ -258,7 +258,7 @@ def prepare_items( else: # TODO handle error asset_bbox_proj = geom_utils.reproject_bounds( - bbox_lonlat, 4326, out_epsg + bbox_lonlat, 4326, out_crs ) item_bbox_proj = asset_bbox_proj # ^ so we can reuse for other assets @@ -268,9 +268,9 @@ def prepare_items( # Auto-compute resolutions if resolution is None: # Prefer computing resolutions from a geotrans, if it exists - if asset_transform is not None and asset_epsg is not None: + if asset_transform is not None and asset_crs is not None: asset_affine = asset_affine or affine.Affine(*asset_transform[:6]) - if asset_epsg == out_epsg: + if asset_crs == out_crs: # Fastpath-ish when asset is already in the output CRS: # pull directly from geotrans coefficients if not asset_affine.is_rectilinear: @@ -291,7 +291,7 @@ def prepare_items( ) transformer = geom_utils.cached_transformer( - asset_epsg, out_epsg, always_xy=True + asset_crs, out_crs, always_xy=True ) out_px_corner_xs, out_px_corner_ys = transformer.transform( px_corner_xs, px_corner_ys, errcheck=True @@ -308,11 +308,11 @@ def prepare_items( f"since asset {id!r} on item {item_i} {item['id']!r} " f"doesn't provide enough metadata to determine its native resolution.\n" f"We'd need at least one of (in order of preference):\n" - f"- The `proj:transform` and `proj:epsg` fields set on the asset, or on the item\n" + f"- The `proj:transform` and `proj:epsg`/`proj:code` fields set on the asset, or on the item\n" f"- The `proj:shape` and one of `proj:bbox` or `bbox` fields set on the asset, " "or on the item\n\n" "Please specify the `resolution=` argument to set the output resolution manually. " - f"(Remember that resolution must be in the units of your CRS (http://epsg.io/{out_epsg})" + f"(Remember that resolution must be in the units of your CRS (http://epsg.io/{out_crs})" "---not necessarily meters." ) @@ -371,12 +371,12 @@ def prepare_items( # At this point, everything has been set (or there was as error) assert out_bounds, f"{out_bounds=}" assert out_resolutions_xy is not None, f"{out_resolutions_xy=}" - assert out_epsg is not None, f"{out_epsg=}" + assert out_crs is not None, f"{out_crs=}" if snap_bounds: out_bounds = geom_utils.snapped_bounds(out_bounds, out_resolutions_xy) spec = RasterSpec( - epsg=out_epsg, + crs=out_crs, bounds=out_bounds, resolutions_xy=out_resolutions_xy, ) @@ -556,14 +556,14 @@ def to_coords( ) ) - # Add `epsg` last in case it's also a field in properties; our data model assumes it's a coordinate - coords["epsg"] = spec.epsg + # Add `crs` last in case it's also a field in properties; our data model assumes it's a coordinate + coords["crs"] = spec.crs return coords, dims def to_attrs(spec: RasterSpec) -> Dict[str, Any]: - attrs = {"spec": spec, "crs": f"epsg:{spec.epsg}", "transform": spec.transform} + attrs = {"spec": spec, "crs": spec.crs, "transform": spec.transform} resolutions = spec.resolutions_xy if resolutions[0] == resolutions[1]: diff --git a/stackstac/raster_spec.py b/stackstac/raster_spec.py index 2a0f190..5d8e6cf 100644 --- a/stackstac/raster_spec.py +++ b/stackstac/raster_spec.py @@ -15,7 +15,7 @@ class RasterSpec: Spatial parameters defining the grid for a raster. """ - epsg: int + _crs: str bounds: Bbox resolutions_xy: Resolutions @@ -28,6 +28,17 @@ def __post_init__(self): assert minx < maxx, f"Invalid bounds: {minx=} >= {maxx=}" assert miny < maxy, f"Invalid bounds: {miny=} >= {maxy=}" + @property + def crs(self) -> str: + return self._crs + + @crs.setter + def crs(self, v: Union[str, int]) -> None: + if isinstance(v, int): + v = f"EPSG:{v}" + + self._crs = v + @cached_property def transform(self) -> affine.Affine: return affine.Affine( @@ -56,7 +67,7 @@ def shape(self) -> Tuple[int, int]: def vrt_params(self) -> dict: height, width = self.shape return { - "crs": self.epsg, + "crs": self.crs, "transform": self.transform, "height": height, "width": width, diff --git a/stackstac/rio_reader.py b/stackstac/rio_reader.py index 13b31fb..229f4bf 100644 --- a/stackstac/rio_reader.py +++ b/stackstac/rio_reader.py @@ -343,7 +343,7 @@ def _open(self) -> ThreadsafeRioDataset: # Only make a VRT if the dataset doesn't match the spatial spec we want if self.spec.vrt_params != { - "crs": ds.crs.to_epsg(), + "crs": ds.crs.to_string(), "transform": ds.transform, "height": ds.height, "width": ds.width, diff --git a/stackstac/show.py b/stackstac/show.py index 7fe376d..ad2772d 100644 --- a/stackstac/show.py +++ b/stackstac/show.py @@ -449,7 +449,7 @@ def register( """ loop = ensure_server() - geom_utils.array_epsg(arr) # just for the error + geom_utils.array_crs(arr) # just for the error if arr.ndim == 2: assert set(arr.dims) == {"x", "y"} arr = arr.expand_dims("band") @@ -656,7 +656,7 @@ def add_to_map( ---------- arr: `~xarray.DataArray` to visualize. Must have ``x`` and ``y``, and optionally ``band`` dims, - and the ``epsg`` coordinate set. + and the ``crs`` coordinate set. ``arr`` must have 1-3 bands. Single-band data can be colormapped; multi-band data will be displayed as RGB. For 2-band arrays, the first band will be duplicated into the third band's spot, @@ -740,7 +740,7 @@ def show( ---------- arr: `~xarray.DataArray` to visualize. Must have ``x`` and ``y``, and optionally ``band`` dims, - and the ``epsg`` coordinate set. + and the ``crs`` coordinate set. ``arr`` must have 1-3 bands. Single-band data can be colormapped; multi-band data will be displayed as RGB. For 2-band arrays, the first band will be duplicated into the third band's spot, @@ -782,14 +782,14 @@ def show( """ map_ = ipyleaflet.Map(**map_kwargs) if center is None or zoom is None: - west, south, east, north = geom_utils.array_bounds(arr, to_epsg=4326) + west, south, east, north = geom_utils.array_bounds(arr, to_crs=4326) if center is None: center = south + (north - south) / 2, west + (east - west) / 2 if zoom is None: west_m, south_m, east_m, north_m = geom_utils.reproject_bounds( - (west, south, east, north), from_epsg=4326, to_epsg=3857 + (west, south, east, north), from_crs=4326, to_crs=3857 ) size_m = max(east_m - west_m, north_m - south_m) target_map_size_px = 800 diff --git a/stackstac/stac_types.py b/stackstac/stac_types.py index 738c4c5..b748fba 100644 --- a/stackstac/stac_types.py +++ b/stackstac/stac_types.py @@ -102,6 +102,7 @@ class EOBand(TypedDict, total=False): { "datetime": Optional[str], "proj:epsg": Optional[int], + "proj:code": Optional[str], "proj:bbox": Tuple[float, float, float, float], "proj:shape": Tuple[int, int], "proj:transform": Union[ diff --git a/stackstac/stack.py b/stackstac/stack.py index dc85fd9..a9f83ae 100644 --- a/stackstac/stack.py +++ b/stackstac/stack.py @@ -30,7 +30,7 @@ def stack( assets: Optional[Union[List[str], AbstractSet[str]]] = frozenset( ["image/tiff", "image/x.geotiff", "image/vnd.stac.geotiff", "image/jp2"] ), - epsg: Optional[int] = None, + crs: Optional[Union[str, int]] = None, resolution: Optional[Union[IntFloat, Resolutions]] = None, bounds: Optional[Bbox] = None, bounds_latlon: Optional[Bbox] = None, @@ -59,7 +59,7 @@ def stack( We'll try to choose the output coordinate reference system, resolution, and bounds based on the metadata in the STAC items. However, if not all items have the necessary metadata, or aren't in the same coordinate reference system, you'll have specify these - yourself---``epsg`` and ``resolution`` are the two parameters you'll set most often. + yourself---``crs`` and ``resolution`` are the two parameters you'll set most often. Examples -------- @@ -71,7 +71,7 @@ def stack( >>> xr_stack = stackstac.stack(items) >>> >>> # Reproject to 100-meter resolution in web mercator - >>> xr_stack = stackstac.stack(items, epsg=3857, resolution=100) + >>> xr_stack = stackstac.stack(items, crs=3857, resolution=100) >>> >>> # Only use specific asset IDs >>> xr_stack = stackstac.stack(items, assets=["B01", "B03", "B02"]) @@ -122,10 +122,10 @@ def stack( are not yet supported. .. _MT: https://github.com/radiantearth/stac-spec/blob/master/best-practices.md#common-media-types-in-stac - epsg: - Reproject into this coordinate reference system, as given by an `EPSG code `_. + crs: + Reproject into this coordinate reference system, as given by an CRS string (authority:code) or `EPSG code `_. If None (default), uses whatever CRS is set on all the items. In this case, all Items/Assets - must have the ``proj:epsg`` field, and it must be the same value for all of them. + must have the ``proj:epsg`` or ``proj:code`` field, and it must be the same value for all of them. resolution: Output resolution. Careful: this must be given in the output CRS's units! For example, with ``epsg=4326`` (meaning lat-lon), the units are degrees of @@ -146,7 +146,7 @@ def stack( bounds: Output spatial bounding-box, as a tuple of ``(min_x, min_y, max_x, max_y)``. This defines the ``(west, south, east, north)`` rectangle the output array will cover. - Values must be in the same coordinate reference system as ``epsg``. + Values must be in the same coordinate reference system as ``crs``. If None (default), the bounding box of all the input items is automatically used. (This only requires the ``bbox`` field to be set on each Item, which is a required @@ -218,7 +218,7 @@ def stack( xy_coords: Whether to add geospatial coordinate labels for the ``x`` and ``y`` dimensions of the DataArray, allowing for spatial indexing. The coordinates will be in the coordinate reference system given - by ``epsg`` + by ``crs`` If ``"topleft"`` (default), the coordinates are for each pixel's upper-left corner, following raster conventions. @@ -290,7 +290,7 @@ def stack( asset_table, spec, asset_ids, plain_items = prepare_items( plain_items, assets=assets, - epsg=epsg, + crs=crs, resolution=resolution, bounds=bounds, bounds_latlon=bounds_latlon,