Skip to content

Commit 87c19cb

Browse files
alpha-beta-soupalpha-beta-soup
committed
update readme
1 parent bd60d48 commit 87c19cb

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
@@ -136,193 +136,97 @@ tree /home/user/example.pq
136136
└── part.0.parquet
137137
```
138138

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

229-
For A5 output, use [`a5`](https://pypi.org/project/pya5/):
144+
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.
230145

231-
```python
232-
import pandas as pd
233-
import geopandas as gpd
234-
import a5
235-
from shapely.geometry import Polygon
146+
Below are some ways to read and visualise it.
236147

237-
o = pd.read_parquet('./tests/data/output/7/sample_tif_a5')
238-
o['geometry'] = o.index.map(lambda a: Polygon(a5.cell_to_boundary(a5.hex_to_u64(a))))
148+
### DuckDB
239149

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

‎image.png‎

747 KB
Loading

0 commit comments

Comments
 (0)