Skip to content

Commit

Permalink
fix filtering of nodata to avoid index misalignment
Browse files Browse the repository at this point in the history
  • Loading branch information
jdries committed Dec 5, 2024
1 parent 611f4fc commit 059e2c6
Show file tree
Hide file tree
Showing 2 changed files with 28 additions and 12 deletions.
6 changes: 4 additions & 2 deletions openeogeotrellis/collections/sentinel3.py
Original file line number Diff line number Diff line change
Expand Up @@ -452,8 +452,10 @@ def variable_name(band_name: str): # actually, the file without the .nc extensi
def readAndReproject(data_mask_,LUT):
band_data, band_settings = read_band(in_file, in_band=variable_name(band_name), data_mask=data_mask_)
flat = band_data.flatten()
flat_mask = data_mask_.where(data_mask_, other=False, drop=True).data.flatten()
filtered = flat[flat_mask]
#inconvenient approach for compatibility with old xarray version
flat_mask = data_mask_.where(data_mask_, drop=True).data.flatten()
flat_mask[np.isnan(flat_mask)] = 0
filtered = flat[flat_mask.astype(np.bool_)]
return band_settings,apply_LUT_on_band(filtered, LUT, band_settings.get('_FillValue', None))

if product_type == SLSTR_PRODUCT_TYPE and band_name in ["solar_azimuth_tn", "solar_zenith_tn", "sat_azimuth_tn", "sat_zenith_tn",]:
Expand Down
34 changes: 24 additions & 10 deletions tests/data_collections/test_sentinel3.py
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,8 @@

from openeogeotrellis.collections.sentinel3 import *
from numpy.testing import assert_allclose
import pytest

from tests.data import get_test_data_file

def test_read_single():

Expand Down Expand Up @@ -93,8 +94,7 @@ def test_read_single_edge_with_some_data():
"key":{
"col":10,
"row":10,
"instant":100,

"instant":1717326516089,
},
"key_extent":{
"xmin":13.86,"xmax":18.10,"ymin":47.096,"ymax":47.597
Expand All @@ -107,17 +107,31 @@ def test_read_single_edge_with_some_data():
result = read_product((product_dir, tiles), SLSTR_PRODUCT_TYPE, ["LST_in:LST", "geometry_tn:solar_zenith_tn"], 1024, True, 0.008928571428571)

assert len(result) == 1

spacetimekey = result[0][0]

#instant should be rounded to the minute
assert spacetimekey.instant == datetime(2024, 6, 2, 11, 8,0)
from rasterio.transform import from_origin, from_bounds

arr = result[0][1].cells

transform = from_bounds(13.86, 47.096, 18.10, 47.597, arr.shape[2], arr.shape[1])

new_dataset = rasterio.open('ref_file_edge3.tif', 'w', driver='GTiff',
height=arr.shape[2], width=arr.shape[2],
count=2, dtype=str(arr.dtype),
crs=CRS.from_epsg(4326),
transform=transform)
# new_dataset = rasterio.open('ref_file_edge3.tif', 'w', driver='GTiff',
# height=arr.shape[1], width=arr.shape[2],
# count=2, dtype=str(arr.dtype),
# crs=CRS.from_epsg(4326),
# transform=transform)
#
# new_dataset.write(arr)
# new_dataset.close()

expected_path = get_test_data_file("binary/Sentinel-3/sentinel3_edge_ref.tif")

new_dataset.write(arr)
new_dataset.close()
with rasterio.open(expected_path) as ds_ref:
expected_result = ds_ref.read()
#with rasterio.open(filename) as ds:
actual_result = arr
assert expected_result.shape == actual_result.shape
assert_allclose(expected_result, actual_result, atol=0.00001)

0 comments on commit 059e2c6

Please sign in to comment.