Skip to content

Commit aa35aea

Browse files
brendancolbrendancol
committed
Button up reproject: README, benchmarks, write_geotiff WKT fix (#1045)
README Reproject section updated: - All 11 projection families listed (added oblique stereographic) - Full pipeline benchmark table (read+reproject+write, all backends) - Datum shift coverage (14 grids, 10 Helmert fallbacks) - Vertical datum support (EGM96/EGM2008, integrated into reproject) - ITRF time-dependent frame transforms - pyproj usage documented (metadata only, Numba does the math) - Merge performance table updated benchmarks/reproject_benchmark.md: - 254-line comprehensive benchmark document with 6 sections - Full pipeline: NumPy 2.7s, CuPy 348ms, rioxarray 418ms - 13 projections tested for accuracy vs pyproj - Datum, vertical, ITRF coverage documented - All numbers from actual benchmark runs write_geotiff WKT fix: - reproject() stores CRS as WKT string in attrs['crs'] - write_geotiff assumed integer EPSG code, crashed with TypeError - Added isinstance check to parse WKT via _wkt_to_epsg()
1 parent e279904 commit aa35aea

File tree

3 files changed

+292
-16
lines changed

3 files changed

+292
-16
lines changed

‎README.md‎

Lines changed: 33 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -212,43 +212,60 @@ write_vrt('mosaic.vrt', ['tile1.tif', 'tile2.tif']) # generate VRT
212212

213213
| Name | Description | Source | NumPy xr.DataArray | Dask xr.DataArray | CuPy GPU xr.DataArray | Dask GPU xr.DataArray |
214214
|:----------:|:------------|:------:|:----------------------:|:--------------------:|:-------------------:|:------:|
215-
| [Reproject](xrspatial/reproject/__init__.py) | Reprojects a raster to a new CRS with Numba JIT / CUDA coordinate transforms and resampling | Standard (inverse mapping) | ✅️ | ✅️ | ✅️ | ✅️ |
215+
| [Reproject](xrspatial/reproject/__init__.py) | Reprojects a raster to a new CRS with Numba JIT / CUDA coordinate transforms and resampling. Supports vertical datums (EGM96, EGM2008) and horizontal datum shifts (NAD27, OSGB36, etc.) | Standard (inverse mapping) | ✅️ | ✅️ | ✅️ | ✅️ |
216216
| [Merge](xrspatial/reproject/__init__.py) | Merges multiple rasters into a single mosaic with configurable overlap strategy | Standard (mosaic) | ✅️ | ✅️ | 🔄 | 🔄 |
217217

218-
Built-in Numba JIT and CUDA projection kernels bypass pyproj for common CRS pairs. Other CRS pairs fall back to pyproj automatically.
218+
Built-in Numba JIT and CUDA projection kernels bypass pyproj for per-pixel coordinate transforms. pyproj is used only for CRS metadata parsing (~1ms, once per call) and output grid boundary estimation (~500 control points, once per call). Any CRS pair without a built-in kernel falls back to pyproj automatically.
219219

220220
| Projection | EPSG examples | CPU Numba | CUDA GPU |
221221
|:-----------|:-------------|:---------:|:--------:|
222222
| Web Mercator | 3857 | ✅️ | ✅️ |
223223
| UTM / Transverse Mercator | 326xx, 327xx, State Plane | ✅️ | ✅️ |
224224
| Ellipsoidal Mercator | 3395 | ✅️ | ✅️ |
225-
| Lambert Conformal Conic | 2154, State Plane | ✅️ | ✅️ |
225+
| Lambert Conformal Conic | 2154, 2229, State Plane | ✅️ | ✅️ |
226226
| Albers Equal Area | 5070 | ✅️ | ✅️ |
227227
| Cylindrical Equal Area | 6933 | ✅️ | ✅️ |
228228
| Sinusoidal | MODIS grids | ✅️ | ✅️ |
229229
| Lambert Azimuthal Equal Area | 3035, 6931, 6932 | ✅️ | ✅️ |
230230
| Polar Stereographic | 3031, 3413, 3996 | ✅️ | ✅️ |
231+
| Oblique Stereographic | custom WGS84 | ✅️ | pyproj fallback |
232+
| Oblique Mercator (Hotine) | 3375 (RSO) | implemented, disabled | pyproj fallback |
231233

232-
**Reproject performance** (end-to-end, bilinear, vs rioxarray):
234+
**Vertical datum support:** `geoid_height`, `ellipsoidal_to_orthometric`, `orthometric_to_ellipsoidal` convert between ellipsoidal (GPS) and orthometric (map/MSL) heights using EGM96 (vendored, 2.6MB) or EGM2008 (77MB, downloaded on first use). Reproject can apply vertical shifts during reprojection via the `vertical_crs` parameter.
233235

234-
| Transform | 1024x1024 | | 4096x4096 | |
235-
|:---|---:|---:|---:|---:|
236-
| | xrspatial | rioxarray | xrspatial | rioxarray |
237-
| WGS84 -> UTM 33N | 33ms | 72ms (2.2x) | 627ms | 1.09s (1.7x) |
238-
| WGS84 -> Web Mercator | 16ms | 44ms (2.9x) | 526ms | 741ms (1.4x) |
239-
| WGS84 -> Albers CONUS | 72ms | 196ms (2.7x) | 649ms | 1.78s (2.7x) |
240-
| WGS84 -> LAEA Europe | 47ms | 74ms (1.6x) | 677ms | 1.03s (1.5x) |
241-
| WGS84 -> Polar Stere S | 34ms | 580ms (17x) | 839ms | 9.13s (11x) |
236+
**Datum shift support:** Reprojection from non-WGS84 datums (NAD27, OSGB36, DHDN, MGI, ED50, BD72, CH1903, D73, AGD66, Tokyo) applies grid-based shifts from PROJ CDN (sub-metre accuracy) with 7-parameter Helmert fallback (1-5m accuracy). 14 grids are registered covering North America, UK, Germany, Austria, Spain, Netherlands, Belgium, Switzerland, Portugal, and Australia.
242237

243-
Times include coordinate transform + bilinear resampling. Speedup in parentheses is rioxarray/xrspatial. The Polar Stereographic advantage comes from rioxarray computing a much larger output grid for the same input extent.
238+
**ITRF frame support:** `itrf_transform` converts between ITRF2000, ITRF2008, ITRF2014, and ITRF2020 using 14-parameter time-dependent Helmert transforms from PROJ data files. Shifts are mm-level.
239+
240+
**Reproject performance** (reproject-only, 1024x1024, bilinear, vs rioxarray):
241+
242+
| Transform | xrspatial | rioxarray |
243+
|:---|---:|---:|
244+
| WGS84 -> Web Mercator | 23ms | 14ms |
245+
| WGS84 -> UTM 33N | 24ms | 18ms |
246+
| WGS84 -> Albers CONUS | 41ms | 33ms |
247+
| WGS84 -> LAEA Europe | 57ms | 17ms |
248+
| WGS84 -> Polar Stere S | 44ms | 38ms |
249+
| WGS84 -> LCC France | 44ms | 25ms |
250+
| WGS84 -> Ellipsoidal Merc | 27ms | 14ms |
251+
| WGS84 -> CEA EASE-Grid | 24ms | 15ms |
252+
253+
**Full pipeline** (read 3600x3600 Copernicus DEM + reproject to EPSG:3857 + write GeoTIFF):
254+
255+
| Backend | Time |
256+
|:---|---:|
257+
| NumPy | 2.7s |
258+
| CuPy GPU | 348ms |
259+
| Dask+CuPy GPU | 343ms |
260+
| rioxarray (GDAL) | 418ms |
244261

245262
**Merge performance** (4 overlapping same-CRS tiles, vs rioxarray):
246263

247264
| Tile size | xrspatial | rioxarray | Speedup |
248265
|:---|---:|---:|---:|
249-
| 512x512 | 11ms | 50ms | **4.5x** |
250-
| 1024x1024 | 82ms | 125ms | **1.5x** |
251-
| 2048x2048 | 347ms | 604ms | **1.7x** |
266+
| 512x512 | 16ms | 29ms | **1.8x** |
267+
| 1024x1024 | 52ms | 76ms | **1.5x** |
268+
| 2048x2048 | 361ms | 280ms | 0.8x |
252269

253270
Same-CRS tiles skip reprojection entirely and are placed by direct coordinate alignment.
254271

‎benchmarks/reproject_benchmark.md‎

Lines changed: 254 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,254 @@
1+
# Reproject Module: Comprehensive Benchmarks
2+
3+
Generated: 2026-03-22
4+
5+
Hardware: AMD Ryzen / NVIDIA A6000 GPU, PCIe Gen4, NVMe SSD
6+
7+
Python 3.14, NumPy, Numba, CuPy, Dask, pyproj, rioxarray (GDAL)
8+
9+
---
10+
11+
## 1. Full Pipeline Benchmark (read -> reproject -> write)
12+
13+
Source file: Copernicus DEM COG (`Copernicus_DSM_COG_10_N40_00_W075_00_DEM.tif`), 3600x3600, WGS84, deflate+floating-point predictor. Reprojected to Web Mercator (EPSG:3857). Median of 3 runs after warmup.
14+
15+
```python
16+
from xrspatial.geotiff import read_geotiff, write_geotiff
17+
from xrspatial.reproject import reproject
18+
19+
dem = read_geotiff('Copernicus_DSM_COG_10_N40_00_W075_00_DEM.tif')
20+
dem_merc = reproject(dem, 'EPSG:3857')
21+
write_geotiff(dem_merc, 'output.tif')
22+
```
23+
24+
| Backend | End-to-end time | Notes |
25+
|:--------|----------------:|:------|
26+
| NumPy | 2,723 ms | Single-threaded Numba JIT resampling |
27+
| CuPy GPU | 348 ms | CUDA kernel for coordinate transform + resampling |
28+
| Dask+CuPy GPU | 343 ms | Chunked (512) GPU pipeline |
29+
| Dask (CPU) | 10,967 ms | Chunked (512) with Dask scheduler overhead |
30+
| rioxarray (GDAL) | 418 ms | C-level warp, highly optimized |
31+
32+
The GPU path (CuPy or Dask+CuPy) is the fastest option for large rasters, running slightly faster than GDAL. The NumPy path is slower due to Python/Numba overhead in the resampling loop. The Dask CPU path has significant scheduler overhead for this single-file workload.
33+
34+
---
35+
36+
## 2. Projection Coverage and Accuracy
37+
38+
Each projection was tested with 5 geographically appropriate points. "Max error vs pyproj" measures the maximum positional difference between the Numba JIT inverse transform and pyproj's reference implementation. Errors are measured as approximate ground distance.
39+
40+
| Projection | EPSG examples | Max error vs pyproj | CPU Numba | CUDA GPU |
41+
|:-----------|:-------------|--------------------:|:---------:|:--------:|
42+
| Web Mercator | 3857 | < 0.001 mm | yes | yes |
43+
| UTM / Transverse Mercator | 326xx, 327xx | < 0.001 mm | yes | yes |
44+
| Ellipsoidal Mercator | 3395 | < 0.001 mm | yes | yes |
45+
| Lambert Conformal Conic | 2154, 2229 | 0.003 mm | yes | yes |
46+
| Albers Equal Area | 5070 | 3.5 m | yes | yes |
47+
| Cylindrical Equal Area | 6933 | 4.8 m | yes | yes |
48+
| Sinusoidal | MODIS | 0.001 mm | yes | yes |
49+
| Lambert Azimuthal Equal Area | 3035 | see note | yes | yes |
50+
| Polar Stereographic (Antarctic) | 3031 | < 0.001 mm | yes | yes |
51+
| Polar Stereographic (Arctic) | 3413 | < 0.001 mm | yes | yes |
52+
| Oblique Stereographic | custom WGS84 | < 0.001 mm | yes | fallback |
53+
| Oblique Mercator (Hotine) | 3375 | N/A | disabled | fallback |
54+
| State Plane (tmerc) | 26983 | 43 cm | yes | yes |
55+
| State Plane (LCC, ftUS) | 2229 | 19 cm | yes | yes |
56+
57+
**Notes:**
58+
- LAEA Europe (3035): The current implementation has a known latitude bias (~700m near Paris, larger at the projection's edges). This is an area for future improvement; for high-accuracy LAEA work, the pyproj fallback is used for unsupported ellipsoids.
59+
- Albers and CEA: Errors of 3-5m stem from the authalic latitude series approximation. Acceptable for most raster reprojection at typical DEM resolutions (30m+).
60+
- State Plane: Sub-metre accuracy in both tmerc and LCC variants. Unit conversion (US survey feet) is handled internally.
61+
- Oblique Stereographic: The Numba kernel exists and works for WGS84-based CRS definitions. EPSG:28992 (RD New) uses the Bessel ellipsoid without a registered datum, so it falls back to pyproj.
62+
- Oblique Mercator: Kernel implemented but disabled pending alignment with PROJ's omerc.cpp variant handling. Falls back to pyproj.
63+
64+
### Reproject-only timing (1024x1024, bilinear)
65+
66+
| Transform | xrspatial | rioxarray |
67+
|:-----------|----------:|----------:|
68+
| WGS84 -> Web Mercator | 23 ms | 14 ms |
69+
| WGS84 -> UTM 33N | 24 ms | 18 ms |
70+
| WGS84 -> Albers CONUS | 41 ms | 33 ms |
71+
| WGS84 -> LAEA Europe | 57 ms | 17 ms |
72+
| WGS84 -> Polar Stere S | 44 ms | 38 ms |
73+
| WGS84 -> LCC France | 44 ms | 25 ms |
74+
| WGS84 -> Ellipsoidal Merc | 27 ms | 14 ms |
75+
| WGS84 -> CEA EASE-Grid | 24 ms | 15 ms |
76+
77+
At 1024x1024, rioxarray (GDAL) is generally faster than the NumPy backend for reproject-only workloads. The GPU backend closes this gap and pulls ahead for larger rasters (see Section 1). The xrspatial advantage is its pure-Python stack with no GDAL dependency, four-backend dispatch (NumPy/CuPy/Dask/Dask+CuPy), and integrated vertical/datum handling.
78+
79+
### Merge timing (4 overlapping same-CRS tiles)
80+
81+
| Tile size | xrspatial | rioxarray | Speedup |
82+
|:----------|----------:|----------:|--------:|
83+
| 512x512 | 16 ms | 29 ms | 1.8x |
84+
| 1024x1024 | 52 ms | 76 ms | 1.5x |
85+
| 2048x2048 | 361 ms | 280 ms | 0.8x |
86+
87+
Same-CRS merge skips reprojection and places tiles by coordinate alignment. xrspatial is faster at small to medium sizes; rioxarray catches up at larger sizes due to its C-level copy routines.
88+
89+
---
90+
91+
## 3. Datum Shift Coverage
92+
93+
The reproject module handles horizontal datum shifts for non-WGS84 source CRS. It first tries grid-based shifts (downloaded from the PROJ CDN on first use), falling back to 7-parameter Helmert transforms when no grid is available.
94+
95+
### Grid-based shifts (sub-metre accuracy)
96+
97+
| Registry key | Grid file | Coverage | Description |
98+
|:-------------|:----------|:---------|:------------|
99+
| NAD27_CONUS | us_noaa_conus.tif | CONUS | NAD27 -> NAD83 (NADCON) |
100+
| NAD27_NADCON5_CONUS | us_noaa_nadcon5_nad27_nad83_1986_conus.tif | CONUS | NAD27 -> NAD83 (NADCON5, preferred) |
101+
| NAD27_ALASKA | us_noaa_alaska.tif | Alaska | NAD27 -> NAD83 (NADCON) |
102+
| NAD27_HAWAII | us_noaa_hawaii.tif | Hawaii | Old Hawaiian -> NAD83 |
103+
| NAD27_PRVI | us_noaa_prvi.tif | PR/USVI | NAD27 -> NAD83 |
104+
| OSGB36_UK | uk_os_OSTN15_NTv2_OSGBtoETRS.tif | UK | OSGB36 -> ETRS89 (OSTN15) |
105+
| AGD66_GDA94 | au_icsm_A66_National_13_09_01.tif | Australia NT | AGD66 -> GDA94 |
106+
| DHDN_ETRS89_DE | de_adv_BETA2007.tif | Germany | DHDN -> ETRS89 |
107+
| MGI_ETRS89_AT | at_bev_AT_GIS_GRID.tif | Austria | MGI -> ETRS89 |
108+
| ED50_ETRS89_ES | es_ign_SPED2ETV2.tif | Spain (E coast) | ED50 -> ETRS89 |
109+
| RD_ETRS89_NL | nl_nsgi_rdcorr2018.tif | Netherlands | RD -> ETRS89 |
110+
| BD72_ETRS89_BE | be_ign_bd72lb72_etrs89lb08.tif | Belgium | BD72 -> ETRS89 |
111+
| CH1903_ETRS89_CH | ch_swisstopo_CHENyx06_ETRS.tif | Switzerland | CH1903 -> ETRS89 |
112+
| D73_ETRS89_PT | pt_dgt_D73_ETRS89_geo.tif | Portugal | D73 -> ETRS89 |
113+
114+
Grids are downloaded from `cdn.proj.org` on first use and cached in `~/.cache/xrspatial/proj_grids/`. Bilinear interpolation within the grid is done via Numba JIT.
115+
116+
### Helmert fallback (1-5m accuracy)
117+
118+
When no grid covers the area, a 7-parameter (or 3-parameter) geocentric Helmert transform is applied:
119+
120+
| Datum / Ellipsoid | Type | Parameters (dx, dy, dz, rx, ry, rz, ds) |
121+
|:------------------|:-----|:-----------------------------------------|
122+
| NAD27 / Clarke 1866 | 3-param | (-8, 160, 176, 0, 0, 0, 0) |
123+
| OSGB36 / Airy | 7-param | (446.4, -125.2, 542.1, 0.15, 0.25, 0.84, -20.5) |
124+
| DHDN / Bessel | 7-param | (598.1, 73.7, 418.2, 0.20, 0.05, -2.46, 6.7) |
125+
| MGI / Bessel | 7-param | (577.3, 90.1, 463.9, 5.14, 1.47, 5.30, 2.42) |
126+
| ED50 / Intl 1924 | 7-param | (-87, -98, -121, 0, 0, 0.81, -0.38) |
127+
| BD72 / Intl 1924 | 7-param | (-106.9, 52.3, -103.7, 0.34, -0.46, 1.84, -1.27) |
128+
| CH1903 / Bessel | 3-param | (674.4, 15.1, 405.3, 0, 0, 0, 0) |
129+
| D73 / Intl 1924 | 3-param | (-239.7, 88.2, 30.5, 0, 0, 0, 0) |
130+
| AGD66 / ANS | 3-param | (-133, -48, 148, 0, 0, 0, 0) |
131+
| Tokyo / Bessel | 3-param | (-146.4, 507.3, 680.5, 0, 0, 0, 0) |
132+
133+
Grid-based accuracy is typically 0.01-0.1m; Helmert fallback accuracy is 1-5m depending on the datum.
134+
135+
---
136+
137+
## 4. Vertical Datum Support
138+
139+
The module provides geoid undulation lookup from EGM96 (vendored, 15-arcminute global grid, 2.6MB) and optionally EGM2008 (25-arcminute, 77MB, downloaded on first use).
140+
141+
### API
142+
143+
```python
144+
from xrspatial.reproject import geoid_height, ellipsoidal_to_orthometric
145+
146+
# Single point
147+
N = geoid_height(-74.0, 40.7) # New York: -32.86m
148+
149+
# Convert GPS height to map height
150+
H = ellipsoidal_to_orthometric(100.0, -74.0, 40.7) # 132.86m
151+
152+
# Batch (array)
153+
N = geoid_height(lon_array, lat_array)
154+
155+
# Raster grid
156+
from xrspatial.reproject import geoid_height_raster
157+
N_grid = geoid_height_raster(dem)
158+
```
159+
160+
### Accuracy vs pyproj geoid
161+
162+
| Location | xrspatial EGM96 (m) | pyproj EGM96 (m) | Difference |
163+
|:---------|---------------------:|------------------:|-----------:|
164+
| New York (-74.0, 40.7) | -32.86 | -32.77 | 0.09 m |
165+
| Paris (2.35, 48.85) | 44.59 | 44.57 | 0.02 m |
166+
| Tokyo (139.7, 35.7) | 35.75 | 36.80 | 1.06 m |
167+
| Null Island (0.0, 0.0) | 17.15 | 17.16 | 0.02 m |
168+
| Rio (-43.2, -22.9) | -5.59 | -5.43 | 0.16 m |
169+
170+
The 1.06m Tokyo difference is due to the 15-arcminute grid resolution in EGM96; the steep geoid gradient near Japan amplifies interpolation differences. Roundtrip accuracy (`ellipsoidal_to_orthometric` then `orthometric_to_ellipsoidal`) is exact (0.0 error).
171+
172+
### Integration with reproject
173+
174+
The `reproject` function accepts a `vertical_crs` parameter to apply vertical datum shifts during reprojection:
175+
176+
```python
177+
from xrspatial.reproject import reproject
178+
179+
# Reproject and convert ellipsoidal heights to orthometric (MSL)
180+
dem_merc = reproject(
181+
dem, 'EPSG:3857',
182+
src_vertical_crs='ellipsoidal',
183+
tgt_vertical_crs='EGM96',
184+
)
185+
```
186+
187+
---
188+
189+
## 5. ITRF Frame Support
190+
191+
Time-dependent transformations between International Terrestrial Reference Frames using 14-parameter Helmert transforms (7 static + 7 rates) from PROJ data files.
192+
193+
### Available frames
194+
195+
- ITRF2000
196+
- ITRF2008
197+
- ITRF2014
198+
- ITRF2020
199+
200+
### Example
201+
202+
```python
203+
from xrspatial.reproject import itrf_transform, itrf_frames
204+
205+
print(itrf_frames()) # ['ITRF2000', 'ITRF2008', 'ITRF2014', 'ITRF2020']
206+
207+
lon2, lat2, h2 = itrf_transform(
208+
-74.0, 40.7, 10.0,
209+
src='ITRF2014', tgt='ITRF2020', epoch=2024.0,
210+
)
211+
# -> (-73.9999999782, 40.6999999860, 9.996897)
212+
# Horizontal shift: 2.4 mm, vertical shift: -3.1 mm
213+
```
214+
215+
### All frame-pair shifts (at epoch 2020.0, location 0E 45N)
216+
217+
| Source | Target | Horizontal shift | Vertical shift |
218+
|:-------|:-------|:----------------:|:--------------:|
219+
| ITRF2000 | ITRF2008 | 33.0 mm | 32.8 mm |
220+
| ITRF2000 | ITRF2014 | 33.2 mm | 30.7 mm |
221+
| ITRF2000 | ITRF2020 | 30.5 mm | 30.0 mm |
222+
| ITRF2008 | ITRF2014 | 1.9 mm | -2.1 mm |
223+
| ITRF2008 | ITRF2020 | 2.6 mm | -2.8 mm |
224+
| ITRF2014 | ITRF2020 | 3.0 mm | -0.7 mm |
225+
226+
Shifts between recent frames (ITRF2014/2020) are at the mm level. Older frames (ITRF2000) show larger shifts (~30mm) due to accumulated tectonic motion.
227+
228+
---
229+
230+
## 6. pyproj Usage
231+
232+
The reproject module uses pyproj for metadata operations only. The heavy per-pixel work is done in Numba JIT or CUDA.
233+
234+
### What pyproj does (runs once per reproject call)
235+
236+
| Task | Cost | Description |
237+
|:-----|:-----|:------------|
238+
| CRS metadata parsing | ~1 ms | `CRS.from_user_input()`, `CRS.to_dict()`, extract projection parameters |
239+
| EPSG code lookup | ~0.1 ms | `CRS.to_epsg()` to check for known fast paths |
240+
| Output grid estimation | ~1 ms | `Transformer.transform()` on ~500 boundary points to determine output extent |
241+
| Fallback transform | per-pixel | Only used for CRS pairs without a built-in Numba/CUDA kernel |
242+
243+
### What Numba/CUDA does (the per-pixel bottleneck)
244+
245+
| Task | Implementation | Notes |
246+
|:-----|:---------------|:------|
247+
| Coordinate transforms | Numba `@njit(parallel=True)` / CUDA `@cuda.jit` | Per-pixel forward/inverse projection |
248+
| Bilinear resampling | Numba `@njit` / CUDA `@cuda.jit` | Source pixel interpolation |
249+
| Nearest-neighbor resampling | Numba `@njit` / CUDA `@cuda.jit` | Source pixel lookup |
250+
| Cubic resampling | `scipy.ndimage.map_coordinates` | CPU only (no Numba/CUDA kernel yet) |
251+
| Datum grid interpolation | Numba `@njit(parallel=True)` | Bilinear interp of NTv2/NADCON grids |
252+
| Geoid undulation interpolation | Numba `@njit(parallel=True)` | Bilinear interp of EGM96/EGM2008 grid |
253+
| 7-param Helmert datum shift | Numba `@njit(parallel=True)` | Geocentric ECEF transform |
254+
| 14-param ITRF transform | Numba `@njit(parallel=True)` | Time-dependent Helmert in ECEF |

‎xrspatial/geotiff/__init__.py‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -380,6 +380,9 @@ def write_geotiff(data: xr.DataArray | np.ndarray, path: str, *,
380380
geo_transform = _coords_to_transform(data)
381381
if epsg is None and crs is None:
382382
epsg = data.attrs.get('crs')
383+
if isinstance(epsg, str):
384+
# attrs['crs'] may be a WKT/PROJ string (e.g. from reproject)
385+
epsg = _wkt_to_epsg(epsg)
383386
if epsg is None:
384387
# Try resolving EPSG from a WKT string in attrs
385388
wkt = data.attrs.get('crs_wkt')
@@ -798,6 +801,8 @@ def write_geotiff_gpu(data, path: str, *,
798801
geo_transform = _coords_to_transform(data)
799802
if epsg is None:
800803
epsg = data.attrs.get('crs')
804+
if isinstance(epsg, str):
805+
epsg = _wkt_to_epsg(epsg)
801806
if nodata is None:
802807
nodata = data.attrs.get('nodata')
803808
if data.attrs.get('raster_type') == 'point':

0 commit comments

Comments
 (0)