Skip to content

Commit 2a56b04

Browse files
alpha-beta-soupalpha-beta-soup
committed
update readme
1 parent 18a1d99 commit 2a56b04

File tree

2 files changed

+86
-180
lines changed

2 files changed

+86
-180
lines changed

‎README.md‎

Lines changed: 86 additions & 180 deletions
Original file line numberDiff line numberDiff line change
@@ -138,193 +138,97 @@ tree /home/user/example.pq
138138
└── part.0.parquet
139139
```
140140

141-
For a quick view of your output, you can read Apache Parquet with pandas, and then use h3-pandas and geopandas to convert this into a GeoPackage for visualisation in a desktop GIS, such as QGIS. The Apache Parquet output is indexed by the DGGS column, so it should be ready for association with other data prepared in the same DGGS.
142-
143-
<details>
144-
<summary>For H3 output...</summary>
145-
146-
```python
147-
>>> import pandas as pd
148-
>>> import h3pandas
149-
>>> o = pd.read_parquet('./tests/data/output/9/Sen2_Test')
150-
>>> o
151-
band B02 B03 B04 B05 B06 B07 B08 B8A B11 B12
152-
h3_09
153-
89bb0981003ffff 9 27 16 62 175 197 228 247 102 36
154-
89bb0981007ffff 10 30 17 66 185 212 238 261 113 40
155-
89bb098100bffff 10 26 15 60 169 190 228 241 103 37
156-
89bb098100fffff 11 29 17 66 181 203 243 257 109 39
157-
89bb0981013ffff 8 26 16 58 172 199 220 244 98 34
158-
... ... ... ... ... ... ... ... ... ... ...
159-
89bb0d6eea7ffff 10 18 15 41 106 120 140 146 102 47
160-
89bb0d6eeabffff 12 19 15 39 95 107 125 131 84 39
161-
89bb0d6eeafffff 12 21 17 43 101 115 134 141 111 51
162-
89bb0d6eeb7ffff 10 20 14 45 120 137 160 165 111 48
163-
89bb0d6eebbffff 15 28 20 56 146 166 198 202 108 47
164-
165-
[5656 rows x 10 columns]
166-
>>> o.h3.h3_to_geo_boundary().to_file('~/Downloads/Sen2_Test_h3-9.gpkg', driver='GPKG')
167-
```
168-
</details>
169-
170-
<details>
171-
<summary>For rHEALPix output...</summary>
172-
173-
For rHEALPix DGGS output, you can use [`rHP-Pandas`](https://github.com/manaakiwhenua/rHP-Pandas):
174-
175-
```python
176-
>>> import pandas as pd
177-
>>> import rhppandas
178-
>>> o = pd.read_parquet('./tests/data/output/11/Sen2_Test_rhp')
179-
>>> o
180-
band B02 B03 B04 B05 B06 B07 B08 B8A B11 B12
181-
rhp_11
182-
R88723652267 11 31 16 65 191 217 263 274 99 36
183-
R88723652268 11 30 15 66 190 214 258 269 96 34
184-
R88723652276 11 27 17 66 179 203 240 255 98 36
185-
R88723652277 13 30 19 68 179 204 246 260 108 41
186-
R88723652278 12 29 20 66 176 199 243 255 110 43
187-
... ... ... ... ... ... ... ... ... ... ...
188-
R88727068804 22 39 41 81 151 167 182 203 166 84
189-
R88727068805 22 40 42 81 150 166 185 203 167 85
190-
R88727068806 23 41 43 83 156 175 188 211 164 83
191-
R88727068807 23 41 42 82 154 171 186 207 164 83
192-
R88727068808 22 39 43 80 146 163 177 198 165 83
193-
194-
[223104 rows x 10 columns]
195-
>>> o.rhp.rhp_to_geo_boundary().to_file('~/Downloads/Sen2_Test_rhp-11.gpkg', driver='GPKG')
196-
```
197-
</details>
198-
199-
<details>
200-
<summary>For S2 output...</summary>
201-
202-
For S2 output, use [`s2sphere`](https://pypi.org/project/s2sphere/):
203-
204-
```python
205-
import pandas as pd
206-
import geopandas as gpd
207-
import s2sphere
208-
from shapely.geometry import Polygon
209-
210-
o = pd.read_parquet('./tests/data/output/7/sample_tif_s2')
211-
212-
def s2id_to_polygon(s2_id_hex):
213-
cell_id = s2sphere.CellId.from_token(s2_id_hex)
214-
cell = s2sphere.Cell(cell_id)
215-
vertices = []
216-
for i in range(4):
217-
vertex = cell.get_vertex(i)
218-
# Convert to lat/lon degrees
219-
lat_lng = s2sphere.LatLng.from_point(vertex)
220-
vertices.append((lat_lng.lng().degrees, lat_lng.lat().degrees)) # (lon, lat)
221-
return Polygon(vertices)
222-
223-
o['geometry'] = o.index.map(s2id_to_polygon)
224-
gpd.GeoDataFrame(o, geometry='geometry', crs='EPSG:4326').to_file('./tests/data/output/7/sample_tif_s2.gpkg')
225-
```
226-
</details>
227-
228-
<details>
229-
<summary>For A5 output...</summary>
141+
Output can also be written to GeoParquet (v1.1.0) by including the `-g/--geo` parameter, which accepts:
142+
- `polygon` for cells represented as boundary polygons
143+
- `point` for cells represented as centre points
144+
- `none` for standard Parquet output (not GeoParquet) ← this is the default if `-g/--geo` is not used
230145

231-
For A5 output, use [`a5`](https://pypi.org/project/pya5/):
146+
GeoParquet output is useful if you want to use the spatial representations of the DGGS cells in traditional spatial analysis, or if you merely want to visualise the output.
232147

233-
```python
234-
import pandas as pd
235-
import geopandas as gpd
236-
import a5
237-
from shapely.geometry import Polygon
148+
Below are some ways to read and visualise it.
238149

