search

ebimodeling / ghgvc

Ecosystem Climate Regulation Services Calculator
http://www.ecosystemservicescalc.org
Other
48 stars 12 forks source link

Create biome map mask/ Selection of biomes map #59

Closed dlebauer closed 8 years ago

dlebauer commented 10 years ago

We currently have the following maps for biome-type (in netcdf/GCS/biomes):

Desert.nc                           TemperateScrubAndWoodland.nc
MarshAndSwampland.nc                TopicalForestAndPeatForest.nc
NorthernPeatlandAndBorealForest.nc  TropicalSavanna.nc
TemperateForest.nc                  Tundra.nc
TemperateGrassland.nc               Unassigned_Ecoregion.nc

Options

  1. Combine these into a single file with "Biome" as an integer value
  2. Use SynMap: http://www.sciencedirect.com/science/article/pii/S0034425706000514
  3. Use this combination of Hurt et al might work (it is on a 0.5 degree grid): ftp://nacp.ornl.gov/synthesis/2009/frescati/model_driver_spinup/global/land_use_change/Hurtt_SYNMAP_Global_HD_base1801_v1.0.nc.gz

    References

For reference, this is Hurtt's Land Use[ 1] map harmonized with the Synergetic Land Cover Map (SYNMAP)[2]". For reference on the combination of these see [3] .

  1. Hurtt GC, Frolking S, Fearon MG, Moore B, Shevliakova E, Malyshev S, Pacala SW, Houghton RA (2006), The underpinnings of land-use history: Three centuries of global gridded land-use transitions, wood-harvest activity, and resulting secondary lands, Global Change Biol., 12, 1-22.
  2. SynMap: Jung M, Henkel K, Herold M, Churkina G (2006), Exploiting synergies of global land cover products for carbon cycle modeling, Remote Sens. Environ., 101, 534-553.
  3. For details, see http://mstmipsynthesis.pbworks.com/w/page/47766343/Global_Spinup_LULC_Readme (either request access or request @dlebauer to access).
teixeirak commented 10 years ago

Sorry for the delayed response on this. I'm not sure I understand the technical challenge. I do think that rethinking biome classification is in order, but will require a bit of rearranging of some parameters, so it's not a very simple task. Lets talk about it today. I'll need to spend some time looking into it.

teixeirak commented 9 years ago

I need to work on this issue more, but wanted to modify this previous post because I realize that my previous posts here were off track.

1- I agree that SYNMAP is a good starting place. I think we should use it as the first cut for determining which ecosystem types may be present in a given location. 2-I think we'll need to complement SYNMAP with something else because its categories don't all map very well onto biome data that we have. One possibility may be FAO ecozones (which we have already procured for this project) because it is a fairly standard way of defining biomes-- used by IPCC and others. This gives the advantage of easy comparison with policy and existing values that could be used as defaults (e.g., IPCC 2006 guidance on national GHG inventories).

What I need to work on is figuring out the best biome categories for default values and how these can map to SYNMAP.

David, does this generally sounds like a good working plan?

ecofloristic zones

teixeirak commented 8 years ago

The following are instructions for the revised biomes mapping. The new vegetation types will be associated with revised default values, which I will provide. Please contact me if you have any questions.

The new biomes mapping scheme will draw upon several maps:

  1. SYNMAP* (including life form assemblage and leaf type/phenology) a. *primary map, used to predict most likely land cover type b. Reference: SynMap: Jung M, Henkel K, Herold M, Churkina G (2006), Exploiting synergies of global land cover products for carbon cycle modeling, Remote Sens. Environ., 101, 534-553. c. David has this map and it is available for download here: SynMap. I’m not sure if it’s currently in Box.
  2. Koppen climate zones a. Reference: Peel MC, Finlayson BL, Mcmahon TA (2007) Updated world map of the Köppen-Geiger climate classification. Hydrology and Earth System Sciences Discussions, 4, 439–473. b. available as netCDF from ORNL DAAC, available for download here: http://webmap.ornl.gov/wcsdown/wcsdown.jsp?dg_id=10012_1
  3. FAO ecoregions a. Shapefile in Box/ gcs/ Maps to be added to the calculator. b. Has not yet been converted to netCDF
  4. IBIS potential vegetation a. Already in the calculator; currently used to define natural biome types
  5. Ramankutty pasture & crop maps a. Already in the calculator; currently used to define regions of crop and pasture.
  6. (Maps for specific crop types- various) a. Already in the calculator; currently used to define regions of specific crop types.

