Skip to content

Commit

Permalink
tests for new logic
Browse files Browse the repository at this point in the history
  • Loading branch information
@abradley60
abradley60 committed Jan 16, 2025
1 parent 79fdbb3 commit 52b1099
Showing 1 changed file with 152 additions and 57 deletions.
209 changes: 152 additions & 57 deletions tests/sar_antarctica/test_dem.py
View file
Original file line number Diff line number Diff line change
@@ -1,28 +1,31 @@
import os
import math
import pytest
import shutil
from shapely.geometry import box
import numpy as np
from pathlib import Path

from sar_antarctica.nci.preparation.dem import (
check_s1_bounds_cross_antimeridian,
split_s1_bounds_at_am_crossing,
get_target_antimeridian_projection,
make_empty_cop30m_profile,
expand_bounds,
get_cop30_dem_for_bounds
get_cop30_dem_for_bounds,
find_required_dem_paths_from_index,
find_required_dem_tile_paths_by_filename
)

from sar_antarctica.nci.preparation.create_dem_vrt import (
find_tiles,
build_tileindex
)
from sar_antarctica.nci.preparation.geoid import remove_geoid

from sar_antarctica.utils.raster import (
read_vrt_in_bounds,
bounds_from_profile
merge_raster_files,
bounds_from_profile,
expand_raster_to_bounds
)

import shapely
from shapely.geometry import box
import pytest
import math
import numpy as np
import os
import shutil
from pathlib import Path

from data.cop30m_profile import TEST_COP30_PROFILE_1, TEST_COP30_PROFILE_2

CURRENT_DIR = Path(os.path.abspath(__file__)).parent.resolve()
Expand All @@ -32,13 +35,18 @@
TEST_GEOID_TIF_PATH = CURRENT_DIR / Path('data/geoid/us_nga_egm2008_1_4326__agisoft_clipped.tif')

# TEST_COP30_FOLDER_PATH = Path('/g/data/v10/eoancillarydata-2/elevation/copernicus_30m_world/')
# TEST_COP30_INDEX_PATH = Path('/g/data/yp75/projects/ancillary/dem/copdem_tindex.gpkg')
# TEST_COP30_INDEX_PATH = Path('/g/data/yp75/projects/ancillary/dem/copdem_tindex.gpkg') # disable test_make_test_gpkg
# TEST_GEOID_TIF_PATH = Path('/g/data/yp75/projects/ancillary/geoid/us_nga_egm2008_1_4326__agisoft.tif')


def test_pytest():
assert True

def test_make_test_gpkg():
"""required to pass for other tests utilising index file"""
tiles = find_tiles(TEST_COP30_FOLDER_PATH, pattern='*.tif')
build_tileindex(tiles, TEST_COP30_INDEX_PATH)
assert os.path.exists(TEST_COP30_INDEX_PATH)

