Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

✏️ Edit docs to use OGC:CRS84 lon/lat instead of EPSG:4326 #45

Merged
merged 2 commits into from
Aug 28, 2022
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
21 changes: 16 additions & 5 deletions docs/vector-segmentation-masks.md
Original file line number Diff line number Diff line change
Expand Up @@ -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)
```
Expand All @@ -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 🔄
Expand Down Expand Up @@ -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}")
```
Expand Down
4 changes: 2 additions & 2 deletions docs/walkthrough.md
Original file line number Diff line number Diff line change
Expand Up @@ -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 📏
Expand Down
10 changes: 5 additions & 5 deletions zen3geo/datapipes/datashader.py
Original file line number Diff line number Diff line change
Expand Up @@ -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.

Expand All @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down
8 changes: 4 additions & 4 deletions zen3geo/tests/test_datapipes_datashader.py
Original file line number Diff line number Diff line change
Expand Up @@ -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


Expand All @@ -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

Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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])
Expand Down