Each of these maps will be associated to major vegetation types as indicated in the file VegTypeAssociations.csv. (in Box/gcs/Revised Biomes Mapping (2016)). Those vegetation types are as follows:

  1. broadleaf evergreen forest
  2. broadleaf deciduous forest
  3. mixed broadleaf-conifer forest
  4. conifer evergreen forest
  5. conifer deciduous forest
  6. savanna
  7. shrub-dominated
  8. grass-dominated 8b. pasture
  9. cropland
  10. barren / very sparse vegetation
  11. Other (specific crop types)

New biome default values will be described/ tracked in issue #64

potterzot commented 8 years ago

Can someone either upload SynMap to box or provide a link to it?

I don't see it on Box right now, and the download link I found through google just downloads a file saying that the service is unavailable.

Thanks.

ValentineHerr commented 8 years ago

I added it in the "Map to be added to Calculator" folder in the box. It is called Hurtt_SYNMAP_Global_HD_2010.nc Let me know if this is not what you needed. Valentine

potterzot commented 8 years ago

First thanks @ValentineHerr for adding the synmap!

Also - I'm not entirely sure I understand what needs to happen here. I'm not really familiar with biome data in general, but I think what I've parse together is the following:

  1. Use VegType_to_BiomeType.csv and VegTypeAssociations.csv to create a link from each map to a biome.
  2. Apply the link data set from (1) to (in descending order) [SynMap, Koppen, FAO, IBIS, Ramkutty, specific crop types].
  3. Outcome should be a single netcdf file with latitude and longitude as dimensions, and each biome as a binary variable (i.e., for a lat/lon pair, there are ~12 biome variables, each with a 0 or a 1).
  4. Modify the app so that it uses the new file to get all biomes with a 1 for the selected lat/lon pair.

Just a thought - this will potentially be slower, since (in R at least) you can only get one variable at a time, so we'll have to get the value for each variable as ~12 separate calls to the file. An alternative would be to set the biome as a dimension, and then for a given lat/lon we could query the netcdf file once and return a vector of 0/1 values that correspond to the biome dimension and create the list from that.

If I'm missing a step or if I have the configuration wrong please let me know. Thanks.

teixeirak commented 8 years ago

I can't advise on the best way to program this. Are my instructions on the desired functionality clear?

potterzot commented 8 years ago

Okay, as of commit 94bf9b5, a working example of biome loading is here. The calculation of said biomes doesn't work yet because various changes in the biome_defaults have to be taken into account in the calculator, but at least you can see the biome selection process.

It would probably be good to walk through this with someone so that we can address any potential errors. Also see the notes below. Needed (see major notes for more info):

  1. integer values that correspond to IBIS names in this csv file. Once I have this I can finish the biophysical_net logic.
  2. where does the map value (HWSD) come from? Once I have this I can finish the OM_SOM calculation.

Major Notes

Minor Notes

teixeirak commented 8 years ago

Thanks. I will review thoroughly as soon as possible and provide feedback as I go. If it would be helpful to walk through with me via phone or videoconference, I could be available most of the afternoon or evening starting in about half an hour.

Regarding the two major needs,

1- @ValentineHerr or I will provide ASAP.

2- HWSD refers to Harmonized World Soil Database (HWSD) SOC estimates to 1 m. depth, and was previously selectable under edit, initial storage of organic matter, Potential Soil Organic Matter loss. There are files for this in box under gcs>European Soil Data Center, but I am not sure if this folder contains the map in the calculator. Beyond this, @dlebauer will need to help.

@ValentineHerr created the FAO/gez file and may be able to help with that.

dlebauer commented 8 years ago

@potterzot

teixeirak commented 8 years ago

Besides IBIS, are the other rules for revised biome mapping fully implemented in this version? That is, is it ready for review?

potterzot commented 8 years ago

I believe so. Let me summarize how it works in the code:

  1. Use synmap value to lookup vegetation types (map_vegtypes.csv)
  2. Use koppen value and vegetation type from #1 to determine biome code (koppen_biomes.csv)
  3. If biome code is GX or APX (grass or pasture), use fao value to determine biome code (fao_biomes.csv)
  4. pull default values for the biome by biome code (biome_defaults.csv)
  5. Calculate OM_SOM based on vegetation type (if "Cropland", 0.43 hwsd, else 0.3hwsd)
  6. if IBIS value for each given vegetation type is 1 (from map_vegtypes.csv), then biophysical_net = latent + sensible. Else 0.
  7. If "Other" vegetation type is indicated, add specific logic for that vegtype (not currently implemented)

