Skip to content

Commit dff6dd3

Browse files
zxdawnzxdawn
committed
Add DEM to RPC transform in the EnMAP coords correction.
- Fix #148
1 parent 53aadd4 commit dff6dd3

File tree

3 files changed

+198
-47
lines changed

3 files changed

+198
-47
lines changed

‎hypergas/dem.py‎

Lines changed: 134 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,134 @@
1+
#!/usr/bin/env python
2+
# -*- coding: utf-8 -*-
3+
# Copyright (c) 2023-2024 HyperGas developers
4+
#
5+
# This file is part of hypergas.
6+
#
7+
# hypergas is a library to retrieve trace gases from hyperspectral satellite data
8+
"""Download DEM data for the region observed by the hyperspectral satellite."""
9+
10+
import os
11+
import logging
12+
import rasterio
13+
from pathlib import Path
14+
from dem_stitcher.stitcher import stitch_dem
15+
16+
17+
LOG = logging.getLogger(__name__)
18+
19+
20+
class DEM():
21+
"""
22+
Download and store a Digital Elevation Model (DEM) for a given bounding box.
23+
"""
24+
25+
def __init__(
26+
self,
27+
filename: str,
28+
bounds: list,
29+
buffer: float = 0.01,
30+
dem_name: str = "glo_30",
31+
dst_area_or_point: str = "Point",
32+
dst_ellipsoidal_height: bool = False,
33+
):
34+
"""
35+
Parameters
36+
----------
37+
filename : str
38+
Path to the input file. The DEM will be saved as:
39+
<filename_root>_dem.tif
40+
bounds : list[float]
41+
Bounding box in EPSG:4326 format:
42+
[xmin, ymin, xmax, ymax] (longitude, latitude).
43+
buffer : float, optional
44+
Buffer (in degrees) added around the bounding box.
45+
Default is 0.01.
46+
dem_name : str, optional
47+
Name of the DEM dataset supported by `dem_stitcher`.
48+
Default is "glo_30".
49+
dst_area_or_point : str, optional
50+
Pixel referencing type:
51+
- "Area": pixel value represents area (upper-left corner reference, GDAL default)
52+
- "Point": pixel value represents the pixel center
53+
Default is "Point".
54+
dst_ellipsoidal_height : bool, optional
55+
If True, output heights are ellipsoidal.
56+
If False, heights are relative to the geoid.
57+
Default is False.
58+
"""
59+
60+
# ---- Validate bounds ----
61+
if len(bounds) != 4:
62+
raise ValueError("`bounds` must be [xmin, ymin, xmax, ymax].")
63+
64+
xmin, ymin, xmax, ymax = bounds
65+
if xmin >= xmax or ymin >= ymax:
66+
raise ValueError("Invalid bounds: xmin < xmax and ymin < ymax required.")
67+
68+
# ---- Validate dst_area_or_point ----
69+
if dst_area_or_point not in {"Area", "Point"}:
70+
raise ValueError("`dst_area_or_point` must be 'Area' or 'Point'.")
71+
72+
# ---- Prepare output path ----
73+
dem_path = Path(filename)
74+
self.file_dem = dem_path.with_name(f"{dem_path.stem}_dem.tif")
75+
76+
self.dem_name = dem_name
77+
self.dst_area_or_point = dst_area_or_point
78+
self.dst_ellipsoidal_height = dst_ellipsoidal_height
79+
80+
# Apply buffer to bounds
81+
self.bounds = [
82+
xmin - buffer,
83+
ymin - buffer,
84+
xmax + buffer,
85+
ymax + buffer,
86+
]
87+
88+
def download(self, overwrite: bool = False) -> Path:
89+
"""
90+
Download the DEM and save it to disk.
91+
92+
Parameters
93+
----------
94+
overwrite : bool, optional
95+
If True, force re-download even if file exists.
96+
Default is False.
97+
98+
Returns
99+
-------
100+
Path
101+
Path to the downloaded DEM file.
102+
"""
103+
104+
if self.file_dem.exists() and not overwrite:
105+
LOG.debug(f"Skipping download. DEM already exists: {self.file_dem}")
106+
return self.file_dem
107+
108+
LOG.debug(
109+
f"Downloading DEM '{self.dem_name}' "
110+
f"for bounds {self.bounds}"
111+
)
112+
113+
# Fetch DEM data
114+
dem_array, profile = stitch_dem(
115+
self.bounds,
116+
dem_name=self.dem_name,
117+
dst_ellipsoidal_height=self.dst_ellipsoidal_height,
118+
dst_area_or_point=self.dst_area_or_point,
119+
)
120+
121+
# Ensure output directory exists
122+
self.file_dem.parent.mkdir(parents=True, exist_ok=True)
123+
124+
# Write GeoTIFF
125+
with rasterio.open(self.file_dem, "w", **profile) as dst:
126+
dst.write(dem_array, 1)
127+
dst.update_tags(
128+
AREA_OR_POINT=self.dst_area_or_point,
129+
source=self.dem_name,
130+
)
131+
132+
LOG.info(f"DEM saved to {self.file_dem}")
133+
134+
return self.file_dem