@pytest.mark.parametrize("bounds, expected", [
((163.121597, -78.632782, 172.382263, -76.383263), False), # Bounds that do not cross the antimeridian
((170, -70, 180, -76), False), # Bounds that do not cross the antimeridian
Expand Down Expand Up @@ -85,6 +93,23 @@ def test_make_empty_cop30m_profile(bounds_to_profile):
continue
assert profile[key] == target_profile[key]

@pytest.mark.parametrize("bounds, trg_shape", [
((-179.9, -79.2, -179.1, -79.1), (1,362,962))
])
def test_make_empty_cop30_dem(bounds, trg_shape):
fill_value = 0
empty_dem_profile = make_empty_cop30m_profile(bounds)
dem_arr, dem_profile = expand_raster_to_bounds(
bounds,
src_profile=empty_dem_profile,
fill_value=fill_value,
buffer_pixels=1,
save_path='empty.tif'
)
dem_bounds = bounds_from_profile(dem_profile)
assert box(*bounds).within(box(*dem_bounds))
assert dem_arr.shape == trg_shape

@pytest.mark.parametrize("bounds, buffer, expanded_bounds", [
((179.1, -79.2, 179.9, -79.1), 0, (179.09161, -79.201308, 179.900922, -79.09867)),
((179.1, -79.2, 179.9, -79.1), 0.01, (179.08161, -79.21130, 179.91092, -79.08867)),
Expand All @@ -99,48 +124,88 @@ def test_expand_bounds(bounds, buffer, expanded_bounds):
assert pytest.approx(new_bounds[2],rel=1e-5) == pytest.approx(expanded_bounds[2],rel=1e-5)
assert pytest.approx(new_bounds[3],rel=1e-5) == pytest.approx(expanded_bounds[3],rel=1e-5)

# @pytest.mark.parametrize("bounds, trg_shape, buffer_pixels", [
# ((-179.9, -79.2, -179.1, -79.1), (1,362,962), 0),
# ((-179.9, -79.2, -179.1, -79.1), (1,370,970), 4),
# ((-179.6, -79.9, -179.4, -79.5), (1, 1442, 242), 0),
# ((179.1, -79.2, 179.9, -79.1), (1,362,962), 0),
# ((179.5, -79.2, 179.6, -79.01), (1,690,126), 2),
# ((179.5, -79.2, 179.6, -79.01), (1,688,124), 1),
# ((178.1, -79.2, 179.95, -79.1), (1, 362, 2222), 0),
# ((178.1, -79.2, 179.95, -79.1), (1, 366, 2226), 2),
# ])
# def test_vrt_dem_read_for_bounds(bounds, trg_shape, buffer_pixels):
# dem_arr, dem_profile = read_vrt_in_bounds(
# vrt_path=TEST_COP30_VRT_PATH, bounds=bounds, buffer_pixels=buffer_pixels)
# dem_bounds = bounds_from_profile(dem_profile)
# assert box(*bounds).within(box(*dem_bounds))
# assert dem_arr.shape == trg_shape

# @pytest.mark.parametrize("bounds, trg_shape, geoid_ref_mean, ellipsoid_ref_mean", [
# ((-179.9, -79.2, -179.1, -79.1), (1,362,962), 44.088665, -9.830775),
# ((178.1, -79.2, 179.95, -79.1), (1, 362, 2222), 38.270348, -15.338912),
# ])
# def test_remove_geoid(bounds, trg_shape, geoid_ref_mean, ellipsoid_ref_mean):
# dem_arr, dem_profile = read_vrt_in_bounds(
# vrt_path=TEST_COP30_VRT_PATH, bounds=bounds, buffer_pixels=0)
# dem_arr_ellipsoid = remove_geoid(
# dem_arr = dem_arr,
# dem_profile = dem_profile,
# geoid_path = TEST_GEOID_TIF_PATH,
# dem_area_or_point = 'Point',
# buffer_pixels = 2,
# save_path='',
# )
# assert dem_arr.shape == dem_arr_ellipsoid.shape
# assert dem_arr.shape == trg_shape
# assert np.mean(dem_arr) == geoid_ref_mean
# assert np.mean(dem_arr_ellipsoid) == ellipsoid_ref_mean
@pytest.mark.parametrize("bounds, search_buffer, expected_tiles", [
((176.1, -79.2, 176.9, -79.1), 0, []),
((179.1, -79.2, 179.9, -79.1), 0, [str(TEST_COP30_FOLDER_PATH / 'Copernicus_DSM_COG_10_S80_00_E179_00_DEM.tif')]),
((179.1, -79.2, 179.9, -79.1), 0.3, [
str(TEST_COP30_FOLDER_PATH / 'Copernicus_DSM_COG_10_S80_00_E178_00_DEM.tif'),
str(TEST_COP30_FOLDER_PATH / 'Copernicus_DSM_COG_10_S80_00_E179_00_DEM.tif')]),
((178.1, -79.2, 179.9, -79.1), 0, [
str(TEST_COP30_FOLDER_PATH / 'Copernicus_DSM_COG_10_S80_00_E178_00_DEM.tif'),
str(TEST_COP30_FOLDER_PATH / 'Copernicus_DSM_COG_10_S80_00_E179_00_DEM.tif')]),
])
def test_find_required_dem_paths(bounds, search_buffer, expected_tiles):
# limited to files in TEST_COP30_FOLDER_PATH
tiles_from_folder = find_required_dem_tile_paths_by_filename(
bounds,
cop30_folder_path=TEST_COP30_FOLDER_PATH,
search_buffer=search_buffer,
tifs_in_subfolder=False)
tiles_from_index = find_required_dem_paths_from_index(
bounds,
cop30_index_path=TEST_COP30_INDEX_PATH,
search_buffer=search_buffer)
assert tiles_from_folder == expected_tiles
assert tiles_from_index == expected_tiles
assert tiles_from_folder == tiles_from_index


@pytest.mark.parametrize("bounds, trg_shape, buffer_pixels", [
((-179.9, -79.2, -179.1, -79.1), (1,362,962), 0),
((-179.9, -79.2, -179.1, -79.1), (1,370,970), 4),
((-179.6, -79.9, -179.4, -79.5), (1, 1442, 242), 0),
((179.1, -79.2, 179.9, -79.1), (1,362,962), 0),
((179.5, -79.2, 179.6, -79.01), (1,690,126), 2),
((179.5, -79.2, 179.6, -79.01), (1,688,124), 1),
((178.1, -79.2, 179.95, -79.1), (1, 362, 2222), 0),
((178.1, -79.2, 179.95, -79.1), (1, 366, 2226), 2),
])
def test_dem_read_for_bounds(bounds, trg_shape, buffer_pixels):
os.makedirs(CURRENT_DIR / Path('TMP'), exist_ok=True)
dem_tiles = find_required_dem_paths_from_index(bounds, cop30_index_path=TEST_COP30_INDEX_PATH)
dem_arr, dem_profile = merge_raster_files(
dem_tiles,
output_path= CURRENT_DIR / Path('TMP') / Path('TMP.tif'),
bounds=bounds,
buffer_pixels=buffer_pixels
)
shutil.rmtree(CURRENT_DIR / Path('TMP'))
dem_bounds = bounds_from_profile(dem_profile)
assert box(*bounds).within(box(*dem_bounds))
assert dem_arr.shape == trg_shape

@pytest.mark.parametrize("bounds, trg_shape, geoid_ref_mean, ellipsoid_ref_mean", [
((-179.9, -79.2, -179.1, -79.1), (1,362,962), 44.088665, -9.830775),
((178.1, -79.2, 179.95, -79.1), (1, 362, 2222), 38.270348, -15.338912),
])
def test_remove_geoid(bounds, trg_shape, geoid_ref_mean, ellipsoid_ref_mean):
os.makedirs(CURRENT_DIR / Path('TMP'), exist_ok=True)
dem_tiles = find_required_dem_paths_from_index(bounds, cop30_index_path=TEST_COP30_INDEX_PATH)
dem_arr, dem_profile = merge_raster_files(
dem_tiles,
output_path= CURRENT_DIR / Path('TMP') / Path('TMP.tif'),
bounds=bounds,
buffer_pixels=0
)
dem_arr_ellipsoid = remove_geoid(
dem_arr = dem_arr,
dem_profile = dem_profile,
geoid_path = TEST_GEOID_TIF_PATH,
dem_area_or_point = 'Point',
buffer_pixels = 2,
save_path='',
)
shutil.rmtree(CURRENT_DIR / Path('TMP'))
assert dem_arr.shape == dem_arr_ellipsoid.shape
assert dem_arr.shape == trg_shape
assert np.mean(dem_arr) == pytest.approx(geoid_ref_mean,rel=1e-5)
assert np.mean(dem_arr_ellipsoid) == pytest.approx(ellipsoid_ref_mean,rel=1e-5)

@pytest.mark.parametrize("bounds, ellipsoid_heights, trg_shape, trg_crs", [
((-179.9, -79.2, -179.1, -79.1), False, (1, 366, 966), 4326),
((178.1, -79.2, 179.95, -79.1), False, (1, 366, 2226), 4326),
((-179.2, -79.2, 179.1, -79.1), False, (1, 1266, 1502), 3031), # across antimeridian
((-179.2, -79.2, 179.1, -79.1), True, (1, 1266, 1502), 3031) # across antimeridian
((-179.2, -79.2, 179.1, -79.1), False, (1, 1396, 1676), 3031), # across antimeridian
((-179.2, -79.2, 179.1, -79.1), True, (1, 1396, 1676), 3031) # across antimeridian
])
def test_get_cop30_dem_for_bounds(bounds, ellipsoid_heights, trg_shape, trg_crs):
os.makedirs(CURRENT_DIR / Path('TMP'), exist_ok=True)
Expand All @@ -149,11 +214,41 @@ def test_get_cop30_dem_for_bounds(bounds, ellipsoid_heights, trg_shape, trg_crs)
save_path= CURRENT_DIR / Path('TMP') / Path('TMP.tif'),
ellipsoid_heights=ellipsoid_heights,
buffer_pixels=2,
#cop30_index_path=TEST_COP30_INDEX_PATH,
cop30_index_path=TEST_COP30_INDEX_PATH,
cop30_folder_path=TEST_COP30_FOLDER_PATH,
geoid_tif_path=TEST_GEOID_TIF_PATH,
adjust_for_high_lat_and_buffer=False
)
#shutil.rmtree(CURRENT_DIR / Path('TMP'))
assert dem_arr.shape == trg_shape
assert dem_profile['crs'].to_epsg() == trg_crs
#shutil.rmtree(CURRENT_DIR / Path('TMP'))

if __name__ == "__main__":

os.makedirs(CURRENT_DIR / Path('TMP'), exist_ok=True)
bounds = (-179.9, -79.2, -179.1, -79.1)
dem_arr, dem_profile = get_cop30_dem_for_bounds(
bounds,
save_path= CURRENT_DIR / Path('TMP') / Path('TMP.tif'),
ellipsoid_heights=True,
buffer_pixels=0,
cop30_index_path=TEST_COP30_INDEX_PATH,
cop30_folder_path=TEST_COP30_FOLDER_PATH,
geoid_tif_path=TEST_GEOID_TIF_PATH,
adjust_for_high_lat_and_buffer=False
)
print(dem_profile)
empty_dem_profile = make_empty_cop30m_profile((0, -90, 1, -86))
print(empty_dem_profile)
dem_arr, dem_profile = expand_raster_to_bounds(
bounds,
src_profile=empty_dem_profile,
fill_value=0,
buffer_pixels=1,
save_path=CURRENT_DIR / Path('TMP') / Path('EMPTY.tif'),
)
dem_bounds = bounds_from_profile(dem_profile)
print(dem_profile)



0 comments on commit 52b1099

Please sign in to comment.