BurnP3 / BurnP3Plus

A SyncroSim package to explore fire risk and susceptibility across a landscape.
https://burnp3.github.io/BurnP3Plus/
MIT License
7 stars 4 forks source link

[Burn-P3+ Bug]: BP3+ mulches projections #15

Open BadgerOnABike opened 1 year ago

BadgerOnABike commented 1 year ago

Contact Details

cptcitrus or badgeronabike

What happened?

When we import rasters, fuels in this case it garbles the CRS to the point of failure.

@cptcitrus brought it up and if he has more can say so here!

What component are you seeing the problem on?

No response

Relevant log output

No response

Approvals Process

cptcitrus commented 1 year ago

The first failure occurs in ignitions.R on line tryCatch(fuelsRaster %>% project("epsg:4326"),...

The raster, when loaded directly from file, has this CRS:

Coordinate System is:
PROJCRS["NAD83 / Canada Atlas Lambert",
    BASEGEOGCRS["NAD83",
        DATUM["North American Datum 1983",
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4269]],
    CONVERSION["Canada Atlas Lambert",
        METHOD["Lambert Conic Conformal (2SP)",
            ID["EPSG",9802]],
        PARAMETER["Latitude of false origin",49,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8821]],
        PARAMETER["Longitude of false origin",-95,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8822]],
        PARAMETER["Latitude of 1st standard parallel",49,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8823]],
        PARAMETER["Latitude of 2nd standard parallel",77,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8824]],
        PARAMETER["Easting at false origin",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8826]],
        PARAMETER["Northing at false origin",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8827]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],

When the fuels raster is loaded from the Syncrosim scenario, the CRS is lost. For example: fuelsRaster <- datasheet(myScenario, "burnP3Plus_LandscapeRasters")[["FuelGridFileName"]] %>% rast()


ENGCRS["NAD83 / Canada Atlas Lambert",
    EDATUM[""],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]```
shreeramsenthi commented 1 year ago

Good catch! That certainly does seem like a bug.

Raphael posted a similar (but not identical) issue on discord and while looking into it I found that terra's help documentation for the crs function suggests that they no longer formally support "proj-string notation" and instead suggest the WKT2 or <authority>:<code> notation instead. (See here, under the Notes section: https://rdrr.io/cran/terra/man/crs.html). Would it be reasonable to overwrite the CRS of your input layers (you only need to for the fuels layer) to one of the other two formats? I don't know if it ultimately worked for Raphael, but I just asked him on the discord server and I'll report back.

edit: oh, haha sorry Quinn, didn't recognize your github handle. Guess you already know about the formatting thing. Did that not work for you?