If any of those steps are wrong, just let me know and it's not worth reviewing. If they all look right, then it's ready for review, bearing in mind (1) the HWSD values aren't in yet, and (2) the IBIS values aren't in yet.

Nicholas Potter http://potterzot.com 9C73 3AAA 6D99 86F1 F8F4 601D 298F 8713 1BE6 5941

On Mon, Aug 22, 2016 at 10:10 AM, Kristina Anderson-Teixeira < notifications@github.com> wrote:

Besides IBIS, are the other rules for revised biome mapping fully implemented in this version? That is, is it ready for review?

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/ebimodeling/ghgvc/issues/59#issuecomment-241481841, or mute the thread https://github.com/notifications/unsubscribe-auth/AAdIbmqi1lLfj9ZrvcWcJBm2MIy2Y2Srks5qidgbgaJpZM4CY1kM .

teixeirak commented 8 years ago

Steps 1 and 5 require modification. 1) Use synmap, KOPPEN, FAO, IBIS, and RAMANKUTTY values to lookup vegetation types (map_vegtypes.csv). There will be multiple options. The use of multiple map layers increases the probability of listing all of the important potential vegetation types for that location. Note that Ramankutty lookup is already implemented. See issue #82.

5) Calculate OM_SOM based on vegetation type (for natural ecosystem types or pasture, if SYNMAP="Cropland", 0.43 * hwsd, else 0.3*hwsd. for crops, OM_SOM=0)

potterzot commented 8 years ago

Okay, awesome. So just to clarify, you'd like me to concatenate all of the vegetation types available from all of the synmap, koppen, fao, ibis, ramankutty values to create the list of vegetation types for which to lookup biomes?

Also - @dlebauer, I also don't know where the ramankutty map file is. So the things I need to fill the above are:

  1. the IBIS values added to map_vegtypes.csv,
  2. the RAMANKUTTY values added to same file
  3. the location/name and variable of the RAMANKUTTY file.
  4. the variable to use from the hwsd.nc file.

I think that's it.

potterzot commented 8 years ago

Also, I think there's maybe a copy/paste error here:

5) Calculate OM_SOM based on vegetation type (for natural ecosystem types or pasture, if SYNMAP="Cropland", 0.43 * hwsd, else 0.3*hwsd. for crops, OM_SOM=0)

Is the only difference the addition of for crops, OM_SOM = 0? If so was already the case so it's all set.

teixeirak commented 8 years ago

Yes, that is correct.

  1. I'll need @ValentineHerr to help with this.
  2. The Ramankutty maps give numerical values of % crop or pasture cover. Any non-zero value should trigger listing of crop or pasture.
  3. The Ramankutty maps are already used to define pasture and cropland area in the previous calculator version (issue #82). Files may be named something like global cropland or global pasture. I can't help with file location.
  4. The variable we want is soil organic carbon (SOC) in the top 1 m of soil. I don't know the variable name offhand.
teixeirak commented 8 years ago

On 5), the difference is the addition of "for natural ecosystem types or pasture, if..." and "for crops, ...". It sounds like it may all be correct already.

dlebauer commented 8 years ago

t_oc is the %SOM in the top 0-30 cm layer; s_oc is %SOM in the 30-70 cm layer

The variable you want SOC = t_oc * 0.3 + s_oc * 0.7

teixeirak commented 8 years ago

NO; we do not want %SOM, we want SOC in units of total mass (tons C ha-1 to 1 m). I believe the file has the variables SOCSUB and SOCTOP. So we want SOC=SOCSUB + SOCTOP. @ValentineHerr can look at the file to verify the variable name.

dlebauer commented 8 years ago

We must be looking at different files - the hwsd.nc file has %SOM as mentioned above and bulk density in units of kg / dm3. With appropriate conversions we can compute the SOC as tons C ha-1 to 1 m. Let me know if you want me to work out the conversion, but the first step will be:

t_oc * 0.3 * t_ref_bulk+ s_oc * 0.7 * s_ref_bulk

dlebauer commented 8 years ago

Percent cover for different crops is here:

netcdf/GCS/Crops/Brazil/Sugarcane/brazil_sugc_fractional_10yr_avg.nc netcdf/GCS/Crops/Brazil/Soybean/brazil_soyb_fractional_10yr_avg.nc netcdf/GCS/Crops/US/SpringWheat/fractioncover/fswh_2.7_us.0.5deg.nc netcdf/GCS/Crops/US/Corn/fractioncover/fcorn_2.7_us.0.5deg.nc netcdf/GCS/Crops/US/Soybean/fractioncover/fsoy_2.7_us.0.5deg.nc

I don't see cover for pasture

dlebauer commented 8 years ago

Actually I think pasture fractional cover is in

netcdf/GCS/Pasture2000_5min.nc

Does that sound right? Are these the Ramankutty files?

teixeirak commented 8 years ago

Regarding SOC, working out the conversion may be easiest. I believe @ValentineHerr has been using a file with units of t/ha (if not we need to fix our calculations, Valentine).

teixeirak commented 8 years ago

Pasture2000_5min.nc sounds right for the Ramankutty files. Is there a crops one too?

The other files you listed are specific crop types We also want generic cropland from Ramankutty.

Are brazil_sugc_fractional_10yr_avg.nc and brazil_soyb_fractional_10yr_avg.nc working? @dlebauer, this relates to a question I sent you this AM via email.

dlebauer commented 8 years ago

This is how to compute SOC from the hwsd.nc file:

X <- (% SOM in top * top bulk density * % of soil profile in top) + (% SOM in subsoil * subsoil bulk density * % of soil profile in subsoil) SOC = X kg / dm3 * (10 dm / m) * (1,000,000 dm2 / ha) * (1 ton / 1000 kg)

SOC <- (t_oc * 0.3 * t_ref_bulk+ s_oc * 0.7 * s_ref_bulk) * 10000
dlebauer commented 8 years ago

I think they are working, will have to defer to @potterzot

potterzot commented 8 years ago

Are brazil_sugc_fractional_10yr_avg.nc and brazil_soyb_fractional_10yr_avg.nc working? @dlebauer, this relates to a question I sent you this AM via email.

Short answer: yes. More details: Specific crops are still handled in the way they used to be. The way it is handled is:

  1. for given lat/lon, load all fractional crop values (including the above 2).
  2. for each file, if the value for the location >= 0.01, then include the relevant crop biome as an agroecosystem biome.

There is some additional specific logic that applies to certain crops, but that's the process. The new biome mapping only applies to native ecosystems and generic pasture and cropland codes: c(AP1, AP2, C1, C2).

teixeirak commented 8 years ago

What is the additional specific logic that applies to sugarcane and soy in Brazil? Earlier, there was trouble getting the sugarcane (and soy?) map to work, and a hack was created whereby these were assigned by state (not desirable).

potterzot commented 8 years ago

I don't know how the map values were created, so there could be something in that that assigns by state, but in terms of the logic after reading the map value, the specifics are:

BR Sugarcane

BR Soy

Actual code

if (!is.na(res$braz_sugarcane_num) & 
      res$braz_sugarcane_num > 0.01 & 
      res$braz_sugarcane_num < 110.0) {
    biome_data$agroecosystem_eco["BR_sugarcane"] <- name_indexed_ecosystems["BR sugarcane"]
    biome_data$agroecosystem_eco[["BR_sugarcane"]]$latent <- custom 
    biome_data$agroecosystem_eco[["BR_sugarcane"]]$sw_radiative_forcing <- custom 
  }
  if (res$braz_fractional_soybean_num == 1 & 
      !is.na(res$br_sugc_latent_heat_flux_diff)) {
    biome_data$agroecosystem_eco["BR_soy"] <- name_indexed_ecosystems["BR soy"]
    biome_data$agroecosystem_eco[["BR_soy"]]$latent <- custom 
    biome_data$agroecosystem_eco[["BR_soy"]]$sw_radiative_forcing <- custom 
  }
potterzot commented 8 years ago

Okay, a minimal working example is up. I believe it only works with selecting one natural ecosystem. Agriculture systems don't work at the moment because of the below.

Question: The format of the original biome defaults json file has source abbreviations like:

{..., F_N2O = {"s000" = ["0.02927426"]}, ...}

Because the new biome defaults are in a csv file, they have the format:

