Merge pull request #6 from charlottecrevier/extract-point-value

Added the read of a single point coordinate to rasterio examples

Author: charlottecrevier

docs/example-cogs.md

@@ -42,6 +42,18 @@ For more details on the rasterio library, see the [rasterio] documentation.
 ```
 <!-- END: Read a portion of a cog using rasterio -->
 
+<!-- START: Read pixel with coord using rasterio -->
+::: how-to-guides.rasterio-point-example
+    options:
+        show_source: false
+        members: no
+        show_root_toc_entry: false # To remove the name of the file in the TOC
+
+``` py linenums="1" hl_lines="34-38"
+--8<-- "how-to-guides/rasterio-point-example.py:code"
+```
+<!-- END: Read pixel with coord cog using rasterio -->
+
 [rasterio installation]: https://rasterio.readthedocs.io/en/stable/installation.html
 [rasterio]: https://rasterio.readthedocs.io/en/latest/quickstart.html
 [Interacting with CCMEO STAC API]: pystac-client.md

how-to-guides/rasterio-point-example.py

@@ -0,0 +1,64 @@
+"""
+## Read a portion of a remote COG - Point
+
+For this example, you will need to also install [shapely](https://shapely.readthedocs.io/en/stable/installation.html):   
+```bash
+pip install shapely
+```
+
+In this code you will : 
+
+- Query a STAC API with pystac-client;
+- Read values of remote COG based on coordinates with the [sample()](https://rasterio.readthedocs.io/en/stable/api/rasterio.io.html#rasterio.io.DatasetReader.sample) functionality;
+
+!!! info
+    This specific example uses the collection **mrdem-30** from CCMEO's datacube
+
+!!! Tip
+    To perform the same request using gdal, refere to the [gdallocationinfo](https://gdal.org/en/stable/programs/gdallocationinfo.html) utility 
+
+"""
+# --8<-- [start:code]
+import pystac_client
+import rasterio 
+import pyproj
+from shapely.geometry import Point
+from shapely.ops import transform
+
+
+# Define the coordinate of the point of interest in EPSG:4326 (for the request to the STAC API)
+lonlat = Point(-104.898838, 69.232434)
+
+# Modify the projection of the point (for the COG extraction)
+proj = pyproj.Transformer.from_crs(4326, 3979, always_xy=True)
+point =  transform(proj.transform, lonlat)
+x1, y1 = point.coords.xy
+
+
+# Link to ccmeo datacube stac-api
+stac_root = "https://datacube.services.geo.ca/stac/api"
+catalog = pystac_client.Client.open(stac_root)
+
+search = catalog.search(
+	collections=['mrdem-30'], 
+    intersects=lonlat,
+	) 
+
+# Get the link to the data asset for mrdem-30 dtm and dsm
+links = {}
+for page in search.pages():
+	for item in page:
+		links['dtm'] = item.assets['dtm'].href
+		links['dsm'] = item.assets['dsm'].href
+
+# Read value from coordinates
+for asset, href in links.items():
+    with rasterio.open(href) as src:
+        for val in src.sample([(x1, y1)]): 
+            print(f'{asset} : {val}')
+   
+# --8<-- [end:code]
+
+
+
+

how-to-guides/rasterio-window-example.py

@@ -1,5 +1,5 @@
 """
-## Read a portion of a remote COG
+## Read a portion of a remote COG - Bounding box
 
 For this example, you will need to also install [shapely](https://shapely.readthedocs.io/en/stable/installation.html):   
 ```bash
@@ -10,7 +10,7 @@
 
 - Query a STAC API with pystac-client to get link to a COG;
 - Read a portion of a remote COG based on an AOI with the window functionality;
-- Optionally write the portion locally inside a .tif file
+- Write the portion locally inside a .tif file (Optional)
 
 !!! info
     This specific example uses the collection **mrdem-30** from CCMEO's datacube
@@ -73,7 +73,7 @@
 
 # Perform analysis ...
 
-# Write the output array to a tiff file
+# (Optional) Write the output array to a tiff file
 output_tiff = r"path/to/output.tif"
 with rasterio.open(output_tiff, 'w', **metadata) as dst:
     dst.write(rst, 1)