‎hypergas/hyper.py‎

Lines changed: 58 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -18,6 +18,7 @@
1818
from rasterio.transform import RPCTransformer
1919

2020
from .angles import Angle2D, compute_raa, compute_sga
21+
from .dem import DEM
2122
from .denoise import Denoise
2223
from .hsi2rgb import Hsi2rgb
2324
from .orthorectification import Ortho
@@ -260,28 +261,50 @@ def _scale_units(self, units):
260261

261262
return scale
262263

263-
def _rpc_transform(self, scn, rpc):
264+
def _rpc_transform(self, scn, rpc, bounds):
265+
"""
266+
Compute a more accurate longitude/latitude transform for a scene
267+
using Rational Polynomial Coefficients (RPC).
268+
269+
Parameters
270+
----------
271+
scn : Scene
272+
Source scene object to be geolocated.
273+
274+
rpc : RPC
275+
Rational Polynomial Coefficients describing the sensor model.
276+
277+
bounds : list[float]
278+
Bounding box in EPSG:4326 coordinates formatted as
279+
[xmin, ymin, xmax, ymax], where x = longitude and y = latitude.
280+
281+
Returns
282+
-------
283+
area : pyresample.geometry.AreaDefinition
284+
Pyresample area definition describing the transformed scene geometry.
285+
"""
286+
# Download DEM data
287+
dem = DEM(self.filename[0], bounds)
288+
dem.download()
289+
264290
# Create RPCTransformer
265-
transformer = RPCTransformer(rpc)
291+
transformer = RPCTransformer(rpc, rpc_dem=dem.file_dem)
266292

267-
# Create 2D coordinate grids for image of shape (1024, 1000)
293+
# Get image dimensions
268294
height, width = scn['radiance'].sizes['y'], scn['radiance'].sizes['x']
269295

270-
# Create meshgrid of pixel coordinates
271-
rows, cols = np.mgrid[0:height, 0:width]
296+
# Create pixel coordinate grids
297+
x_coords = np.arange(width)
298+
y_coords = np.arange(height)
299+
xx, yy = np.meshgrid(x_coords, y_coords)
272300

273-
# Flatten to 1D arrays
274-
cols_flat = cols.ravel() # sample/column coordinates (x)
275-
rows_flat = rows.ravel() # line/row coordinates (y)
276-
heights_flat = np.zeros_like(cols_flat, dtype=float) # elevation (ground level = 0)
301+
# Transform pixel coordinates to lon/lat using xy() method
302+
# xy() takes (row, col) and returns (x, y) which is (lon, lat)
303+
lons, lats = transformer.xy(yy, xx)
277304

