Skip to content

Removes GeoTIFF mask for HSIA production coverage#714

Merged
charparr merged 1 commit intomainfrom
remove_hsia_geotiff
Apr 9, 2026
Merged

Removes GeoTIFF mask for HSIA production coverage#714
charparr merged 1 commit intomainfrom
remove_hsia_geotiff

Conversation

@BobTorgerson
Copy link
Copy Markdown
Contributor

This PR removes the GeoTIFF mask for the hsia_arctic_production Rasdaman coverage as it was masking out valid pixels of communities on the shoreline in our https://seaiceatlas.org website.
Screenshot 2026-04-09 at 10 17 18 AM

You can see in the picture above that anywhere that color is being shown through the mask is a valid pixel that was incorrectly being marked as having no data. The code that was being used for this is in the check_geotiff function found here:

def check_geotiffs(lat, lon, coverages):
"""Load a binary GeoTIFF mask corresponding to the coverage(s) requested, then use this to check if lat/lon has data available.
Args:
lat (int or float): latitude
lon (int or float): longitude
coverages (list): list of coverages to check for data availability
Returns:
True if valid, or HTTP 404 status code if no data was found
"""
for coverage in coverages:
# Use the same GeoTIFF for all CMIP6 downscaled coverages.
if coverage.startswith("cmip6_downscaled_"):
coverage = "cmip6_downscaled"
reference_geotiff = "geotiffs/" + coverage + ".tif"
# Do not perform GeoTIFF check if the file does not exist.
if not os.path.isfile(reference_geotiff):
return True
# Do not perform GeoTIFF check if the file does not open properly.
# This seems safer than the alternative of hiding data due to a corrupt file.
try:
with rasterio.open(reference_geotiff) as dataset:
if coverage in geotiff_projections:
crs = geotiff_projections[coverage]
else:
crs = "EPSG:3338"
if crs == "EPSG:4326":
# check if lat/lon is within the dataset's geographic bounds
if (
dataset.bounds.left <= lon <= dataset.bounds.right
and dataset.bounds.bottom <= lat <= dataset.bounds.top
):
row, col = dataset.index(lon, lat)
if dataset.read(1)[row, col] == 1:
return True
else:
x, y = project_latlon(lat, lon, crs)
row, col = dataset.index(x, y)
if 0 <= row < dataset.height and 0 <= col < dataset.width:
if dataset.read(1)[row, col] == 1:
return True
except:
return True
return 404

NOTE: When this goes into production, we are going to want to clear the Earth Maps cache of the /seaice endpoint so we can get all valid data back again.

@BobTorgerson BobTorgerson requested a review from charparr April 9, 2026 18:51
@charparr
Copy link
Copy Markdown
Member

charparr commented Apr 9, 2026

I was able to reproduce Bob's findings by taking the mask GeoTIFF that was removed in this PR, and comparing to a GeoTIFF of the production data (December 2024), and I observed the mask layer knocking out pixels of valid data.

Copy link
Copy Markdown
Member

@charparr charparr left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Image

Easy to test with a side-by-side requests for a location where we'd expect data, but we don't get it in prod: https://seaiceatlas.org/report/58.760/-69.500#report

@charparr charparr merged commit a316ca0 into main Apr 9, 2026
1 check failed
@charparr charparr deleted the remove_hsia_geotiff branch April 9, 2026 19:12
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants