- A suggestion for creating your own VIC soil, vegetation and vegetation parameters for VIC
Nathalie Voisin - August 2006
Disclaimer: no guarantee whatsoever...
Note that all the steps described here required a Windows machine for the SoilProgram, and a unix/linux/etx machine in order to run the C-sell scripts. You also need awk, c, and f90.
Based on: Laura Bowling (Laura's readme)and Fengge Su (Fengge's readme) readme files, and my own experience.
This page is certainly and unfortunately not complete. Some digging is expected
( windows to dos and vice versa, code changes for your own use, format changes etc).
See also Laura's page.
1/ The Soil File
You need:
- Soilprogram in Global Soil Data Products CD-ROM (IGBP-DIS)
The CD is presently on my desk but will soon find a more official place, like on one of the shelves.
- Read the soilfile description in the VIC documentation web page here.
1.1 The Soil Program
Use the soil program the retrieve the sand, clay, bulk density, field capacity, Ksat, WaterN, WiltPoint, ThetaS
as described in the VIC documentation
Before using it, if you do not want to go through this process several times, think about the following:
o how many soil layers do you want? ( 2 or 3 in general)
o think about the corresponding depths. Usually layer 1 is 10cm, layer 2 and 3 will be calibrated
but give an estimate, say 30cm and 150cm for example.
o what are the boundaries of the mask file that you will be using. For example in the Zambeze basin,
xmin=18, xmax=36, ymin=-21, ymax=-9.
o know the resolution, say 0.0833 degrees
The SoilProgram runs a the Windows machine... sit to one, put the CD, etc ...
o Create map unit files ( .pop and .map)
Plan to do it for each soil layer. (I originally wanted to do
layers deeper than 150cm but too many data were missing deeper than 150 cm.)
soilprogram -> Generate -> Map Unit Data
In the box to the bottom right, enter the depth of the soil ( 0 -10cm, then 10-30cm,
then 30cm-150cm).In the scroll down, select the variables
(sand, clay, bulk density, field capacity, Ksat), some are in the secondary variables. Select (.pop)
AND (.map) .
o Create the surface data files (.srf) AND one *.doc file
- (optional step; do it for small basins else go to next bullet).
Extract an image file for the basin of interest. ( go to tools -> extract image, then choose
Maps ->YourContinent.img, choose output image name (basin.img), enter. Then choose the box to be
clipped ( (-9,18) and (-21,36), resolution is 5 minutes, so that I did not need to change it).
NOTE that F. Su had a bigger domain than I had and she chose not to change the resolution
here but to aggregate the parameters later, after the Soil Program step. If you need
the entire continent or world, you do not need to do this step.
- Then create data surface using the extracted image; go to Data -> create data surface ;
image file is the .img you just created or YourContinent.img,
map unit file are all the *.map files
that you created in the previous step. This will create the *.srf files. Do this step
as many times as you have *.map files. Use the Soil Data option
- Then repeat the same step (create a datasurface) for ONLY ONE variable only (say Clay10.map) using the IDRISI option this time, creating another clay10.img file and
most importantly clay10.doc (ascii) file that describes the *.srf files and
will be used afterwards.
The resulting files can be displayed by the soilprogram: soilprogram -> View Map -> Data Surface.
You are hopefully done with the soil program. The Soil program also offers products from the UMD but I haven't been into that.
See also the readme file coming from Laura, soil.readme . In particular, she chose to get clay, sand, silt and Buld Density from the soil program only, the rest is coming from a lookup table
and Van Genuchten relationships.
1.2 Disaggregate
Note that each .srf file needed to be treated for the undesired carriage return.
I copied each file to .txt, then vi , then :1,$s/.$//, then save :wq .
IMPORTANT: The following scripts and codes are ONLY VALID for a ending spatial resolution that is a multiple of 5 minute ( 1/8th of a degree for example, 0.5 degree, etc). If you are aiming at say a 0.3 degree resolution, use the SoilProgram aggregation step.
Do the following if you DID NOT change the resolution in the soilprogram before. See F. Su large domain case.
make_arc30_ll_data.scr is Fengge's script do dissagregate *.srf file from 5 minutes to 30 arc seconds.
Edit the script in order to make it for each variable and depth layer : foreach (BulkDens10 BulkDens30 BulkDens150 Sand10 Sand30 etc) with the names you gave to your *.srf files. You also need to give to right *.doc file created as explained above. Finally, edit for the number of rows and columns you expect to have in you "box".
output: lonlat30_* (Sand/Clay/BulkDens)
run_soil_aggr.scr will aggregate 30 arc seconds to 1/8 degree.
I did not need this one so you still may need to update it to your basin boundaries, hard coded names
and cell size
/**manually change NaN to some values**/ and need to download and compile aggr_soil_2vic.c
output: vic_plata_BuldDens SandClayBuld.txt
Note that Laura used Scalesrf.c to perform this step ( see her readme file).Your call at this point.
Do the following if your *.srf files have the desired resolution ( you ask the soil program to put it to the right resolution).
make_5min_ll_data.scr will create for each variable and depth a file with the lat lon
and the corresponding variable value.
In the script, you need to change manually to which resolution you want it .
1.3 Make SandClay.txt
Based on the amount of clay and sand, the soil texture will be chosen later ( texture based soil file ...). More on this later.
Here is what to do:
tail +2 Sand10.txt >! junk
tail +2 Clay10.txt >! junk2
paste junk junk2 > sandclay10.txt
and do it for each soil layer
1.4 Assign texture to your soil based on % sand and % clay
You need fortan90 for this one so either you have it or go to plane2.
The usda.triangle program calls a function that classifies soil in the USDA textural triangle using sand and clay %.
(Created by: aris gerakis)
Download the program here.
Unzip and recompile:
f90 -o usda.triangle INPOLY.F90 TRIANGLE.F90 WHAT_TEX.F90
usda.triangle
input: sandclay.txt (no "NaN")
output: by default this is soilclas.out . Change name before doing it for each soil layer
1.5 Assign hydraulic parameter to your soil texture
Create the hydraulic parameters according the "cosby.1984.csv" look up table.
You need run_awk_nv.scr ,texture2prop.awk, cosby.1984.csv
Edit run_awk_nv.scr to make sure that the path is correct ( and do it foreach layer).
Output is para_hydro.txt ( change name for each layer afterwards)
Colum: Ks, Porosity, Field Cap, WiltPoint, b
1.6 Retrieve lonlat_PT.txt
You need to get for each cell the annual average precipitation and if you plan on using
the energy balance mode, the annual mean temperature as well, alse can be kep at -9999.
At this point, I just used the 0.5 degree global soil file I had in hands and use the regrid program
(Symap) to downscale it to my desired resolution.
awk '{ if ($1==1) printf(" %.2f " , $26)}' soilfile >! T_tobegridded.txt
awk '{ if ($1==1) printf(" %.2f " , $49)}' soilfile >! P_tobegridded.txt
awk '{ if ($1==1) printf(" %.2f %.2f\n",$3,$4 )} soilfile >! station_list ( then add number of row in first line
Then regrid using the station_list, T_tobegridded.txt then P_tobegridded.txt, the mask file of your new basin
Paste on lonlat_PT.txt.
1.7 Finish up the soil file
Retrieve last parameters and put it in the right format:
a) At this point, paste all the soil layer file for each parameters together:
paste KSlayer1 KSlayer2 Kslayer3 >! Ks.txt, etc
As an example here is my script make_parametersfiles.scr.
One last note: if you do not want to have to put 0 and 1 in the first column later on, you can do it here.
I used the following scripts :clipSoilVariable.c,
clipSoilVariableClass.c, and
make_clipsoil.scr. Input files are the pasted files ( 3 columns in each parameters file.
As before, some edits have to be done for your personal use...
b) Two ways for b):
I personnaly retrieved many parameters from the SoilProgram (BulkDens, WiltPoint, sandclay,
Ks ( in cm/day)), the porosity from the look up table from Cosby, and the exponent c using the
get_van_GN.c ( using ThetaS and sandclay) to output van_GN.txt.
do appropriate changes in path and hard coded names
gcc -lm -o get_van_GN get_van_GN.c
Usage: get_van_GN sandclay1.txt sandclay2.txt sandclay3.txt latlon.txt van_GN.txt
Note: assume soil texture are the same in one column
But Laura only used the VanGenuchten parameters and only used ((BulkDens, sand, clay, silt?).
Depending on if you used SoilProgram only or the Van Genuchten parameters, you mainly just need to change the hardcoded name in create_soilf.c.
Are hardcoded in the code :
-the path to the files
-the number of cells in the basin
-the units of the soilprogram variable
-the depth of the three layers
gcc -o create_soilf create_soilf.c -lm
create_soilf lonlatelev.txt yoursoilfile
NOTE: the ARNOS PARAMS option in the global file should be set up to FALSE. Read Nijssen 2001 to see
how the ARNO params are used in the soil file and the VIC update log to see if this is still a default or not. Else find the calc_arno.scr script.
NOTE2:( T Bohn contribution)
At this point you may wish to apply the following two steps:
fix_resid.pl
- This
does the quality control on your soil file.
set_mois
t.pl -
This sets your initial soil moisture based on soicharacteristics.
These are both available in the Utility
Programs directory of the ftp area.
That's it! Remains to calibrate it ...
Back to top
Nathalie Voisin, last updated 06/2008