Repository navigation
load_degrad(): two pre-existing CRS/datum bugs found while fixing operation_crs #350
Description
While replacing the hardcoded operation_crs in R/degrad.R (the same area-bias fix already made in load_deter(): the old +proj=poly ... +ellps=aust_SA string is not equal-area and has a mistyped ellipsoid), I found two pre-existing problems in how load_degrad() handles the CRS of the raw DEGRAD files. Neither is caused by that fix: both behave identically with the current master code and with the fixed code. Nothing is changed for them for now; this issue only documents them.
Related: #283 reported the same arguments have different crs error for load_degrad() and was closed without a diagnosis in the thread. Bug 1 below is its cause.
Measured on 2026-10-12 (R 4.3.3, sf 1.0.19, GDAL 3.8.2, PROJ 9.3.1; R/degrad.R as on master at c64f0de).
1. The 10 raw shapefiles do not share a CRS, but the code assumes they all lack one
CRS of each year's shapefile as read by sf (load_degrad(..., raw_data = TRUE)):
| Years | CRS read from the file |
|---|---|
| 2007-2009 | none (NA) |
| 2010-2014 | SAD69 (EPSG:4618), embedded .prj carrying a TOWGS84 |
| 2015-2016 | SIRGAS 2000 |
load_degrad() row-binds the requested years right after downloading them, and sf refuses to bind different CRSs, so any time_period that crosses a group boundary fails. Ranges inside one group work, which is why the roxygen example (2010:2012) does:
library(datazoom.amazonia)
load_degrad("degrad", time_period = 2010:2012) # OK
load_degrad("degrad", time_period = 2010:2014) # OK
load_degrad("degrad", time_period = 2009:2010) # Error: arguments have different crs
load_degrad("degrad", time_period = 2014:2015) # Error: arguments have different crs
load_degrad("degrad", time_period = 2007:2016) # Error: arguments have different crsAll five reproduce with master's R/degrad.R and with the operation_crs fix. Loading the ten years one at a time works.
Also, the block if (is.na(sf::st_crs(dat$geometry))) (which assigns a SAD69 approximation with an explicit +towgs84) is written as if the files never carried a CRS; it can only ever apply to a load made solely of 2007-2009. actions/scripts/check_all_loads.R loads 2016 alone, so none of this is exercised in CI.
2. The SAD69 -> SIRGAS 2000 datum shift is never applied
Neither the current operation_crs (+proj=poly ... +ellps=aust_SA) nor the IBGE Albers string that replaces it (+proj=aea ... +ellps=GRS80) names a datum (PROJ reads the new one as "Unknown based on GRS 1980 ellipsoid"). sf::st_transform() to such a CRS is a zero-shift ballpark conversion, so the input datum is ignored.
Measured at three points ((-60,-5), (-52,-3), (-70,-9)), for both targets: the same longitude/latitude read as SAD69 (the 2012 and 2014 files' own CRS), as SIRGAS 2000, or through the +towgs84 string of the missing-CRS fallback block lands on exactly the same projected point (0.000 m apart). For scale, transforming the 2012 file's CRS to EPSG:4674 (SIRGAS 2000, which has a datum) moves the points by 64.6 to 74.5 m.
So the polygons of the SAD69-era years (2007-2014; for 2007-2009 SAD69 is an assumption, those files carry no CRS) are overlaid on SIRGAS 2000 municipality boundaries without that correction. Areas are unaffected (a translation does not change an area); the assignment of polygons to municipalities near boundaries is.
Measured effect on 2012, using IBGE's own published WKT for the same Albers projection (it names D_SIRGAS2000; IBGE, Malha Municipal Digital e Areas Territoriais 2024, Notas metodologicas 01/2025, Apendice B) as operation_crs instead of the proj4 string:
- total area 8,532.66 -> 8,532.66 km2 (-0.0001%); output rows 9,111 -> 9,113;
- per-municipality area, over the 290 municipalities with more than 1 km2 of DEGRAD area: median change 0.03%; largest relative change 13.2% (Santa Maria das Barreiras, PA: 1.97 -> 1.71 km2); largest absolute change 0.72 km2 (Barra do Corda, MA: 475.12 -> 474.40 km2, -0.15%).
2016 (already SIRGAS 2000) is identical under both definitions. load_deter() is not affected: its raw files are already SIRGAS 2000.
Where this is documented
In the # FRAGILE: note above if (is.na(sf::st_crs(dat$geometry))) in R/degrad.R (branch update-degrad), and in the "Known limitations" section of the project's shared SKILL.md maintenance notes.
Activity
Fix imprecise municipality stat in R/degrad.R FRAGILE tag
Fix DEGRAD cross-group CRS crash and unapplied datum shift (closes #350)
phlebailly6 commentedon Oct 17, 2026
Fixed on
update-degradin 52bcc7d (pushed, not yet merged intomaster).Both problems above are addressed by one change in
R/degrad.R: a newharmonize_degrad_year_crs()runs right after the download and beforebind_rows(), and applies onesf::st_transform()per year to SIRGAS 2000 (EPSG:4674). Years with no CRS (2007-2009) first get the pre-existing SAD69+towgs84approximation assigned, then are transformed like every other year. The downstreamoperation_crs(IBGE Albers) is unchanged.Bug 1 (crash). Re-run on 2026-10-16 against the real INPE files, the three failing cases from the report now succeed:
2009:2010(32,153 rows),2014:2015(59,033 rows) and2007:2016(231,453 rows, 145,554.45 km2). In each case rows and km2 are exactly equal to the sum of the corresponding single-year loads, and all geometries are valid.Bug 2 (datum shift). Measured polygon by polygon on the real data, the SAD69-era years (2007-2014) move 59.5-75.0 m (median 65-69 m depending on the year) compared with reading the raw coordinates as SIRGAS 2000; 2015-2016 move 0.0 m. For 2007-2009 the reference is the assigned fallback, since those files carry no CRS. The shift comes from the TOWGS84 embedded in the 2010-2014
.prjfiles (a BOUNDCRS), not from the separately registered EPSG:4618 -> 4674 operation; the two are 0.3 m apart (at (-60, -5): 69.49 m vs 69.83 m).Effect of the datum shift alone (same Albers target, before vs after this change, all 10 years):
Erratum: the message of 52bcc7d gives the datum-only area effect as -0.0009%; the sign is wrong, it is +0.0009% (new / old - 1).
Combined effect against the original polyconic code (Albers + datum): +0.41% total area, unchanged from the 2026-10-12 measurement (+0.30% to +0.74% by year), and 231,453 output rows vs 227,218 raw polygons (polygons split at municipal boundaries or dropped outside the Legal Amazon).
raw_data = TRUEis untouched: the harmonization runs after its early return, and the output is byte-identical to the previous behaviour on2010:2012,2007:2009,2015:2016and2010:2014.Side effect to be aware of (not introduced by this change). The 2007-2009 raw files have no
view_date. A multi-year load that crosses them used to crash; it now succeeds and returnsNAindate(datawithlanguage = "pt") for those years (7,628NAfor2009:2010, 92,760 for2007:2016). Single-year loads of 2007-2009 simply have nodatecolumn. The fix and this behaviour are documented inNEWS.mdand in the three# FRAGILE:notes inR/degrad.R.Why this issue stays open. Two things are still pending:
update-degradintomaster;tests/testthat/test-degrad-crs.R: three synthetic polygons at the same lon/lat declared with no CRS, SAD69 + TOWGS84 and SIRGAS 2000,external_download()mocked; it fails on the pre-fix code witharguments have different crs, passes now, runs in about 1 s with no network) and a wideractions/scripts/check_all_loads.Rcheck that loadsload_degradforc(2009, 2010, 2016)instead of 2016 alone (one call crossing both CRS group boundaries; about 75 s against the real servers; fails on the pre-fix code with the same error). They will arrive in a follow-up commit on the same branch, linked here once pushed.Add offline regression test and widen check_all_loads for load_degrad…