diff --git a/docs/vector-segmentation-masks.md b/docs/vector-segmentation-masks.md index 3f04de1..4f6e1d3 100644 --- a/docs/vector-segmentation-masks.md +++ b/docs/vector-segmentation-masks.md @@ -84,18 +84,18 @@ dp_rioxarray ``` The Sentinel-1 image from Planetary Computer comes in longitude/latitude 🌐 -geographic coordinates by default (EPSG:4326). To make the pixels more equal 🔲 +geographic coordinates by default (OGC:CRS84). To make the pixels more equal 🔲 area, we can project it to a 🌏 local projected coordinate system instead. ```{code-cell} def reproject_to_local_utm(dataarray: xr.DataArray, resolution: float=100.0) -> xr.DataArray: """ - Reproject an xarray.DataArray grid from EPSG;4326 to a local UTM coordinate + Reproject an xarray.DataArray grid from OGC:CRS84 to a local UTM coordinate reference system. """ # Estimate UTM coordinate reference from a single pixel pixel = dataarray.isel(y=slice(0, 1), x=slice(0,1)) - new_crs = dataarray.rio.reproject(dst_crs="EPSG:4326").rio.estimate_utm_crs() + new_crs = dataarray.rio.reproject(dst_crs="OGC:CRS84").rio.estimate_utm_crs() return dataarray.rio.reproject(dst_crs=new_crs, resolution=resolution) ``` @@ -113,6 +113,16 @@ matter much when we're looking at spatial resolutions over several 10s of metres 🙂. ``` +```{hint} +For those wondering what `OGC:CRS84` is, it is the longitude/latitude version +of [`EPSG:4326`](https://epsg.io/4326) 🌐 (latitude/longitude). I.e., it's a +matter of axis order, with `OGC:CRS84` being x/y and `EPSG:4326` being y/x. + +🔖 References: +- https://gis.stackexchange.com/questions/54073/what-is-crs84-projection +- https://github.com/opengeospatial/geoparquet/issues/52 +``` + ### Transform and visualize raster data 🔎 Let's visualize 👀 the Sentinel-1 image, but before that, we'll transform 🔄 @@ -201,10 +211,11 @@ geodataframe.dropna(axis="columns") Cool, and we can also visualize the polygons 🔷 on a 2D map. To align the coordinates with the 🛰️ Sentinel-1 image above, we'll first use {py:meth}`geopandas.GeoDataFrame.to_crs` to reproject the vector from 🌐 -EPSG:4326 (longitude/latitude) to 🌏 EPSG:32648 (UTM Zone 48N). +EPSG:9707 (WGS 84 + EGM96 height, latitude/longitude) to 🌏 EPSG:32648 (UTM +Zone 48N). ```{code-cell} -print(f"Original bounds in EPSG:4326:\n{geodataframe.bounds}") +print(f"Original bounds in EPSG:9707:\n{geodataframe.bounds}") gdf = geodataframe.to_crs(crs="EPSG:32648") print(f"New bounds in EPSG:32648:\n{gdf.bounds}") ``` diff --git a/docs/walkthrough.md b/docs/walkthrough.md index e99635d..e2380d3 100644 --- a/docs/walkthrough.md +++ b/docs/walkthrough.md @@ -84,14 +84,14 @@ This is how the Sentinel-2 image looks like over Singapore on 15 Jan 2022. ![Sentinel-2 image over Singapore on 20220115](https://planetarycomputer.microsoft.com/api/data/v1/item/preview.png?collection=sentinel-2-l2a&item=S2A_MSIL2A_20220115T032101_R118_T48NUG_20220115T170435&assets=visual&asset_bidx=visual%7C1%2C2%2C3&nodata=0) -## 1️⃣ Construct [DataPipe](https://github.com/pytorch/data/tree/v0.3.0#what-are-datapipes) 📡 +## 1️⃣ Construct [DataPipe](https://github.com/pytorch/data/tree/v0.4.0#what-are-datapipes) 📡 A torch `DataPipe` is a way of composing data (rather than inheriting data). Yes, I don't know what it really means either, so here's some extra reading. 🔖 References: - https://pytorch.org/blog/pytorch-1.11-released/#introducing-torchdata -- https://github.com/pytorch/data/tree/v0.3.0#what-are-datapipes +- https://github.com/pytorch/data/tree/v0.4.0#what-are-datapipes - https://realpython.com/inheritance-composition-python ### Create an Iterable 📏 diff --git a/zen3geo/datapipes/datashader.py b/zen3geo/datapipes/datashader.py index 8143b9e..5aa8b52 100644 --- a/zen3geo/datapipes/datashader.py +++ b/zen3geo/datapipes/datashader.py @@ -84,9 +84,9 @@ class DatashaderRasterizerIterDataPipe(IterDataPipe): AttributeError If either the canvas in ``source_datapipe`` or vector geometry in ``vector_datapipe`` is missing a ``.crs`` attribute. Please set the - coordinate reference system (e.g. using ``canvas.crs = 'EPSG:4326'`` + coordinate reference system (e.g. using ``canvas.crs = 'OGC:CRS84'`` for the :py:class:`datashader.Canvas` input or - ``vector = vector.set_crs(epsg=4326)`` for the + ``vector = vector.set_crs(crs='OGC:CRS84')`` for the :py:class:`geopandas.GeoSeries` or :py:class:`geopandas.GeoDataFrame` input) before passing them into the datapipe. @@ -113,7 +113,7 @@ class DatashaderRasterizerIterDataPipe(IterDataPipe): ... "https://github.com/geopandas/pyogrio/raw/v0.4.0/pyogrio/tests/fixtures/test_gpkg_nulls.gpkg", ... read_geometry=True, ... ) - >>> assert geodataframe.crs == "EPSG:4326" # longitude/latitude coords + >>> assert geodataframe.crs == "EPSG:4326" # latitude/longitude coords >>> dp_vector = IterableWrapper(iterable=[geodataframe]) ... >>> # Setup blank raster canvas where we will burn vector geometries onto @@ -196,7 +196,7 @@ def __iter__(self) -> Iterator[xr.DataArray]: raise AttributeError( "Missing crs information for datashader.Canvas with " f"x_range: {canvas.x_range} and y_range: {canvas.y_range}. " - "Please set crs using e.g. `canvas.crs = 'EPSG:4326'`." + "Please set crs using e.g. `canvas.crs = 'OGC:CRS84'`." ) # Reproject vector geometries to coordinate reference system @@ -208,7 +208,7 @@ def __iter__(self) -> Iterator[xr.DataArray]: raise AttributeError( f"Missing crs information for input {vector.__class__} object " f"with the following bounds: \n {vector.bounds} \n" - f"Please set crs using e.g. `vector = vector.set_crs(epsg=4326)`." + f"Please set crs using e.g. `vector = vector.set_crs(crs='OGC:CRS84')`." ) from e # Convert vector to spatialpandas format to allow datashader's diff --git a/zen3geo/tests/test_datapipes_datashader.py b/zen3geo/tests/test_datapipes_datashader.py index 6947a59..bac6542 100644 --- a/zen3geo/tests/test_datapipes_datashader.py +++ b/zen3geo/tests/test_datapipes_datashader.py @@ -19,7 +19,7 @@ def fixture_canvas(): canvas = datashader.Canvas( plot_width=14, plot_height=10, x_range=(1, 8), y_range=(0, 5) ) - canvas.crs = "EPSG:4326" + canvas.crs = "OGC:CRS84" return canvas @@ -38,7 +38,7 @@ def fixture_geodataframe(): shapely.geometry.Polygon([(6, 5), (3.5, 2.5), (6, 0), (6, 2.5), (5, 2.5)]), ] geodataframe = gpd.GeoDataFrame(data={"geometry": geometries}) - geodataframe = geodataframe.set_crs(epsg=4326) + geodataframe = geodataframe.set_crs(crs="OGC:CRS84") return geodataframe @@ -99,7 +99,7 @@ def test_datashader_rasterize_vector_geometry(canvas, geodataframe, geom_type, s assert dataarray.data.sum() == sum_val assert dataarray.dims == ("y", "x") - assert dataarray.rio.crs == "EPSG:4326" + assert dataarray.rio.crs == "OGC:CRS84" assert dataarray.rio.shape == (10, 14) assert dataarray.rio.transform().e == -0.5 @@ -161,7 +161,7 @@ def test_datashader_rasterize_vector_geometrycollection(canvas, geodataframe): # Merge points, lines and polygons into a single GeometryCollection geocollection = gpd.GeoSeries(data=geodataframe.unary_union) - geocollection = geocollection.set_crs(epsg=4326) + geocollection = geocollection.set_crs(crs="OGC:CRS84") dp = IterableWrapper(iterable=[canvas]) dp_vector = IterableWrapper(iterable=[geocollection])