278-
# Transform from image coordinates to lon, lat
279-
# Try xy() method instead - it transforms pixel to geographic
280-
lons_flat, lats_flat = transformer.xy(rows_flat, cols_flat, heights_flat)
281-
282-
# Reshape back to 2D arrays
283-
lons = np.array(lons_flat).reshape(height, width)
284-
lats = np.array(lats_flat).reshape(height, width)
305+
# Reshape back to 2D
306+
lons = lons.reshape(height, width)
307+
lats = lats.reshape(height, width)
285308

286309
area = SwathDefinition(lons, lats)
287310
area.name = '_'.join([scn['radiance'].attrs['platform_name'], str(scn['radiance'].attrs['start_time'])])
@@ -329,13 +352,26 @@ def load(self, drop_waterbands=True):
329352

330353
# recalculate the lon and lat for EnMAP
331354
# because the default yc and xc coords of TIFF files are not accurate
332-
if self.reader == 'hsi_l1b':
333-
rpc_swir = scn['rpc_coef_swir'].isel(bands_swir=0).item()
334-
new_area = self._rpc_transform(scn, rpc_swir)
355+
if self.reader == "hsi_l1b":
356+
# Get default lon/lat from radiance area
357+
lons_default, lats_default = scn["radiance"].attrs["area"].get_lonlats()
358+
359+
bounds = [
360+
lons_default.min(),
361+
lats_default.min(),
362+
lons_default.max(),
363+
lats_default.max(),
364+
]
365+
366+
LOG.info("Calculating lon/lat using RPC and DEM data")
367+
368+
rpc_swir = scn["rpc_coef_swir"].isel(bands_swir=0).item()
369+
new_area = self._rpc_transform(scn, rpc_swir, bounds)
370+
371+
# Replace area definition for all datasets that contain an area attribute
335372
for key in scn.keys():
336-
if 'area' in scn[key].attrs:
337-
# Replace with the new area definition
338-
scn[key].attrs['area'] = new_area
373+
if "area" in scn[key].attrs:
374+
scn[key].attrs["area"] = new_area
339375

340376
# get attrs
341377
self.start_time = scn['radiance'].attrs['start_time']

‎hypergas/orthorectification.py‎

Lines changed: 6 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -25,6 +25,8 @@
2525
from rasterio.warp import transform
2626
from rasterio.control import GroundControlPoint as GCP
2727

28+
from .dem import DEM
29+
2830
LOG = logging.getLogger(__name__)
2931

3032

@@ -113,32 +115,11 @@ def _get_default_transform(self):
113115

114116
def _download_dem(self):
115117
"""Download GLO-30 DEM data."""
116-
self.dem_name = 'glo_30'
117-
118118
# download DEM data and update config
119-
dem_path = Path(self.data.attrs['filename'])
120-
self.file_dem = dem_path.with_name(dem_path.stem+'_dem.tif')
121-
122-
if not os.path.exists(self.file_dem):
123-
dst_area_or_point = 'Point'
124-
dst_ellipsoidal_height = False
125-
126-
# get the DEM data (Copernicus GLO-30)
127-
LOG.debug('Downloading GLO-30 DEM data using stitch_dem')
128-
X, p = stitch_dem(self.bounds,
129-
dem_name=self.dem_name,
130-
dst_ellipsoidal_height=dst_ellipsoidal_height,
131-
dst_area_or_point=dst_area_or_point)
132-
133-
# export to tif file
134-
with rasterio.open(self.file_dem, 'w', **p) as ds:
135-
ds.write(X, 1)
136-
ds.update_tags(AREA_OR_POINT=dst_area_or_point, source=self.dem_name)
137-
else:
138-
LOG.debug(f'Reading the existed {self.file_dem}')
139-
with rasterio.open(self.file_dem) as ds:
140-
self.dem_name = ds.tags().get('source')
141-
119+
dem = DEM(self.data.attrs['filename'], self.bounds)
120+
self.file_dem = dem.file_dem
121+
self.dem_name = dem.dem_name
122+
dem.download()
142123

143124
def _utm_epsg(self):
144125
"""Find the suitable UTM epsg code based on lons and lats

0 commit comments

Comments
 (0)