Reconcile coordinate systems before extracting, and say what got nothing - #20
Merged
Conversation
attach_bathymetry() passed st_coordinates(dat) straight to terra::extract(). Those coordinates are in whatever system the object carries, and extract() reads them as the raster's own. Observations in a projected CRS were therefore read as though their metres were degrees: every point landed outside a lon/lat grid and came back NA. Demonstrated on three points, identical but for their CRS: EPSG:4326 DEPTH: 100 150 200 EPSG:32619 DEPTH: NA NA NA coords sent: 500000, 4649776 That is observations with no depth in an area that plainly has depth, with nothing said about it. sf input is now reprojected to the raster's CRS first, and no CRS at all warns rather than being assumed right. A plain data frame carries no CRS to check, so the only available signal is magnitude - and it is unambiguous, since a longitude cannot exceed 180. Those are refused with the conversion to make. NA depths are now also reported with their cause, because the two causes have different remedies and neither is guessable from a column of NAs. Outside the grid means the bounding box was drawn too small. Inside it means ETOPO calls that cell land, which at 4 arc-minutes happens readily to inshore stations, since a cell roughly 7 km across is land if most of it is. That second one is easy to mistake for a bug and is not. matchData() already reprojected correctly, but failed on a missing CRS with sf's "crs not found: is it missing?" - true, and silent about which object and what to do. It now names both. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Diagnosing "areas with observations but no depth".
The bug
attach_bathymetry()passedst_coordinates(dat)straight toterra::extract(). Those coordinates are in whatever system the object carries, andextract()reads them as the raster's own. Observations in a projected CRS were read as though their metres were degrees — every point landed outside a lon/lat grid and came backNA.Three points, identical but for their CRS:
That is exactly "observations with no depth in an area that plainly has depth", with nothing said about it.
sfinput is now reprojected to the raster's CRS first. No CRS at all warns rather than being assumed correct.A plain data frame carries no CRS to check, so magnitude is the only available signal — an unambiguous one, since a longitude cannot exceed 180. Those are refused, with the conversion to make.
NA depths now say why
The two causes have different remedies, and neither is guessable from a column of
NAs:That second cause may well be what you are seeing if the missing points are coastal — it is not a bug, and no CRS fix will move it. A finer
resolutionwill.matchData
Already reprojected correctly — verified projected and geographic observations match the same cells. But a missing CRS failed with sf's
"crs not found: is it missing?", which is true and says neither which object nor what to do. It now names both sides.Verification
R CMD check— Status: OK. Full suite green. Installed into both R 4.6.1 and 4.3.2.🤖 Generated with Claude Code