Open cmbarton opened 4 years ago
I believe the version that multiplies by 1000 is correct as we want g/m2 at this stage. I've now fixed this in the nonanthropogenic version.
We should formalize the whole charcoal conversion process here to include in the documentation without having to refer to Grant's cheat sheet.
I believe the version that multiplies by 1000 is correct as we want g/m2 at this stage. I've now fixed this in the nonanthropogenic version.
I don't think so. That was the way it was to correct from kg/m3 to g/cm3. But we are correcting in all the other charcoal equations from kg/m2 to g/cm2. That is a division by 10. I caught this because it produced charcoal 1000X higher for anthropogenic land use than for non-anthropogenic.
I just realized we still did it wrong, but the answer happily comes out OK by accident. Here is the equation for non-anthropogenic charcoal (simpler than anthropogenic, but the calculations are the same).
grass.mapcalc("${charcoal}=eval(biomass=graph(${lcov}, 0,0, 7,0.1, 18.5,0.66, 35,0.74, 50,1.95), pcntcharc=graph(${lcov}, 0,0, 5,0.0048, 18.5,0.0101, 50,0.0325), ((biomass pcntcharc 0.5) / (0.00011304 * 10)))", quiet=True, overwrite=True, lcov=lc, charcoal=charcoalmap, firemap=firemap)
The key part is: (biomass pcntcharc 0.5) / (0.00011304 * 10)
I assumed, incorrectly, that the goal was to end up with g/cm2, and incorrectly calculated the conversion factor as 1/10. But kg/m2 -> g/cm2 is 100cm x 100cm/1000gm = 10 not 1/10
But the goal is to come up with g/cell. So kg/m2 -> g/cell is 10m x 10m/1000m2 = 1/10 which is what we did for the 10x10m cells used. So we got it right but for the wrong reason.
Nonetheless, it would be better for different size cells to make the equation: (nsres() ewres() biomass pcntcharc 0.5) / (0.00011304 * 1000)
I should not multi-task (though I'm doing it again). Rethinking, my original correction for kg/m2->g/cm2 was correct = 1/10 kg/m2 = 1000g/10000cm2 = 1/10 g/cm2
My conversion to g/cell was likewise incorrect because the cell size is 5x5m. But it seems like the script already takes this into account in the accumulation loop.
There may be an issue in the phytolith calculations. But I'll start a new issue for that.
More working through the code. See if I understand.
After accumulating the initial proxy table, with one entry per model year in g/cm2, the loop for accumprox multiplies all values by 10 to get g/cm3 for 10cm levels. Then it adds all non eroded layers together into 10cm levels. I can't see where it is using the depth of each yearly layer in this, but I'm assuming it does somewhere.
So this all makes me think that phytoliths, starting out as g/cell, divided by cm2/cell to get g/cm2, and multiplied by 10 to get g/cm3 is right.
This means that charcoal is wrong because it starts out as g/cm2. The charcoal generating equation, like the phytolith generating equation should include multiplication by the cell size in meters, and conversion of kg to gm. Assuming that charcoal weight is in grams, this would make it:
(nsres()ewres() biomass pcntcharc 0.5 * 1000) / (0.00011304)
All current values need to be multiplied by 2500 = (5 5 1000/10)
I also just now realized that the current equation is calculating fragments/cm2 not grams. So to get a final grams/cell, it further needs to be changed to:
nsres()ewres() biomass pcntcharc 0.5 * 1000
To convert my current values, they need to be multiplied by 2500 and 0.00011304
I had a chat with Grant about the charcoal. It is as I suspected. The current code calculates large charcoal fragments per cm2. It should be
nsres()ewres() biomass pcntcharc 0.5 1000 # for farmed land or lightning fires or nsres()ewres() biomass pcntcharc 0.5 0.05 * 1000 # for grazed land
He also said that a better conversion from grams to fragments would be 0.000033493. We don't want to use it in the code, which is supposed to output g/cm3 in the end. But I will use it to convert my current results
Great, this makes sense and brings it in closer alignment with how the other proxies are calculated. I'll update the script on here accordingly.
It is producing good values too in the final analysis too.
The mapcalc statement for anthropogenic landscape didn't get corrected. Part still incorrectly transforms units by multiplying by 1000 instead of dividing by 10.