239-
o = pd.read_parquet('./tests/data/output/7/sample_tif_a5')
240-
o['geometry'] = o.index.map(lambda a: Polygon(a5.cell_to_boundary(a5.hex_to_u64(a))))
150+
### DuckDB
241151

242-
gpd.GeoDataFrame(o, geometry="geometry", crs="EPSG:4326").to_file('tests/data/output/7/sample_tif_a5.gpkg')
152+
```bash
153+
$ duckdb
154+
DuckDB v1.4.1 (Andium) b390a7c376
155+
Enter ".help" for usage hints.
156+
Connected to a transient in-memory database.
157+
Use ".open FILENAME" to reopen on a persistent database.
158+
D INSTALL spatial;
159+
D LOAD spatial;
160+
D SELECT * FROM read_parquet('se_island.pq') LIMIT 7;
161+
┌┌────────┬────────┬────────┬────────────────────────────────────────────────────────────────────────────────┬─────────────┬─────────┐
162+
│ band_1 │ band_2 │ band_3 │ geometry │ s2_19 │ s2_08 │
163+
│ float │ float │ float │ geometry │ varchar │ varchar │
164+
├────────┼────────┼────────┼────────────────────────────────────────────────────────────────────────────────┼─────────────┼─────────┤
165+
│ 0.0 │ 0.0 │ 0.0 │ POLYGON ((-176.17946725380486 -44.33542073938414, -176.17946725380486 -44.33… │ 72b47e01e24 │ 72b47 │
166+
│ 0.0 │ 0.0 │ 0.0 │ POLYGON ((-176.18439390505398 -44.33543749229784, -176.18439390505398 -44.33… │ 72b47e02a14 │ 72b47 │
167+
│ 0.0 │ 0.1 │ 0.1 │ POLYGON ((-176.18550630891403 -44.33547457195554, -176.18550630891403 -44.33… │ 72b47e1d54c │ 72b47 │
168+
│ 0.0 │ 0.0 │ 0.0 │ POLYGON ((-176.17819578278952 -44.33537828938332, -176.17819578278952 -44.33… │ 72b47e01d64 │ 72b47 │
169+
│ 0.1 │ 0.1 │ 0.3 │ POLYGON ((-176.18344039674218 -44.335553297533835, -176.18344039674218 -44.3… │ 72b47e0282c │ 72b47 │
170+
│ 0.0 │ 0.0 │ 0.0 │ POLYGON ((-176.17899045588274 -44.335404822417665, -176.17899045588274 -44.3… │ 72b47e01dfc │ 72b47 │
171+
│ 0.1 │ 0.1 │ 0.3 │ POLYGON ((-176.1832814769592 -44.33554799806149, -176.1832814769592 -44.3356… │ 72b47e02824 │ 72b47 │
172+
└────────┴────────┴────────┴────────────────────────────────────────────────────────────────────────────────┴─────────────┴─────────┘
243173
```
244174
245-
</details>
246-
247-
<details>
248-
<summary>For Geohash output...</summary>
249-
250-
For Geohash output, you can use [`python-geohash`](https://github.com/hkwi/python-geohash) or other similar Geohash library. Example:
251-
252-
```python
253-
import pandas as pd
254-
import geohash
255-
from shapely.geometry import Point, box
256-
import geopandas as gpd
257-
o = pd.read_parquet('./tests/data/output/8/sample_geohash')
258-
259-
260-
def geohash_to_geometry(gh, mode="polygon"):
261-
lat, lon, lat_err, lon_err = geohash.decode_exactly(gh)
262-
263-
if mode == "point":
264-
return Point(lon, lat)
265-
elif mode == "polygon":
266-
return box(lon - lon_err, lat - lat_err, lon + lon_err, lat + lat_err)
267-
else:
268-
raise ValueError("mode must be 'point' or 'polygon'")
269-
270-
o["geometry"] = o.index.map(lambda gh: geohash_to_geometry(gh, mode="polygon"))
271-
272-
'''
273-
band geohash_08 1 2 3 geometry
274-
0 u170f2sq 0 0 0 POLYGON ((4.3238067626953125 52.16686248779297...
275-
1 u170f2sr 0 0 0 POLYGON ((4.3238067626953125 52.16703414916992...
276-
2 u170f2sw 0 0 0 POLYGON ((4.324150085449219 52.16686248779297,...
277-
3 u170f2sx 0 0 0 POLYGON ((4.324150085449219 52.16703414916992,...
278-
4 u170f2sy 0 0 0 POLYGON ((4.324493408203125 52.16686248779297,...
279-
... ... .. .. .. ...
280-
232720 u171mc2g 0 0 0 POLYGON ((4.472808837890625 52.258358001708984...
281-
232721 u171mc2h 0 0 0 POLYGON ((4.471778869628906 52.25852966308594,...
282-
232722 u171mc2k 0 0 0 POLYGON ((4.4721221923828125 52.25852966308594...
283-
232723 u171mc2s 0 0 0 POLYGON ((4.472465515136719 52.25852966308594,...
284-
232724 u171mc2u 0 0 0 POLYGON ((4.472808837890625 52.25852966308594,...
285-
286-
[232725 rows x 5 columns]
287-
'''
288-
289-
gpd.GeoDataFrame(o, geometry="geometry", crs="EPSG:4326").to_file('tests/data/output/8/sample_geohash.gpkg')
175+
### GDAL
176+
177+
```bash
178+
ogrinfo -so -al ./se_island.pq
179+
INFO: Open of `se_island.pq'
180+
using driver `Parquet' successful.
181+
182+
Layer name: se_island
183+
Geometry: Polygon
184+
Feature Count: 18390
185+
Extent: (-176.185824, -44.356933) - (-176.159915, -44.335364)
186+
Layer SRS WKT:
187+
GEOGCRS["WGS 84",
188+
ENSEMBLE["World Geodetic System 1984 ensemble",
189+
MEMBER["World Geodetic System 1984 (Transit)"],
190+
MEMBER["World Geodetic System 1984 (G730)"],
191+
MEMBER["World Geodetic System 1984 (G873)"],
192+
MEMBER["World Geodetic System 1984 (G1150)"],
193+
MEMBER["World Geodetic System 1984 (G1674)"],
194+
MEMBER["World Geodetic System 1984 (G1762)"],
195+
MEMBER["World Geodetic System 1984 (G2139)"],
196+
MEMBER["World Geodetic System 1984 (G2296)"],
197+
ELLIPSOID["WGS 84",6378137,298.257223563,
198+
LENGTHUNIT["metre",1]],
199+
ENSEMBLEACCURACY[2.0]],
200+
PRIMEM["Greenwich",0,
201+
ANGLEUNIT["degree",0.0174532925199433]],
202+
CS[ellipsoidal,2],
203+
AXIS["geodetic latitude (Lat)",north,
204+
ORDER[1],
205+
ANGLEUNIT["degree",0.0174532925199433]],
206+
AXIS["geodetic longitude (Lon)",east,
207+
ORDER[2],
208+
ANGLEUNIT["degree",0.0174532925199433]],
209+
USAGE[
210+
SCOPE["Horizontal component of 3D system."],
211+
AREA["World."],
212+
BBOX[-90,-180,90,180]],
213+
ID["EPSG",4326]]
214+
Data axis to CRS axis mapping: 2,1
215+
Geometry Column = geometry
216+
band_1: Real(Float32) (0.0)
217+
band_2: Real(Float32) (0.0)
218+
band_3: Real(Float32) (0.0)
219+
s2_19: String (0.0)
220+
s2_08: String (0.0)
290221
```
291-
</details>
292-
293-
<details>
294-
<summary>For Maidenhead output...</summary>
295-
296-
For Maidenhead output, you can use [`maidenhead`](https://https://github.com/space-physics/maidenhead) or other similar Maidenhead library. Example:
297-
298-
```python
299-
import pandas as pd
300-
import maidenhead
301-
from shapely.geometry import shape
302-
import geopandas as gpd
303-
o = pd.read_parquet('./tests/data/output/5/sample_maidenhead.pq')
304-
305-
o['geometry'] = o.index.map(lambda mh: shape(maidenhead.to_geoJSONObject(mh, center=True)['features'][1]['geometry']))
306-
307-
'''
308-
band 1 2 3 geometry
309-
maidenhead_5
310-
JO22de80UB 0 0 0 POLYGON ((4.323611111111111 52.16684027777778,...
311-
JO22de80UC 0 0 0 POLYGON ((4.323611111111111 52.16701388888889,...
312-
JO22de80UD 0 0 0 POLYGON ((4.323611111111111 52.1671875, 4.3239...
313-
JO22de80UE 0 0 0 POLYGON ((4.323611111111111 52.16736111111111,...
314-
JO22de80UF 0 0 0 POLYGON ((4.323611111111111 52.16753472222222,...
315-
... .. .. .. ...
316-
JO22fg62PB 0 0 0 POLYGON ((4.471875 52.25850694444444, 4.472222...
317-
JO22fg62QA 0 0 0 POLYGON ((4.472222222222222 52.25833333333333,...
318-
JO22fg62QB 0 0 0 POLYGON ((4.472222222222222 52.25850694444444,...
319-
JO22fg62RA 0 0 0 POLYGON ((4.472569444444445 52.25833333333333,...
320-
JO22fg62RB 0 0 0 POLYGON ((4.472569444444445 52.25850694444444,...
321-
322-
[227470 rows x 4 columns]
323-
'''
324-
325-
gpd.GeoDataFrame(o, geometry="geometry", crs="EPSG:4326").to_file('tests/data/output/5/sample_maidenhead.gpkg')
222+
223+
### QGIS
224+
225+
```bash
226+
qgis sample.pq
326227
```
327-
</details>
228+
229+
With some styling applied:
230+
231+
![Example output shown in QGIS](image.png)
328232
329233
## Installation
330234
@@ -391,6 +295,8 @@ Two sample files have been uploaded to an S3 bucket with `s3:GetObject` public p
391295
392296
You may use these for experimentation. However you can also use local files too, which will be faster. A good, small (5 MB) sample image is available [here](https://github.com/mommermi/geotiff_sample).
393297
298+
A small test file is also available at [`tests/data/se-island.tif`] (tests/data/se-island.tif).
299+
394300
## Example commands
395301
396302
```bash
@@ -402,7 +308,7 @@ raster2dggs rhp --resolution 11 -d 0 s3://raster2dggs-test-data/Sen2_Test.tif ./
402308
```
403309
404310
```bash
405-
raster2dggs h3 --resolution 13 --compression zstd --resampling nearest -a median -d 1 -u 2 s3://raster2dggs-test-data/TestDEM.tif ./tests/data/output/13/TestDEM
311+
raster2dggs h3 --resolution 13 --compression zstd --resampling nearest -a median -d 1 -u 2 --geo polygon s3://raster2dggs-test-data/TestDEM.tif ./tests/data/output/13/TestDEM
406312
```
407313
408314
## Citation

‎image.png‎

747 KB
Loading

0 commit comments

Comments
 (0)