/* masks.c
* *
* Emilio 12/10/2000
* Moved InitSiteIndex() prototype from emu.h to here, to make it a
* local function b/c it's used only here. In ReadMask(), changed
* name of szszCurrentFile to szCurrentFile. Made some layout changes.
* *
* Function file for all the C programs written to analyze binary and NetCDF
* input and output for modelling using the Environmental Modeling Utilities,
* or "EMU" package.
******************************************************************************/
#include "emu.h"
#define MASKIDSFILE "./maskids.cfg"
#define MASKIDSMAXLINLEN 300
#define M 7
#define NSTACK 50
/* functions and variables defined in function.c */
extern char * pszConfigFile;
extern int GetLine(FILE * fpin, char * pszLine, char * pszMarker, BOOL bFromStart);
/* LOCAL FUNCTIONS: Important "system" functions*/
void ReadMask(void);
/* LOCAL FUNCTIONS: "utilities" used to support the above functions */
int InitSiteIndex(int *ids, int maskid, int *idsLength);
void Parse(char * pszLine, int n);
void extractElement(char * linesegment, int n, int * pnArray);
void indexx(unsigned short n, int arr[], unsigned short indx[]);
/* ================================================================= */
BOOL masktest(int maskid)
{
int m=0;
/* When MASKTYPE is "rangefileID1" and the Aggregate ID config file hasn't been *
* read yet (ie, in the first masktest() call in ReadMask() ), MskSitesN == 0 *
* and the 2nd if-then block will be used instead of the third */
if (istype(szMaskType, "singleposval"))
return((maskid == nRegionID) ? TRUE:FALSE);
else if (istype(szMaskType, "rangeposval") || istype(szMaskType, "rangefileIDs")
|| (istype(szMaskType, "rangefileID1") && MskSitesN == 0))
return((maskid > 0) ? TRUE:FALSE);
else if (istype(szMaskType, "rangefileID1"))
{
for (m=0; m < Masks[0].nMaskids; m++)
if (maskid == Masks[0].maskids[m])
return(TRUE);
return(FALSE);
}
else
return(FALSE);
}
/* ================================================================= */
int siteindex(int *ids, int maskid)
{
int s;
/* look for ids[] already assigned */
for (s=0; s < FlMskSitesN; s++)
if (ids[s] == maskid)
return(s);
/* return NONE if maskid doesn't match an existing id */
return(NONE);
}
/* ================================================================= */
/* ======================= LOCAL FUNCTIONS ========================= */
/* ================================================================= */
/* Read in the mask file, write the value selected as the mask key to
* the global variable regionid;
* request "-1" from varsfile if no specific key is desired; */
void ReadMask(void)
{
FILE * fpc, * fpmask, * maskidsfile;
int nsite, n, i, m;
int *tmpSitesIDs, *tmpSitesCellcount;
double *tmpSitesArea;
int firstvalidindx=0, nsort;
unsigned short *sortindx;
char szLine[MASKIDSMAXLINLEN];
char szMaskPath[FILENAMELEN];
char szMaskDataType[10];
char szCurrentFile[FILENAMELEN];
BYTE *pbyMask;
SHORT *psMask;
unsigned short *pusMask;
if (!bQuiet) printf("\n---- Reading Mask in ReadMask() \n");
if ((fpc = fopen(pszConfigFile, "rt")) == NULL)
ExitError("Couldn't open file %s for reading in ReadMask\n", pszConfigFile);
if (GetLine(fpc, szLine, "# MASK", TRUE))
sscanf(szLine, "%s %s %d", szMaskPath, szMaskDataType, &nRegionID);
else
Warning("Couldn't find \"# MASK\" in file %s in ReadMask\n", pszConfigFile);
if (GetLine(fpc, szLine, "# MASKTYPE", TRUE))
sscanf(szLine, "%s", szMaskType);
/* szMaskType: singleposval OR rangeposval OR rangefileIDs OR rangefileID1 */
else
Warning("Couldn't find \"# MASKTYPE\" in file %s in ReadMask\n", pszConfigFile);
strcpy(szCurrentFile, szMaskPath);
if ((fpmask = fopen(szCurrentFile, "r")) == NULL)
ExitError("Couldn't open pMask file %s for reading\n", szCurrentFile);
pMask = alloc1d_l(0,nRows*nCols-1);
/* convert input file data type to long-type */
if (istype(szMaskDataType, "byte"))
{
pbyMask = alloc1d_uc(0,nRows*nCols-1);
if (fread(pbyMask, sizeof(BYTE), nRows*nCols, fpmask) != (size_t) (nRows*nCols))
ExitError("Could not do fread in %s \n", szCurrentFile);
for (i=0; i < nRows*nCols; i++)
pMask[i] = (LONG) pbyMask[i];
free1d_uc(pbyMask, 0,nRows*nCols-1);
}
else if (istype(szMaskDataType, "short"))
{
psMask = alloc1d_s(0,nRows*nCols-1);
if (fread(psMask, sizeof(SHORT), nRows*nCols, fpmask) != (size_t) (nRows*nCols))
ExitError("Could not do fread in %s \n", szCurrentFile);
for (i=0; i < nRows*nCols; i++)
pMask[i] = (LONG) psMask[i];
free1d_s(psMask, 0,nRows*nCols-1);
}
else if (istype(szMaskDataType, "ushort"))
{
pusMask = alloc1d_us(0,nRows*nCols-1);
if (fread(pusMask, sizeof(unsigned short), nRows*nCols, fpmask) != (size_t) (nRows*nCols))
ExitError("Could not do fread in %s \n", szCurrentFile);
for (i=0; i < nRows*nCols; i++)
pMask[i] = (LONG) pusMask[i];
free1d_us(pusMask, 0,nRows*nCols-1);
}
else if (istype(szMaskDataType, "long"))
{
if (fread(pMask, sizeof(LONG), nRows*nCols, fpmask) != (size_t) (nRows*nCols))
ExitError("Could not do fread in %s \n", szCurrentFile);
}
else
ExitError("Mask data type %s not recognized \n", szMaskDataType);
fclose(fpc);
fclose(fpmask);
/* ===== READ ALL UNIQUE IDS FROM MASK FILE INTO GLOBAL ARRAY SITEIDS ======= */
MskSitesN=0; FlMskSitesN=0;
/* allocate and initialize tmpSites* arrays */
tmpSitesIDs = alloc1d_i(1,SITESMAX);
tmpSitesCellcount = alloc1d_i(1,SITESMAX);
tmpSitesArea = alloc1d_d(1,SITESMAX);
for(n=1; n <= SITESMAX; n++)
{
tmpSitesIDs[n] = 0;
tmpSitesCellcount[n] = 0;
tmpSitesArea[n] = 0.0;
}
/* read the mask file and make initial (tmpSites*) assingments */
for (i = 0; i < nRows*nCols; i++)
{
if (masktest((int)pMask[i]))
{
n = InitSiteIndex(&tmpSitesIDs[1], (int)pMask[i], &FlMskSitesN);
tmpSitesCellcount[n+1]++;
tmpSitesArea[n+1] += cellarea(i);
}
}
if (bDebug) printf(" FlMskSitesN = %d\n", FlMskSitesN);
/* allocate and assign values to FlMskSitesIDs[], FlMskSitesCellcount[], *
* and FlMskSitesArea[]; also assign value to FlMskTotCellcount and FlMskTotArea */
/* sort tmpSitesIDs in ascending order by creating an index array */
sortindx = alloc1d_us(1,SITESMAX);
indexx(SITESMAX, tmpSitesIDs, sortindx);
while (tmpSitesIDs[sortindx[firstvalidindx]] == 0)
firstvalidindx++;
/* now allocate and populate the sorted arrays into FlMskSites* arrays */
FlMskSitesIDs = alloc1d_i(0,FlMskSitesN-1);
FlMskSitesCellcount = alloc1d_i(0,FlMskSitesN-1);
FlMskSitesArea = alloc1d_d(0,FlMskSitesN-1);
FlMskTotCellcount = 0;
FlMskTotArea = 0.0;
for (n=0; n < FlMskSitesN; n++)
{
nsort = sortindx[firstvalidindx+n];
FlMskSitesIDs[n] = tmpSitesIDs[nsort];
FlMskSitesCellcount[n] = tmpSitesCellcount[nsort];
FlMskSitesArea[n] = tmpSitesArea[nsort];
FlMskTotCellcount += FlMskSitesCellcount[n];
FlMskTotArea += FlMskSitesArea[n];
if (bDebug)
printf(" FlMskSitesIDs[%d]=%d, FlMskSitesCellcount[%d]=%d, FlMskSitesArea[%d]=%.1f \n",
n,FlMskSitesIDs[n],n,FlMskSitesCellcount[n],n,FlMskSitesArea[n]/1000000.0);
}
free1d_i(tmpSitesIDs,1,SITESMAX);
free1d_i(tmpSitesCellcount,1,SITESMAX);
free1d_d(tmpSitesArea,1,SITESMAX);
free1d_us(sortindx,1,SITESMAX);
if (bDebug) printf(" < alloc1d_MASK if-blocks \n");
/* ===== populate aggregate ID Masks struct array and Msk* variables ===== */
if (istype(szMaskType, "rangefileID"))
{
/* read the aggregate ID file */
if ((maskidsfile = fopen(MASKIDSFILE, "r")) == NULL)
nrerror("Couldn't open file %s for reading\n", MASKIDSFILE);
/* MskSitesN = number of IDs (aggregate or otherwise) in the pMask */
while(fgetline(szLine, MASKIDSMAXLINLEN, maskidsfile) != 0)
MskSitesN++;
rewind(maskidsfile);
/* allocate space for array Masks of size MskSitesN and type MASKS */
Masks = alloc1d_MASK(0,MskSitesN-1);
for (n=0; n < MskSitesN; n++)
{
/* read individual IDs for nth Aggregate ID, then allocate space for array *
* Masks[n].maskids; Then read IDs into that array for each Aggregate ID */
/* parse szLine, read and write id elements */
(void) fgetline(szLine, MASKIDSMAXLINLEN, maskidsfile);
Parse(szLine, n);
if (bDebug)
{
printf("\nIN rangefileIDs for loop, Masks[%d].ID=%d: ", n,Masks[n].ID);
for(m=0; m < Masks[n].nMaskids; m++)
printf(" .maskids[%d]=%d: ", m,Masks[n].maskids[m]);
}
Masks[n].cellcount = 0;
Masks[n].area = 0.0;
for(m=0; m < Masks[n].nMaskids; m++)
{
nsite = siteindex(FlMskSitesIDs, Masks[n].maskids[m]);
Masks[n].cellcount += FlMskSitesCellcount[nsite];
Masks[n].area += FlMskSitesArea[nsite];
}
} /* END n FOR-LOOP */
fclose(maskidsfile);
}
else if (istype(szMaskType, "rangeposval"))
{
/* allocate space for array Masks of size MskSitesN and type MASKS */
MskSitesN = FlMskSitesN;
Masks = alloc1d_MASK(0,MskSitesN-1);
for (n=0; n < MskSitesN; n++)
{
Masks[n].ID = FlMskSitesIDs[n];
Masks[n].nMaskids = 1;
Masks[n].maskids = alloc1d_i(0,0);
Masks[n].maskids[0] = FlMskSitesIDs[n];
Masks[n].cellcount = FlMskSitesCellcount[n];
Masks[n].area = FlMskSitesArea[n];
if (bDebug) printf("IN rangeposval for-loop, Masks[%d].ID=%d \n", n,Masks[n].ID);
}
}
if (bDebug) printf(" > alloc1d_MASK if-blocks; MskSitesN=%d \n", MskSitesN);
/* allocate memory and assign values to MskSitesIDs[]; *
* also assign value to MskTotCellcount and MskTotArea */
MskSitesIDs = alloc1d_i(0,MskSitesN-1);
MskTotCellcount = 0;
MskTotArea = 0.0;
for (n=0; n < MskSitesN; n++)
{
IDs">MskSitesIDs[n] = Masks[n].ID;
MskTotCellcount += Masks[n].cellcount;
MskTotArea += Masks[n].area;
}
if (bDebug)
printf(" > MskSitesIDs[n] for-block; MskSitesIDs[1]=%d, MskTotCellcount=%d, MskTotArea=%.1f \n",
MskSitesIDs[1], MskTotCellcount, MskTotArea/1000000.0);
}
/* ================================================================= */
/* ================================================================= */
int InitSiteIndex(int *ids, int maskid, int *idsLength)
{
int sindex=0;
/* look for ids[] already assigned */
while (ids[sindex] != 0)
{
if (ids[sindex] == maskid)
return(sindex);
sindex++;
}
/* new ids[]; also calculates number of unique IDs in mask file, FlMskSitesN */
(*idsLength)++;
ids[sindex] = maskid;
return(sindex);
}
/* ================================================================= */
void Parse(char * pszLine, int n)
{
char * pCh;
char s[MASKIDSMAXLINLEN], s1[MASKIDSMAXLINLEN];
int m, array[SITESMAX];
Masks[n].nMaskids = 0;
strcpy(s, pszLine);
pCh = strchr(s, ':'); /* check s for :, no error checking */
*pCh = '\0'; /* clip the string s at : */
sscanf(s, "%d", &Masks[n].ID);
strcpy(s1, pszLine);
pCh = strchr(s1, ':'); /* find : again */
pCh++; /* skip to next char */
strcpy(s, pCh); /* copy string into s from that point */
while ((pCh = strchr(s, ',')) != NULL)
{
*pCh = '\0';
pCh++;
extractElement(s, n, array);
strcpy(s, pCh); /* copy string into s from that point */
}
/* get the last entry - note that this will fail if there's no entry. */
extractElement(s, n, array);
/* now allocate space for Masks[n].maskids array,
* then write array elements from array[i] to Masks[n].maskids[i] */
Masks[n].maskids = alloc1d_i(0,Masks[n].nMaskids - 1);
for (m=0; m < Masks[n].nMaskids; m++)
Masks[n].maskids[m] = array[m];
}
/* ================================================================= */
void extractElement(char * linesegment, int n, int * pnArray)
{
char * pOper;
char s[NAMELEN], s1[NAMELEN];
int i, m;
int nID, nFirstID, nLastID;
strcpy(s, linesegment);
if ((pOper = strchr(s, '-')) != NULL) /* handle "-" operator */
{
*pOper = '\0';
pOper++;
sscanf(s, "%d", &nFirstID);
strcpy(s1, pOper);
sscanf(s1, "%d", &nLastID);
/* make id assignments to aggregate mask */
for (m=0; m < FlMskSitesN; m++)
if (FlMskSitesIDs[m] <= nLastID)
{
if (FlMskSitesIDs[m] >= nFirstID)
{
pnArray[Masks[n].nMaskids] = FlMskSitesIDs[m];
Masks[n].nMaskids++;
}
}
else
break;
}
else if ((pOper = strchr(s, '<')) != NULL) /* handle "<" operator (less or equal) */
{
*pOper = '\0';
pOper++;
strcpy(s1, pOper);
sscanf(s1, "%d", &nLastID);
/* make id assignments to aggregate mask */
for (m=0; m < FlMskSitesN; m++)
if (FlMskSitesIDs[m] <= nLastID)
{
pnArray[Masks[n].nMaskids] = FlMskSitesIDs[m];
Masks[n].nMaskids++;
}
else
break;
}
else if ((pOper = strchr(s, '>')) != NULL) /* handle ">" operator (greater or equal) */
{
*pOper = '\0';
pOper++;
strcpy(s1, pOper);
sscanf(s1, "%d", &nFirstID);
/* make id assignments to aggregate mask */
for (m=0; m < FlMskSitesN; m++)
if (FlMskSitesIDs[m] >= nFirstID)
{
pnArray[Masks[n].nMaskids] = FlMskSitesIDs[m];
Masks[n].nMaskids++;
}
}
else /* no special operator present, just a single number */
{
sscanf(s, "%d", &nID);
pnArray[Masks[n].nMaskids] = nID;
Masks[n].nMaskids++;
}
}
/* ================================================================== */
/* indexes an array arr[1..n], ie, outputs the array indx[1..n] such that
arr[indx[j]] is in ascending order for j=1,2,...,N. The input quantities
n and arr are not changed */
void indexx(unsigned short n, int arr[], unsigned short indx[])
{
unsigned short i,indxt,ir=n,j,k,l=1;
int jstack=0,*istack;
int a;
istack=alloc1d_i(1,NSTACK);
for (j=1;j<=n;j++) indx[j]=j;
for (;;) {
if (ir-l < M) {
for (j=l+1;j<=ir;j++) {
indxt=indx[j];
a=arr[indxt];
for (i=j-1;i>=1;i--) {
if (arr[indx[i]] <= a) break;
indx[i+1]=indx[i];
}
indx[i+1]=indxt;
}
if (jstack == 0) break;
ir=istack[jstack--];
l=istack[jstack--];
} else {
k=(l+ir) >> 1;
SWAP(indx[k],indx[l+1]);
if (arr[indx[l+1]] > arr[indx[ir]]) {
SWAP(indx[l+1],indx[ir])
}
if (arr[indx[l]] > arr[indx[ir]]) {
SWAP(indx[l],indx[ir])
}
if (arr[indx[l+1]] > arr[indx[l]]) {
SWAP(indx[l+1],indx[l])
}
i=l+1;
j=ir;
indxt=indx[l];
a=arr[indxt];
for (;;) {
do i++; while (arr[indx[i]] < a);
do j--; while (arr[indx[j]] > a);
if (j < i) break;
SWAP(indx[i],indx[j])
}
indx[l]=indx[j];
indx[j]=indxt;
jstack += 2;
if (jstack > NSTACK) nrerror("NSTACK too small in indexx.","");
if (ir-i+1 >= j-l) {
istack[jstack]=ir;
istack[jstack-1]=i;
ir=j-1;
} else {
istack[jstack]=j-1;
istack[jstack-1]=l;
l=i;
}
}
}
free1d_i(istack,1,NSTACK);
}