{..., F_N2O = ["0.02927426"], ...}

Essentially one level removed, with no "s000" key. This is a key for the source of the data from what I can tell. Right now this causes problems because the app will parse either one or the other, but not both. Do you want me to:

a. convert all to the format that includes the "s000" b. convert all to the format that does not include "s000" c. add logic to the app to allow for both formats?

dlebauer commented 8 years ago

I don't recognize "s000" notation - if it is not needed, please convert to the format that does not include it

potterzot commented 8 years ago

The demo is a complete working demo now, excepting the new issues that were created today and the remaining things below.

Here's the current list of stuff that needs to be done for this issue:

teixeirak commented 8 years ago

The soil carbon numbers are wrong (way too high); most values should be SOC= 50-250 t C ha-1. Then the variable that goes into the calculator is "vulnerable soil organic matter", OM_SOM=SOC/.58*0.3 (or 0.43 when SYNMAP = crops). The .58 converts from C to organic matter (in retrospect, a stupid convention in my code!) and the 0.3 or 0.43 defines the fraction of that SOM that is vulnerable. Also, you need to be sure to convert % to fraction before applying David's formula. So here's a revised formula:

SOC <- (t_oc/100 * 0.3 * t_ref_bulk+ s_oc/100 * 0.7 * s_ref_bulk) * 10000 OM_SOM<- SOC /0.58 * 0.3 (or 0.43).

teixeirak commented 8 years ago

Regarding this, " have not seen a list of biome names for each biome code like T1, T3, S1, etc... Right now I'm using vegetation types, but this means multiple different biomes will show up in the app as "Broadleaf Evergreen Forest", for example. If we have a list of names and codes for each biome, it would be easy to change this to reflect the specific biome name."... I intentionally chose to define forest type based only on the tree type; I think that makes it easier for users who wouldn't find the biome names very meaningful. So this item is fine.

teixeirak commented 8 years ago

Regarding this: "the Revised Biome Mapping docx file has this item: " There are a few biome types for which default values do not (yet) exist, as indicated by “NaN”. For these, the calculator should give an error message. Something like: “Insufficient data for this vegetation type in this location. Please enter values manually or de-select this ecosystem type." This is probably not something that it makes sense for me to do. Right now these values are coded as NA in R."... Would it be complicated to make an error message pop up if the user tries to check the box? The error message can read, "Insufficient data for this vegetation type in this location. Please enter values manually or select a different ecosystem type."

teixeirak commented 8 years ago

Regarding this, " there are not currently instructions for other "MAP" biome settings, for sw_radiative_forcing, latent, and sensible. These need to be implemented."…

Basically, these will be treated exactly the same as they were in the previous version, with the exception that they now apply only some of the time to the natural vegetation types (trees, savanna, shrubs, grassland). The basic logic is: if [selected ecosystem type]==[IBIS vegetation type at selected location], include the biophysical forcings (sw_radiative_forcing, latent). else, no data available.

These maps (or at least most of them) are already loaded in the calculator, and I do not know what they are called. Perhaps @dlebauer can help. The variable definitions and formulas below describe the calculations you will need:

R_net- net radiation or sw radiative forcing, in W/m2.
LE- latent heat, in W/m2

For both, we have paired maps for vegetation (e.g., Natural Vegetation and Natural Vegetation_Bare). There are maps for natural vegetation, maize, soy, sugarcane, switchgrass, and miscanthus.

∆R_net= R_natural vegetation - Rnet bare (will usually be a negative number) ∆LE= LE natural vegetation - LE bare

for both, the number that goes into the calculator is converted as follows: ∆R_net /51007200000 x 1000000000 ∆LE /51007200000 x 1000000000

I’m not familiar with all the inner workings of the calculator, so I hope this all makes sense. Let me know if you have any questions.

ValentineHerr commented 8 years ago

I added the IBIS values in map_vegtypes.csv. Let me know if that works.

ValentineHerr commented 8 years ago

About the Soil carbon units. Units seem to be in t ha-1 according to instructions that I have about the layers. Also the values stored don't look like percentages. The files I am talking about are:

dlebauer commented 8 years ago

@ValentineHerr I don't recall where those files came from, but I was referring to the file hwsd.nc produced by MsTMIP. That should be easier because it has all of the necessary data in a single file.

potterzot commented 8 years ago

Closing this issue, I think it is done.