11
2- # From : https://geog-312.gishub.org/book/geospatial/rasterio.html#reading-raster-data
3- import rasterio
4- import rasterio .plot
2+ # # From : https://geog-312.gishub.org/book/geospatial/rasterio.html#reading-raster-data
3+ # import rasterio
4+ # import rasterio.plot
55
6- import geopandas as gpd
7- import numpy as np
8- import matplotlib .pyplot as plt
6+ # import geopandas as gpd
7+ # import numpy as np
8+ # import matplotlib.pyplot as plt
99
1010
11- raster_path = (
12- "https://github.com/opengeos/datasets/releases/download/raster/dem_90m.tif"
13- )
14- src = rasterio .open (raster_path )
11+ # raster_path = (
12+ # "https://github.com/opengeos/datasets/releases/download/raster/dem_90m.tif"
13+ # )
14+ # src = rasterio.open(raster_path)
1515
16- print (src .name )
16+ # print(src.name)
1717
18- # rasterio.plot.show((src, 1))
19- rasterio .plot .show (src )
18+ # # rasterio.plot.show((src, 1))
19+ # rasterio.plot.show(src)
2020
21- fig , ax = plt .subplots (figsize = (8 , 8 ))
22- rasterio .plot .show (src , cmap = "terrain" , ax = ax , title = "Digital Elevation Model (DEM)" )
23- plt .show ()
21+ # fig, ax = plt.subplots(figsize=(8, 8))
22+ # rasterio.plot.show(src, cmap="terrain", ax=ax, title="Digital Elevation Model (DEM)")
23+ # plt.show()
24+
25+ from pystac_client import Client
26+ from odc .stac import load
27+ import odc .geo
28+
29+ # use publically available stac link such as
30+ client = Client .open ("https://earth-search.aws.element84.com/v1" )
31+
32+ # ID of the collection
33+ collection = "sentinel-2-l2a"
34+
35+ # Geometry of AOI
36+ geometry = {
37+ "coordinates" : [
38+ [
39+ [74.66218437999487 , 19.46556170905807 ],
40+ [74.6629598736763 , 19.466339343697722 ],
41+ [74.6640371158719 , 19.4667885366414 ],
42+ [74.66395296156406 , 19.46614872872264 ],
43+ [74.66376889497042 , 19.466150941501425 ],
44+ [74.66369077563286 , 19.46577508478787 ],
45+ [74.6635865047574 , 19.465278788212864 ],
46+ [74.66282073408365 , 19.46540270444271 ],
47+ [74.66218437999487 , 19.46556170905807 ],
48+ ]
49+ ],
50+ "type" : "Polygon" ,
51+ }
52+
53+ # Complete month
54+ date_YYMM = "2023-01"
55+ # run pystac client search to see available dataset
56+ search = client .search (
57+ collections = [collection ], intersects = geometry , datetime = date_YYMM
58+ )
59+ # spit out data as GeoJSON dictionary
60+ print (search .item_collection_as_dict ())
61+ # loop through each item
62+ for item in search .items_as_dicts ():
63+ print (item )
0 commit comments