/* mkcasastaticnc.c
* *
* Emilio 4/1/99
* Updated program.
* Emilio 2/11-24/99
* Wrote and debugged program.
* COMPILE AS: gccemu mkcasastaticnc
******************************************************************************/
/* STEPS
- read all the static maps: AMAZON BASIN MASK, elev, veg map, soil type map
(need a specified order for reading the variables: mask, elev, veg, soil)
- read in the 2 soil data tables into two arrays, using a structure
- step through each pixel (do entire continent in the window, not just within the
mask), get the soil ID at current pixel, then get the %sand and %clay for all
depths, store in an array.
** consider writing into the attributes the entire "legend", or table giving names
to each ID, for veg maps and soil types
*/
#include "emu.h"
const char * cpszUsage = {"usage: mkcasastaticnc \n\
.\n\
"};
/* local functions */
void Setup(void);
void ProcessAndOutput(void);
void WriteSMap(short v, char *name, char *lname, char *units, short sFillValue);
void AssignRootDepth(void);
short GetDepthNumber(int depth);
short GetProfileIndex(int soilID);
short GetSoilGroup(short pcurr);
/* local variables */
#define ZINTERVALS 13
int errnc;
int ncID, tID, varIDsand, varIDclay;
static size_t count[] = {0,0,0}, start[] = {0,0,0};
float *SoilTypeMap, *map, *maptmp;
short **sand, **clay;
static double Zvalues[] = {50.0,100.0,150.0,200.0,250.0,300.0,
400.0,500.0,600.0,700.0,800.0,900.0,1000.0};
typedef struct _soilprof
{
short soilID;
char soilType[10];
char soilGroup[10];
float depthbot[ZINTERVALS];
float clay[ZINTERVALS];
float sand[ZINTERVALS];
} SOILPROF;
SOILPROF soildb[300];
short TotProfiles = 0;
int main(int argc, char *argv[])
{
AssignString(&pszNCtitle, "EOS-Amazon Model static input data sets, Ver. 2.0");
AssignString(&pszNChist, "Created by EMU/mkcasastaticnc.c. Feb. 24, 1999");
Initialize(&argc, argv, NOMASK);
Setup();
ProcessAndOutput();
CleanUp();
return(0);
}
/* ================================================================= */
void Setup(void)
{
int p, z;
maptmp = alloc1d_f(0,nRows*nCols-1);
map = alloc1d_f(0,nRows*nCols-1);
sand = alloc2d_s(0,ZINTERVALS-1,0,nRows*nCols-1);
clay = alloc2d_s(0,ZINTERVALS-1,0,nRows*nCols-1);
for (p=0; p < 300; p++)
{
soildb[p].soilID = 0;
for (z=0; z < ZINTERVALS; z++)
{
soildb[p].depthbot[z] = -999.0;
soildb[p].clay[z] = -999.0;
soildb[p].sand[z] = -999.0;
}
}
/* CREATE OUTPUT NETCDF FILE */
SetupDefaultAxes(FALSE, FALSE);
FillAxisProp("z", "Soil Depth (bottom)", "cm", FALSE, 0.0, 1000.0, ZINTERVALS, Zvalues, TRUE);
count[0] = 1;
count[1] = nRows;
count[2] = nCols;
ncID = newNCsetup(&tID, pszOutputFile, pszNCtitle, pszNChist);
}
/* ================================================================= */
void ProcessAndOutput(void)
{
int nLineArgs, i, z, soilID, depth; int bug=0;
float claycurr, sandcurr;
short p=0, pcurr, sFillValue = -999;
char line[200], *VarArgs[20];
FILE *fin;
/* READ MAPS AND WRITE THEM INTO NC FILE */
WriteSMap(0, "amzmask", "Amazon Basin Mask", "ID", 0);
WriteSMap(1, "elev30", "GTOPO30 DEM", "m", 0);
WriteSMap(2, "landcover_d", "Defries 8km Land Cover Map", "ID", 0);
AssignRootDepth();
WriteSMap(-1, "rootdepth_d", "Rooting Depth from Defries Land Cover", "cm", 0);
WriteSMap(3, "faosoiltype", "FAO Soil Type", "ID", 0);
WriteSMap(4, "soiltype", "Soil Type, FAO+EMBRAPA", "ID", 0);
SoilTypeMap = map; /* SoilTypeMap will be a float pointer */
/* ====READ SOIL TEXTURE TABLES, FILL STRUCTURE ARRAY (BOTH INTO A SINGLE ARRAY)==== */
if ((fin = fopen("/oc3/am_arcdir/casa/soils2/faosoiltexture.csv","r")) == NULL)
nrerror("Couldn't open file %s for writing\n", "faosoiltexture.csv");
(void)fgetline(line, 200, fin); /* omit the first line (field headers) */
while (fgetline(line, 200, fin) != 0)
{
nLineArgs = GetLineElements(VarArgs, line, ',');
soilID = atoi(VarArgs[0]);
pcurr = GetProfileIndex(soilID);
if (pcurr != p)
{
p = pcurr;
strcpy(soildb[p].soilType, VarArgs[2]);
strcpy(soildb[p].soilGroup, VarArgs[3]);
}
depth = atoi(VarArgs[1]);
claycurr = atof(VarArgs[4]);
sandcurr = atof(VarArgs[5]);
z = GetDepthNumber(depth);
soildb[p].depthbot[z] = (float)depth;
soildb[p].clay[z] = claycurr;
soildb[p].sand[z] = sandcurr;
/* clear the allocated memory for VarArgs */
for (i=0; i<nLineArgs; i++)
free(VarArgs[i]);
}
fclose(fin);
/* fill in the missing horizons with data from the SoilGroup, or upper horizons */
/* THIS BLOCK WILL ONLY EXTRAPOLATE MISSING HORIZONS IN *EXISTING* PROFILES! IF THE WHOLE PROFILE
IS MISSING (the soil type), IT DOESN'T DO ANYTHING */
for (pcurr=0; pcurr < TotProfiles; pcurr++)
{
p = GetSoilGroup(pcurr);
for (z=0; z < ZINTERVALS; z++)
if ((int)soildb[pcurr].depthbot[z] == -999)
{
if ((p == pcurr) || ((int)soildb[p].depthbot[z] == -999))
{
/* if z == 0 when p == pcurr, we have a problem! */
soildb[pcurr].depthbot[z] = Zvalues[z];
soildb[pcurr].clay[z] = soildb[p].clay[z-1];
soildb[pcurr].sand[z] = soildb[p].sand[z-1];
}
else
{
soildb[pcurr].depthbot[z] = Zvalues[z];
soildb[pcurr].clay[z] = soildb[p].clay[z];
soildb[pcurr].sand[z] = soildb[p].sand[z];
}
}
}
if ((fin = fopen("/oc3/am_arcdir/casa/soils2/embrapasoiltexture.csv","r")) == NULL)
nrerror("Couldn't open file %s for writing\n", "faosoiltexture.csv");
(void)fgetline(line, 200, fin); /* omit the first line (field headers) */
while (fgetline(line, 200, fin) != 0)
{
nLineArgs = GetLineElements(VarArgs, line, ',');
soilID = atoi(VarArgs[0]);
pcurr = GetProfileIndex(soilID);
if (pcurr != p)
{
p = pcurr;
strcpy(soildb[p].soilType, VarArgs[2]);
strcpy(soildb[p].soilGroup, VarArgs[3]);
}
depth = atoi(VarArgs[1]);
claycurr = atof(VarArgs[4]);
sandcurr = atof(VarArgs[5]);
z = GetDepthNumber(depth);
soildb[p].depthbot[z] = (float)depth;
soildb[p].clay[z] = claycurr;
soildb[p].sand[z] = sandcurr;
/* clear the allocated memory for VarArgs */
for (i=0; i<nLineArgs; i++)
free(VarArgs[i]);
}
fclose(fin);
if (bDebug)
printf("\n**bug=%d: Past read soil texture tables block; TotProfiles=%d\n", bug++,TotProfiles);
/* ============== CREATE THE MAPS OF SAND AND CLAY FOR EACH HORIZON ============== */
for (i = 0; i < nRows*nCols; i++)
{
soilID = (int) SoilTypeMap[i];
p = GetProfileIndex(soilID);
for (z = 0; z < ZINTERVALS; z++)
{
if (bDebug && (i % 20000 == 0) && (z % 3 == 0))
printf("\n i=%d, soilID=%d, p=%d, z=%d, clay=%f", i,soilID,p,z,soildb[p].clay[z]);
/* for profile with index "p" in soildb array, get the depth "z";
BUT if data for that depth doesn't exist, use FillValue !!!?? (-999) */
/* Use a scale factor of 0.01, offset of 0 */
clay[z][i] = (soildb[p].depthbot[z] > 0.0) ? (short) floor(100.0 * soildb[p].clay[z] + 0.5) : sFillValue;
sand[z][i] = (soildb[p].depthbot[z] > 0.0) ? (short) floor(100.0 * soildb[p].sand[z] + 0.5) : sFillValue;
}
}
/* ==================== WRITE THE DATA INTO NETCDF VARIABLES ==================== */
clay">varIDclay = DefNCvar(ncID, "clay", EMU_SHORT, "xyz", "Percent Clay, to this depth", "%");
AddNCCompressionAttributes(ncID, varIDclay, 0.01, 0.0);
sand">varIDsand = DefNCvar(ncID, "sand", EMU_SHORT, "xyz", "Percent Sand, to this depth", "%");
AddNCCompressionAttributes(ncID, varIDsand, 0.01, 0.0);
for (z=0; z < ZINTERVALS; z++)
{
start[0] = z;
errnc = nc_put_vara_short(ncID, clay">varIDclay, start, count, clay[z]);
errnc = nc_put_vara_short(ncID, sand">varIDsand, start, count, sand[z]);
}
errnc = nc_close(ncID);
}
/* ================================================================= */
void ProcessCommandLineArgs(int * pArgc, char * argv[])
{
}
/* ================================================================= */
void LocalCleanUp(void)
{
free1d_f(maptmp,0,nRows*nCols-1);
free1d_f(map,0,nRows*nCols-1);
free2d_s(sand, 0,ZINTERVALS-1,0,nRows*nCols-1);
free2d_s(clay, 0,ZINTERVALS-1,0,nRows*nCols-1);
}
/* ================================================================= */
/* ================================================================= */
void WriteSMap(short v, char *name, char *lname, char *units, short sFillValue)
{
int varID;
int i, irev, col, row;
if (v >= 0)
{
ReadXY(maptmp, pstVars[v], NONE, NONE);
for (i = 0; i < nRows*nCols; i++)
{
irev = getRowRevCol(i, &col, &row);
map[irev] = maptmp[i];
}
}
varID = DefNCvar(ncID, name, EMU_SHORT, "xy", lname, units);
AddAttribute(ncID, "none", varID, "none", "_FillValue", EMU_SHORT, 1, &sFillValue);
errnc = nc_put_var_float(ncID, varID, map);
}
/* ================================================================= */
void AssignRootDepth(void)
{
int i, landcoverID;
for(i=0; i < nRows*nCols; i++)
{
landcoverID = (int) map[i];
switch(landcoverID)
{
case 0: /* water */
map[i] = 0.0;
break;
case 1: /* Evergreen Needleleaf Forests */
map[i] = 500.0;
break;
case 2: /* Evergreen Broadleaf Forests */
map[i] = 1000.0;
break;
case 3: /* Deciduous Needleleaf Forests ???? */
map[i] = 500.0;
break;
case 4: /* Deciduous Broadleaf Forests */
map[i] = 1000.0;
break;
case 5: /* Mixed Forests */
map[i] = 800.0;
break;
case 6: /* Woodlands */
map[i] = 700.0;
break;
case 7: /* Wooded Grasslands/Shrubs */
map[i] = 400.0;
break;
case 8: /* Closed Bushlands or Shrublands */
map[i] = 300.0;
break;
case 9: /* Open Shrublands */
map[i] = 200.0;
break;
case 10: /* Grasses */
map[i] = 200.0;
break;
case 11: /* Croplands */
map[i] = 100.0;
break;
case 12: /* Bare */
map[i] = 50.0;
break;
case 13: /* Mosses and Lichens */
map[i] = 50.0;
break;
default:
map[i] = 0.0;
break;
}
}
}
/* ================================================================= */
short GetDepthNumber(int depth)
{
short z=0;
while(depth != (int)Zvalues[z])
z++;
return(z);
}
/* ================================================================= */
short GetProfileIndex(int soilID)
{
int p=0;
/* look for soildb[].soilID already assigned */
while (soildb[p].soilID != 0)
{
if (soildb[p].soilID == soilID)
return(p);
p++;
}
/* new soildb[].soilID; also calculates number of unique profiles in dataset, TotProfiles */
if (soilID > 0)
{
TotProfiles++;
soildb[p].soilID = soilID;
return(p);
}
else
return(TotProfiles);
}
/* ================================================================= */
short GetSoilGroup(short pcurr)
{
short p=0;
while (!istype(soildb[p].soilType, soildb[pcurr].soilGroup) && ((int)strlen(soildb[p].soilType) != 1))
p++;
return(p);
}