Nick Craig-Wood skrev:
Run
gdb python
Then type
run msgppscomp.py
at the gdb prompt.
When it crashes, type "bt" for a back trace. This will pinpoint
exactly where it crashes.
Hopefully that will make things clearer. If not post the output.
I will post the entire extention if it helps:
/* These must be defined and included before all other standard
libraries */
#define ACPG_PYMODULE_WITH_IMPORT_ARRAY
#include "acpg_arrayobject_wrap.h"
/************************************************************************/
#define NumberOfTiles 12
#define res 15
#define AreaRow 4860
#define AreaCol 3645
#define LenTot 17714700
#define TileTot 1476225
#define ResRow 324
#define ResCol 243
#define LenRes 78732
#define PHYS
"/data/proj/safworks/sheldon/pps_v11/import/AUX_data/remapped/physiography."
#include <stdio.h>
#include <stdlib.h>
#include <stdarg.h>
#include <string.h>
#include "read_vhlhdf.h"
#include "read_hlhdf.h"
#include "hlhdf.h"
#include <sys/types.h>
#include <dirent.h>
#include "vhlhdf.h"
#include <math.h>
#include "smhi_saf.h"
#include <stdarg.h>
#include <mcheck.h>
/* function declaration */
static int inithdf = 1; /* 1 om initHlHdf ej har anropats, 0 annars.
*/
/* Functions for C wrapper */
//static PyObject *ErrorObject;
static int createPythonObject(void);
static int readPythonObject(void);
PyMODINIT_FUNC initmsgppsarea(void); /* msgppsarea will be the name of
the SO file */
/* Copied from Sarah */
static int arrShortFromHdfNode(HL_NodeList *nodelist, char *nodename);
static int arrIntFromHdfNode(HL_NodeList *nodelist, char *nodename,
char *satp);
/* *********************** */
static int readHDFPPS(void);
static int readHDFMSG(void);
static int readHDFQual(void);
static int readHDFPHYS(void);
static int evaluateBIT(unsigned short int value, int bitnumber);
static int Statistics(void);
static int FillArea(void);
static int Percent(void);
static int InterpAvgF(float** array, float** array15km);
static int InterpAvgI(int** array, int** array15km);
static int Set2Zero(void);
static int Reset(void);
static int Bias(void);
static int LightCon(void);
static int PrintInt(int **array, int R, int C);
static int Memory(void);
static int MemoryFree1D(void);
static int MemoryFreeOrd(void);
static int MemoryFree15km(void);
static int MemoryFreeA(void);
/* Data structure */
static struct Region {
int ysize;
int xsize;
} reg;
static struct Variables {
int** dayA;
int** twilightA;
int** msgFCA;
int** ppsFCA;
int** msgcloudyA;
int** ppscloudyA;
float** lightconA;
float** bias100A;
float** bias75A;
float** bias50A;
float** bias25A;
int** valconcenA;
float** lightcon15km;
float** bias10015km;
float** bias7515km;
float** bias5015km;
float** bias2515km;
int** valconcen15km;
int** day;
int** twilight;
int** msgcloudy;
int** msgFC;
int** ppsFC;
int** ppscloudy;
int** msgdata;
int** ppsdata;
int** valconcen;
int** frac;
int* stats;
int Row;
int Col;
int Lentot;
int p1,p2,p3,p4;
int sumscenes;
int* nscenes;
int Cat;
unsigned short** qual;
int datein;
float* percentage;
char date[10];
char** tiles;
char** msg_scenes;
char** pps_scenes;
char tile[5];
char Region[20];
char Reprodir[200];
char PPSfile[200];
char MSGfile[200];
char PHYSfile[200];
PyObject* nscenesobj;
PyObject* tileobj;
PyObject* msgobj;
PyObject* ppsobj;
PyObject* outobj1;
PyObject* outobj2;
} work;
char *tileinfo[NumberOfTiles][6] = {
"05","ibla_68n29w","0","1215","0","1215",
"06","ibla_68n00e","0","1215","1215","2430",
"07","ibla_68n29e","0","1215","2430","3645",
"0B","ibla_57n20w","1215","2430","0","1215",
"0C","ibla_57n00e","1215","2430","1215","2430",
"0D","ibla_57n20e","1215","2430","2430","3645",
"13","ibla_46n16w","2430","3645","0","1215",
"14","ibla_46n00e","2430","3645","1215","2430",
"15","ibla_46n16e","2430","3645","2430","3645",
"1D","ibla_35n13w","3645","4860","0","1215",
"1E","ibla_35n00e","3645","4860","1215","2430",
"1F","ibla_35n13e","3645","4860","2430","3645"
};
PyObject *Rlightcon=NULL;
PyObject *Rbias100=NULL;
PyObject *Rbias75=NULL;
PyObject *Rbias50=NULL;
PyObject *Rbias25=NULL;
PyObject *Rvalcon=NULL;
PyObject *Rstats=NULL;
PyObject *Rpercentage=NULL;
static PyObject* msgppsarea(PyObject* self, PyObject* args, PyObject
*kwargs) { /* Start main function */
mtrace ();
mcheck(NULL);
static char *kwlist[] =
{"nscenesobj","tileobj","msgobj","ppsobj","dateint","sumscenes",NULL};
int i,j, xp;
int q,p;
int snum;
int counter;
int oldpos;
/*Take care of the input. Remember that the & sign is need for all
types of variables*/
if (!PyArg_ParseTupleAndKeywords(args,kwargs,"OOOOii:nothing",
kwlist,
&work.nscenesobj,
&work.tileobj,
&work.msgobj,
&work.ppsobj,
&work.datein,
&work.sumscenes)) {
printf("Error while receiving data in C\n");
goto fail;
}
/* Converts the date to a string */
sprintf(work.date,"%d",work.datein);
printf("Date: %s\n", work.date);
/* Allocate memory to python objects */
work.pps_scenes = malloc(sizeof *work.pps_scenes*work.sumscenes);
work.msg_scenes = malloc(sizeof *work.msg_scenes*work.sumscenes);
if (work.pps_scenes == NULL) {
printf("fail to allocate memory!");
exit(EXIT_FAILURE);
}
if (work.msg_scenes == NULL) {
printf("fail to allocate memory!");
exit(EXIT_FAILURE);
}
work.tiles = malloc(sizeof*work.tiles*NumberOfTiles);
for (i = 0; i < work.sumscenes; i++) {
work.pps_scenes
= malloc(sizeof 300);
work.msg_scenes = malloc(sizeof 300);
}
for (i = 0; i < NumberOfTiles; i++) {
work.tiles = malloc(sizeof 5);
}
work.nscenes = malloc(sizeof(int)*NumberOfTiles);
if (work.pps_scenes==NULL) {
printf("Failed to allocate memory for pps and msg scenes\n");
goto fail;
}
if (work.msg_scenes==NULL) {
printf("Failed to allocate memory for pps and msg scenes\n");
goto fail;
}
if (work.tiles==NULL) {
printf("Failed to allocate memory for pps and msg scenes\n");
goto fail;
}
if (work.nscenes==NULL) {
printf("Failed to allocate memory for pps and msg scenes\n");
goto fail;
}
work.Row = 1215;
work.Col = 1215;
reg.xsize = 1215;
reg.ysize = 1215;
work.Cat = 60;
/* Reading the Python lists into character arrays */
if (!readPythonObject()) {
printf("readPythonObject_int error\n");
goto fail;
}
if (!Memory()) {
printf("Failed to assign memory\n");
goto fail;
}
if (!(Set2Zero())) {
printf("Failed to initialize arrays\n");
goto fail;
}
snum = 0;
counter = 0;
oldpos = 0;
for (i = 0; i < NumberOfTiles; i++) {
snum += work.nscenes;
strcpy(work.tile,work.tiles);
/*printf("%d\n%d\n%s\n",work.Row,work.Col,work.Region);*/
/*printf("\nDoing tile: %s\n", work.tile);*/
for (xp = 0; xp < NumberOfTiles; xp++) {
if (strncmp(work.tile,tileinfo[xp][0],2) == 0) {
strcpy(work.Region,tileinfo[xp][1]);
work.p1 = atoi(tileinfo[xp][2]);
work.p2 = atoi(tileinfo[xp][3]);
work.p3 = atoi(tileinfo[xp][4]);
work.p4 = atoi(tileinfo[xp][5]);
}
}
strcpy(work.PHYSfile,PHYS);
strcat(work.PHYSfile,work.Region);
strcat(work.PHYSfile,".hdf");
printf("PHYSfile: %s\n", work.PHYSfile);
/* fetching the latlon data */
if (!(readHDFPHYS())) {
printf("Failed reading frac data\n");
goto fail;
}
//printf("\nsnum: %d\n", snum);
/* Popuplate the arrays */
printf("\nProcessing %d scenes for tile %s\n",
work.nscenes,work.tile);
for (q = oldpos; q < snum; q++) {
oldpos = snum;
/*printf("\n Processing scene number %d \n", ++counter);*/
/* printf("Copying strings\n"); */
strcpy(work.MSGfile,work.msg_scenes[q]);
/*printf("%s\n", work.MSGfile);*/
strcpy(work.PPSfile,work.pps_scenes[q]);
/*printf("%s\n", work.PPSfile);*/
/*printf("Copying complete\n"); */
/* Reading the data into the structure */
if (!(readHDFPPS())) {
printf("Failed reading pps data\n");
goto fail;
}
if (!(readHDFMSG())) {
printf("Failed reading msg data\n");
goto fail;
}
if (!(readHDFQual())) {
printf("Failed reading quality data\n");
goto fail;
}
/* Doing the statistics */
if (!(Statistics())) {
printf("Failed collecting the stats\n");
goto fail;
}
} /* End loop over scenes */
if (!(FillArea())) {
printf("Failed collecting the stats\n");
goto fail;
}
printf("Moving on to next tile\n");
/* Reinitialize arrays */
if (!(Reset())) {
printf("Failed to reinitialize arrays\n");
goto fail;
}
} /* End loop over tiles */
for (i = 0; i < work.sumscenes; i++) {
free(work.pps_scenes);
free(work.msg_scenes);
}
free(work.pps_scenes);
free(work.msg_scenes);
for (i = 0; i < NumberOfTiles; i++) {
free(work.tiles);
}
free(work.tiles);
free(work.nscenes);
if(!Percent()) {
printf("Failed calculate percentages\n");
goto fail;
}
/* for (i = 0; i < work.Cat; i++) {
printf("percentage[%d]: %f\n",i,work.percentage);
} */
if (!(Bias())) {
printf("Failed to calculate biases\n");
goto fail;
}
if (!(LightCon())) {
printf("Failed to calculate lightcon\n");
goto fail;
}
/* Create Python Objects */
/* Free up memory from the ordinary arrays */
if (!(MemoryFreeOrd())) {
printf("Allocated free memory for ordinary arrays\n");
goto fail;
}
work.lightcon15km = malloc(sizeof(float*)*ResRow);
work.bias10015km = malloc(sizeof(float*)*ResRow);
work.bias7515km = malloc(sizeof(float*)*ResRow);
work.bias5015km = malloc(sizeof(float*)*ResRow);
work.bias2515km = malloc(sizeof(float*)*ResRow);
work.valconcen15km = malloc(sizeof(int*)*ResRow);
for (i = 0; i < ResRow; i++) {
work.lightcon15km = malloc(sizeof(float)*ResCol);
work.bias10015km = malloc(sizeof(float)*ResCol);
work.bias7515km = malloc(sizeof(float)*ResCol);
work.bias5015km = malloc(sizeof(float)*ResCol);
work.bias2515km = malloc(sizeof(float)*ResCol);
work.valconcen15km = malloc(sizeof(int)*ResCol);
}
if ((work.lightcon15km)==NULL) {
printf("Failed to allocating lightcon memory\n");
exit(EXIT_FAILURE);
}
if ((work.bias10015km)==NULL) {
printf("Failed to allocating bias100 memory\n");
exit(EXIT_FAILURE);
}
if ((work.bias7515km)==NULL) {
printf("Failed to allocating bias75 memory\n");
exit(EXIT_FAILURE);
}
if ((work.bias5015km)==NULL) {
printf("Failed to allocating bias50 memory\n");
exit(EXIT_FAILURE);
}
if ((work.bias2515km)==NULL) {
printf("Failed to allocating bias25 memory\n");
exit(EXIT_FAILURE);
}
if ((work.valconcen15km)==NULL) {
printf("Failed to allocating val memory\n");
exit(EXIT_FAILURE);
}
/* Interpolating data array */
if(!InterpAvgF(work.lightconA,work.lightcon15km)) {
printf("Failed interpolation lightcon\n");
goto fail;
}
if (!InterpAvgI(work.valconcenA,work.valconcen15km)) {
printf("Failed interpolation val\n");
goto fail;
}
if(!InterpAvgF(work.bias100A,work.bias10015km)) {
printf("Failed interpolation b100\n");
goto fail;
}
if (!InterpAvgF(work.bias75A,work.bias7515km)) {
printf("Failed interpolation b75\n");
goto fail;
}
if(!InterpAvgF(work.bias50A,work.bias5015km)) {
printf("Failed interpolation b50\n");
goto fail;
}
if (!InterpAvgF(work.bias25A,work.bias2515km)) {
printf("Failed interpolation b25\n");
goto fail;
}
/* Free Area arrays */
if ((MemoryFreeA())) {
printf("Allocated free memory for area arrays\n");
goto fail;
}
/*printf("VA[%d]: %f\n",i, work.va_area);
printf("biasdata[%d]: %f\n",i, work.biasdata);
printf("Number[%d]: %f\n",i, work.number);*/
printf("Creating python objects\n");
Rlightcon=PyTuple_New(ResRow*ResCol);
Rbias100=PyTuple_New(ResRow*ResCol);
Rbias75=PyTuple_New(ResRow*ResCol);
Rbias50=PyTuple_New(ResRow*ResCol);
Rbias25=PyTuple_New(ResRow*ResCol);
Rvalcon=PyTuple_New(ResRow*ResCol);
Rstats=PyTuple_New(work.Cat);
Rpercentage=PyTuple_New(work.Cat);
if (!(createPythonObject())) {
printf("Failed to create python arrays objects\n");
goto fail;
}
printf("Returning data back to Python!\n");
muntrace ();
return Py_BuildValue("OOOOOOOO",Rbias100,
Rbias75,
Rbias50,
Rbias25,
Rlightcon,
Rvalcon,
Rstats,
Rpercentage);
fail:
exit(EXIT_FAILURE);
}
/* Functions for the C-wrapper */
/* This structure is needed for the python-C-interface */
static struct PyMethodDef msgppsareaMethods[] = {
{"msgppsarea", (PyCFunction)msgppsarea,
METH_VARARGS|METH_KEYWORDS,"HELP for msgppsarea\n"},
{NULL,NULL,0,NULL}
};
PyMODINIT_FUNC initmsgppsarea(void) {
(void) Py_InitModule("msgppsarea",msgppsareaMethods);
import_array(); /* access to Numeric PyArray functions */
}
/*Read an list of strings from a python object*/
static int readPythonObject(void) {
int i;
PyObject *msgop;
PyObject *ppsop;
PyObject *tileop;
PyObject *sceneop;
for (i = 0; i < work.sumscenes; i++) {
msgop = PyList_GetItem(work.msgobj, i);
work.msg_scenes = PyString_AsString(msgop);
ppsop = PyList_GetItem(work.ppsobj, i);
work.pps_scenes = PyString_AsString(ppsop);
}
for (i = 0; i < NumberOfTiles; i++) {
tileop = PyList_GetItem(work.tileobj, i);
work.tiles = PyString_AsString(tileop);
sceneop = PyList_GetItem(work.nscenesobj, i);
work.nscenes = PyInt_AsLong(sceneop);
}
return 1;
} /*end readPythonObject*/
static int createPythonObject(void) {
int i,j,k;
PyObject* op1;
PyObject* op2;
PyObject* ligop;
PyObject* valop;
PyObject* b100op;
PyObject* b75op;
PyObject* b50op;
PyObject* b25op;
for (i = 0; i < work.Cat; i++) {
op1 = PyInt_FromLong((long)work.stats);
op2 = PyFloat_FromDouble((double)work.percentage);
if (PyTuple_SetItem(Rstats,i,op1) !=0) {
printf("Error in creating python object Rstats\n");
exit(EXIT_FAILURE);
}
if (PyTuple_SetItem(Rpercentage,i,op2) !=0) {
printf("Error in creating python object Rpercentage\n");
exit(EXIT_FAILURE);
}
}
printf("Completed 1D arrays\n");
k=0;
for (i = 0; i < ResRow; i++) {
for (j = 0; j < ResCol; j++) {
ligop = PyFloat_FromDouble((double)work.lightcon15km[j]);
if (PyTuple_SetItem(Rlightcon,k,ligop) !=0) {
printf("Error in creating python light object\n");
exit(EXIT_FAILURE);
}
valop = PyInt_FromLong((long)work.valconcen15km[j]);
if (PyTuple_SetItem(Rvalcon,k,valop) !=0) {
printf("Error in creating python val object\n");
exit(EXIT_FAILURE);
}
b100op = PyFloat_FromDouble((double)work.bias10015km[j]);
if (PyTuple_SetItem(Rbias100,k,b100op) !=0) {
printf("Error in creating python bias100 object\n");
exit(EXIT_FAILURE);
}
b75op = PyFloat_FromDouble((double)work.bias7515km[j]);
if (PyTuple_SetItem(Rbias75,k,b75op) !=0) {
printf("Error in creating python bias75 object\n");
exit(EXIT_FAILURE);
}
b50op = PyFloat_FromDouble((double)work.bias5015km[j]);
if (PyTuple_SetItem(Rbias50,k,b50op) !=0) {
printf("Error in creating python bias50 object\n");
exit(EXIT_FAILURE);
}
b25op = PyFloat_FromDouble((double)work.bias2515km[j]);
if (PyTuple_SetItem(Rbias25,k,b25op) !=0) {
printf("Error in creating python bias25 object\n");
exit(EXIT_FAILURE);
}
k++;
}
}
printf("Completed 2D arrays.\nNow freeing 1D arrays\n");
if (!(MemoryFree1D())) {
printf("Couldn't free memory for area arrays\n");
goto fail;
}
printf("freed 1D arrays\n");
printf("freeing 15km arrays\n");
if (!(MemoryFree15km())) {
printf("Couldn't free memory for 15km arrays\n");
goto fail;
}
printf("Complete\n");
return 1;
fail:
exit(EXIT_FAILURE);
} /* end createPythonObject */
static int readHDFPHYS(void) {
HL_NodeList *nodelist = NULL; /* These are take from read_vlhdf.h
where some explanation is given */
char *field = "/fracofland";
if (inithdf) { /* Körs bara första gången rutinen anropas. Måste
vara med */
initHlHdf();
inithdf = 0;
printf("Finished setting inithdf\n");
}
if (!(nodelist=readHL_NodeList(work.PHYSfile))) {
printf("Error reading nodelist! from files: %s\n", work.PHYSfile);
goto fail;
}
/*printf("Checking dataset's existens\n"); */
if (!selectAllNodes(nodelist)) {
printf("Fail to select items\n");
goto fail;
}
if(!fetchMarkedNodes(nodelist)) {
printf("Failed to fetch data\n");
goto fail;
}
/* Reading cloudmask data into arrays */
if (!arrIntFromHdfNode(nodelist,field,"f\0")) {
printf("Failed getting fracofland data\n");
goto fail;
}
if (nodelist!=NULL) {
freeHL_NodeList(nodelist);
}
return 1;
fail:
if (nodelist!=NULL) {
freeHL_NodeList(nodelist);
}
exit(EXIT_FAILURE);
}
static int readHDFPPS(void) {
HL_NodeList *nodelist = NULL; /* These are take from read_vlhdf.h
where some explanation is given */
char *field = "/cloudmask";
if (inithdf) { /* Körs bara första gången rutinen anropas. Måste
vara med */
initHlHdf();
inithdf = 0;
printf("Finished setting inithdf\n");
}
if (!(nodelist=readHL_NodeList(work.PPSfile))) {
printf("Error reading nodelist! from files: %s\n", work.PPSfile);
goto fail;
}
/*printf("Checking dataset's existens\n"); */
if (!selectAllNodes(nodelist)) {
printf("Fail to select items\n");
goto fail;
}
if(!fetchMarkedNodes(nodelist)) {
printf("Failed to fetch data\n");
goto fail;
}
/* Reading cloudmask data into arrays */
if (!arrIntFromHdfNode(nodelist,field,"p\0")) {
printf("Failed getting the cloud mask\n");
goto fail;
}
if (nodelist!=NULL) {
freeHL_NodeList(nodelist);
}
return 1;
fail:
if (nodelist!=NULL) {
freeHL_NodeList(nodelist);
}
exit(EXIT_FAILURE);
}
static int readHDFMSG(void) {
HL_NodeList *nodelist = NULL; /* These are take from read_vlhdf.h
where some explanation is given */
char *field = "/cloudmask";
if (inithdf) { /* Körs bara första gången rutinen anropas. Måste
vara med */
initHlHdf();
inithdf = 0;
printf("Finished setting inithdf\n");
}
if (!(nodelist=readHL_NodeList(work.MSGfile))) {
printf("Error reading nodelist! from files: %s\n", work.MSGfile);
goto fail;
}
/*printf("Checking dataset's existens\n"); */
if (!selectAllNodes(nodelist)) {
printf("Fail to select items\n");
goto fail;
}
if(!fetchMarkedNodes(nodelist)) {
printf("Failed to fetch data\n");
goto fail;
}
/* Reading cloudmask data into arrays */
if (!arrIntFromHdfNode(nodelist,field,"m\0")) {
printf("Failed getting the cloud mask\n");
goto fail;
}
if (nodelist!=NULL) {
freeHL_NodeList(nodelist);
}
return 1;
fail:
if (nodelist!=NULL) {
freeHL_NodeList(nodelist);
}
exit(EXIT_FAILURE);
}
static int readHDFQual(void) {
HL_NodeList *nodelist = NULL; /* These are take from read_vlhdf.h
where some explanation is given */
char *field = "/quality_flag";
if (inithdf) { /* Körs bara första gången rutinen anropas. Måste
vara med */
initHlHdf();
inithdf = 0;
printf("Finished setting inithdf\n");
}
if (!(nodelist=readHL_NodeList(work.PPSfile))) {
printf("Error reading nodelist from files: %s\n", work.PPSfile);
goto fail;
}
if (!selectAllNodes(nodelist)) {
printf("Fail to select items\n");
goto fail;
}
if(!fetchMarkedNodes(nodelist)) {
printf("Failed to fetch data\n");
goto fail;
}
if (!(arrShortFromHdfNode(nodelist,field))) {
printf("Failed getting the quality flag\n");
goto fail;
}
if (nodelist!=NULL) {
freeHL_NodeList(nodelist);
}
return 1;
fail:
if (nodelist!=NULL) {
freeHL_NodeList(nodelist);
}
exit(EXIT_FAILURE);
}
/* Statistics function */
static int Statistics(void) { /* Start */
int clear = 1;
int cloudy = 3;
int i, j;
int island, iscoast, isnight, istwilight, isday;
/*printf("Calculating the stats\n");*/
for(i = 0; i < work.Row; i++) { /* Start loop */
for(j = 0; j < work.Col; j++) { /* Start loop */
if (work.ppsdata[j] == 0) work.ppsdata[j] = 255;
if (work.msgdata[j] == 0) work.msgdata[j] = 255;
if (work.ppsdata[j] == 5) work.ppsdata[j] = 255;
if (work.msgdata[j] == 5) work.msgdata[j] = 255;
if (work.ppsdata[j] == 4) work.ppsdata[j] = 1;
if (work.msgdata[j] == 4) work.msgdata[j] = 1;
if (work.msgdata[j] != 255 && work.ppsdata[j] != 255) {
work.valconcen[j]++;
if (work.ppsdata[j] == 2) {
work.ppsFC[j]++;
}
if (work.ppsdata[j] == 3) {
work.ppscloudy[j]++;
}
if (work.msgdata[j] == 2) {
work.msgFC[j]++;
}
if (work.msgdata[j] == 3) {
work.msgcloudy[j]++;
}
if (work.frac[j] > 0 && work.frac[j] < 255) {
iscoast = 1;
} else {
iscoast = 0;
}
island = evaluateBIT((unsigned short int)work.qual[j],0);
/*iscoast = evaluateBIT((unsigned short int)work.qual[j],1);
coast not set in DWD files */
isnight = evaluateBIT((unsigned short int)work.qual[j],2);
istwilight = evaluateBIT((unsigned short int)work.qual[j],3);
/* Test */
/*if (iscoast == 1) printf("Found coast\n"); */
/* Creating day logical and doing some stats */
if(isnight) {
isday = 0;
work.stats[5]++; /*6*/
} else if (istwilight) {
isday = 0;
work.twilight[j]++;
work.stats[10]++; /*11*/
} else {
isday = 1;
work.day[j]++;
work.stats[0]++; /*1*/
}
/* Clear conditions */
if (work.msgdata[j] == clear && work.ppsdata[j] == clear &&
isday) {
work.stats[1]++; /*2*/
} else if (work.msgdata[j] == clear && work.ppsdata[j] == clear
&& isnight) {
work.stats[6]++; /*7*/
} else if (work.msgdata[j] == clear && work.ppsdata[j] == clear
&& istwilight) {
work.stats[11]++; /*12*/
}
/* Cloudy conditions */
if (work.msgdata[j] == cloudy && work.ppsdata[j] == cloudy &&
isday) {
work.stats[2]++; /*3*/
} else if (work.msgdata[j] == cloudy && work.ppsdata[j] ==
cloudy && isnight) {
work.stats[7]++; /*8*/
} else if (work.msgdata[j] == cloudy && work.ppsdata[j] ==
cloudy && istwilight) {
work.stats[12]++; /*13*/
}
/* Clear work.msgdata & cloudy work.ppsdata */
if (work.msgdata[j] == clear && work.ppsdata[j] == cloudy &&
isday) {
work.stats[3]++; /*4*/
} else if (work.msgdata[j] == clear && work.ppsdata[j] == cloudy
&& isnight) {
work.stats[8]++; /*9*/
} else if (work.msgdata[j] == clear && work.ppsdata[j] == cloudy
&& istwilight) {
work.stats[13]++; /*14*/
}
/* Cloudy work.msgdata & clear pps */
if (work.msgdata[j] == cloudy && work.ppsdata[j] == clear &&
isday) {
work.stats[4]++; /*5*/
} else if (work.msgdata[j] == cloudy && work.ppsdata[j] == clear
&& isnight) {
work.stats[9]++; /*10*/
} else if (work.msgdata[j] == cloudy && work.ppsdata[j] == clear
&& istwilight) {
work.stats[14]++; /*15*/
}
/* Counting the pixels over land */
if (isday && island) {
work.stats[15]++; /*16*/
} else if (isnight && island) {
work.stats[20]++; /*21*/
} else if (istwilight && island) {
work.stats[25]++; /*26*/
}
/* Land clear */
if (work.msgdata[j] == clear && work.ppsdata[j] == clear &&
isday && island) {
work.stats[16]++; /*17*/
} else if (work.msgdata[j] == clear && work.ppsdata[j] == clear
&& isnight && island) {
work.stats[21]++; /*22*/
} else if (work.msgdata[j] == clear && work.ppsdata[j] == clear
&& istwilight && island) {
work.stats[26]++; /*27*/
}
/* Land cloudy */
if (work.msgdata[j] == cloudy && work.ppsdata[j] == cloudy &&
isday && island) {
work.stats[17]++; /*18*/
} else if (work.msgdata[j] == cloudy && work.ppsdata[j] ==
cloudy && isnight && island) {
work.stats[22]++; /*23*/
} else if (work.msgdata[j] == cloudy && work.ppsdata[j] ==
cloudy && istwilight && island) {
work.stats[27]++; /*28*/
}
/* work.msgdata clear pps cloudy */
if (work.msgdata[j] == clear && work.ppsdata[j] == cloudy &&
isday && island) {
work.stats[18]++; /*19*/
} else if (work.msgdata[j] == clear && work.ppsdata[j] == cloudy
&& isnight && island) {
work.stats[23]++; /*24*/
} else if (work.msgdata[j] == clear && work.ppsdata[j] == cloudy
&& istwilight && island) {
work.stats[28]++; /*29*/
}
/* work.msgdata cloudy pps clear */
if (work.msgdata[j] == cloudy && work.ppsdata[j] == clear &&
isday && island) {
work.stats[19]++; /*20*/
} else if (work.msgdata[j] == cloudy && work.ppsdata[j] == clear
&& isnight && island) {
work.stats[24]++; /*25*/
} else if (work.msgdata[j] == cloudy && work.ppsdata[j] == clear
&& istwilight && island) {
work.stats[29]++; /*30*/
}
/* Counting pixels over water */
if (isday && !island) {
work.stats[30]++; /*31*/
} else if (isnight && !island) {
work.stats[35]++; /*36*/
} else if (istwilight && !island) {
work.stats[40]++; /*41*/
}
/* both clear */
if (work.msgdata[j] == clear && work.ppsdata[j] == clear &&
isday && !island) {
work.stats[31]++; /*32*/
} else if (work.msgdata[j] == clear && work.ppsdata[j] == clear
&& isnight && !island) {
work.stats[36]++; /*37*/
} else if (work.msgdata[j] == clear && work.ppsdata[j] == clear
&& istwilight && !island) {
work.stats[41]++; /*42*/
}
/* both cloudy */
if (work.msgdata[j] == cloudy && work.ppsdata[j] == cloudy &&
isday && !island) {
work.stats[32]++; /*33*/
} else if (work.msgdata[j] == cloudy && work.ppsdata[j] ==
cloudy && isnight && !island) {
work.stats[37]++; /*38*/
} else if (work.msgdata[j] == cloudy && work.ppsdata[j] ==
cloudy && istwilight && !island) {
work.stats[42]++; /*43*/
}
/*work.msgdata clear pps cloudy */
if (work.msgdata[j] == clear && work.ppsdata[j] == cloudy &&
isday && !island) {
work.stats[33]++; /*34*/
} else if (work.msgdata[j] == clear && work.ppsdata[j] == cloudy
&& isnight && !island) {
work.stats[38]++; /*39*/
} else if (work.msgdata[j] == clear && work.ppsdata[j] == cloudy
&& istwilight && !island) {
work.stats[43]++; /*44*/
}
/* work.msgdata cloudy pps clear */
if (work.msgdata[j] == cloudy && work.ppsdata[j] == clear &&
isday && !island) {
work.stats[34]++; /*35*/
} else if (work.msgdata[j] == cloudy && work.ppsdata[j] == clear
&& isnight && !island) {
work.stats[39]++; /*40*/
} else if (work.msgdata[j] == cloudy && work.ppsdata[j] == clear
&& istwilight && !island) {
work.stats[44]++; /*45*/
}
/* Coastal effects */
if (isday && iscoast) {
work.stats[45]++; /*46 Number of valid pixs*/
} else if (isnight && iscoast) {
work.stats[46]++; /* 47 Number of valid pixs */
} else if (istwilight && iscoast) {
work.stats[47]++; /* 48 Number of valid pixs */
}
/* Clear conditions */
if (work.msgdata[j] == clear && work.ppsdata[j] == clear &&
isday && iscoast) {
work.stats[48]++; /*49*/
} else if (work.msgdata[j] == clear && work.ppsdata[j] == clear
&& isnight && iscoast) {
work.stats[49]++; /* 50*/
} else if (work.msgdata[j] == clear && work.ppsdata[j] == clear
&& istwilight && iscoast) {
work.stats[50]++; /*51*/
}
/* Cloudy conditions */
if (work.msgdata[j] == cloudy && work.ppsdata[j] == cloudy &&
isday && iscoast) {
work.stats[51]++; /*52*/
} else if (work.msgdata[j] == cloudy && work.ppsdata[j] ==
cloudy && isnight && iscoast) {
work.stats[52]++; /*53*/
} else if (work.msgdata[j] == cloudy && work.ppsdata[j] ==
cloudy && istwilight && iscoast) {
work.stats[53]++; /*54*/
}
/* Mixed conditions */
if (work.msgdata[j] == cloudy && work.ppsdata[j] == clear &&
isday && iscoast) {
work.stats[54]++; /*55*/
} else if (work.msgdata[j] == cloudy && work.ppsdata[j] == clear
&& isnight && iscoast) {
work.stats[55]++; /*56*/
} else if (work.msgdata[j] == cloudy && work.ppsdata[j] == clear
&& istwilight && iscoast) {
work.stats[56]++; /*57*/
}
if (work.msgdata[j] == clear && work.ppsdata[j] == clear &&
isday && iscoast) {
work.stats[57]++; /*58*/
} else if (work.msgdata[j] == clear && work.ppsdata[j] == clear
&& isnight && iscoast) {
work.stats[58]++; /*59*/
} else if (work.msgdata[j] == clear && work.ppsdata[j] == clear
&& istwilight && iscoast) {
work.stats[59]++; /* 60 */
}
}
}
}
/*for (i = 0; i < work.Cat; i++) {
printf("Stats[%d]: %d\t", i, work.stats);
}*/
return 1;
}
static int evaluateBIT(unsigned short int value, int bitnumber) {
int one=1;
return ( one & (value >> bitnumber) ); /* This does a & logical
calculation after the bitwise shift */
}
/*Read an array of integers, unknow type (i.e. int or ushort or uchar),
from node.
Give it as output in form of an array of integers*/
static int arrIntFromHdfNode(HL_NodeList *nodelist, char *nodename,
char *satp) {
HL_Node *node;
unsigned char *hdfdata=NULL;
int ii,i,j;
/*Read from node*/
if(!(node=getNode(nodelist,nodename))){
printf("Failed reading node %s\n",nodename);
goto fail;
}
/*if (node->format == "uchar or U8LE") {*/
/*put data first in array of "unsigned char", then in array of int*/
if ((hdfdata = malloc(sizeof(int)*TileTot))==NULL) {
printf("Failed allocating data\n");
goto fail;
}
/*put data first in array of "unsigned char", then in array of int*/
memcpy(hdfdata,node->data,node->dSize*TileTot);
ii = 0;
for (i = 0; i < work.Row; i++) {
for (j = 0; j < work.Col; j++) {
if ((strcmp(satp,"m\0")) == 0) {
work.msgdata[j] = (int)hdfdata[ii];
} else if ((strcmp(satp,"p\0")) == 0) {
work.ppsdata[j] = (int)hdfdata[ii];
} else if ((strcmp(satp,"f\0")) == 0) {
work.frac[j] = (int)hdfdata[ii];
}
ii++;
}
}
if (hdfdata!=NULL)
free(hdfdata);
return 1;
fail:
exit(EXIT_FAILURE);
} /*end arrIntFromHdfNode*/
/*Read an array of short from node.
Give it as output in form of an array of short */
static int arrShortFromHdfNode(HL_NodeList *nodelist, char *nodename) {
HL_Node *node;
unsigned short *hdfdata=NULL;
int ii;
int i,j;
/*Read from node*/
if(!(node=getNode(nodelist,nodename))) {
printf("Failed getting nodename %s\n",nodename);
goto fail;
}
/*else if (node->format == "ushort") {*/
/*put data first in array of "unsigned short", then in array of int*/
if ((hdfdata = malloc(sizeof(unsigned short)*TileTot))==NULL) {
printf("Failed allocating short array data\n");
goto fail;
}
/*printf("dSize short: %d", node->dSize);*/
memcpy(hdfdata,node->data,node->dSize*TileTot);
ii = 0;
for (i = 0; i < work.Row; i++) {
for (j = 0; j < work.Col; j++) {
work.qual[j] = (unsigned short)hdfdata[ii];
ii++;
}
}
if(hdfdata!=NULL)
free(hdfdata);
return 1;
fail:
exit(EXIT_FAILURE);
} /*end arrShortFromHdfNode*/
static int FillArea(void) {
/* Function to create 2D array and fill in with data */
int i,j,x,y;
x=0;
y=0;
for (i = work.p1; i < work.p2; i++) {
for (j = work.p3; j < work.p4; j++) {
work.valconcenA[j] = work.valconcen[x][y];
work.msgFCA[j] = work.msgFC[x][y];
work.ppsFCA[j] = work.ppsFC[x][y];
work.ppscloudyA[j] = work.ppscloudy[x][y];
work.msgcloudyA[j] = work.msgcloudy[x][y];
work.dayA[j] = work.day[x][y];
work.twilightA[j] = work.twilight[x][y];
y++;
}
y=0;
x++;
}
return 1;
}
/* Calculating the percentages */
static int Percent(void) {
int i;
printf("Calculating the percentage\n");
for (i = 0; i < work.Cat; i++) {
if (i > 0 && i < 5) { /* Day only */
if (work.stats[0] > 0) {
work.percentage = (float)((float)work.stats /
(float)work.stats[0] * 100.0);
if (work.percentage > 100.0) {
printf("Day error!\n");
goto error;
}
}
}
if (i > 5 && i < 10) { /* Night only */
if (work.stats[5] > 0) {
work.percentage = (float)((float)work.stats /
(float)work.stats[5] * 100.0);
if (work.percentage > 100.0) {
printf("Night error!\n");
goto error;
}
}
}
if (i > 10 && i < 15) { /* Twilight only */
if (work.stats[10] > 0) {
work.percentage = (float)((float)work.stats /
(float)work.stats[10] * 100.0);
if (work.percentage > 100.0) {
printf("Twilight error!\n");
goto error;
}
}
}
if (i > 15 && i < 20) { /* Day land */
if (work.stats[15] > 0) {
work.percentage = (float)((float)work.stats /
(float)work.stats[15] * 100.0);
if (work.percentage > 100.0) {
printf("Day land error!\n");
goto error;
}
}
}
if (i > 20 && i < 25) { /* Night land */
if (work.stats[20] > 0) {
work.percentage = (float)((float)work.stats /
(float)work.stats[20] * 100.0);
if (work.percentage > 100.0) {
printf("Night land error!\n");
goto error;
}
}
}
if (i > 25 && i < 30) { /* Twilight land */
if (work.stats[25] > 0) {
work.percentage = (float)((float)work.stats /
(float)work.stats[25] * 100.0);
if (work.percentage > 100.0) {
printf("Twilight land error!\n");
goto error;
}
}
}
if (i > 30 && i < 35) { /* Water day */
if (work.stats[30] > 0) {
work.percentage = (float)((float)work.stats /
(float)work.stats[30] * 100.0);
if (work.percentage > 100.0) {
printf("Water day error!\n");
goto error;
}
}
}
if (i > 35 && i < 40) { /* Water night */
if (work.stats[35] != 0) {
work.percentage = (float)((float)work.stats /
(float)work.stats[35] * 100.0);
if (work.percentage > 100.0) {
printf("Water night error!\n");
goto error;
}
}
}
if (i > 40 && i < 45) {
if (work.stats[40] != 0) { /* Twilight water */
work.percentage = (float)((float)work.stats /
(float)work.stats[40] * 100.0);
if (work.percentage > 100.0) {
printf("Twilight water error!\n");
goto error;
}
}
}
if (work.stats[45] > 0) { /*day coast*/
if (i == 48) work.percentage = (float)((float)work.stats /
(float)work.stats[45] * 100.0);
if (i == 51) work.percentage = (float)((float)work.stats /
(float)work.stats[45] * 100.0);
if (i == 54) work.percentage = (float)((float)work.stats /
(float)work.stats[45] * 100.0);
if (i == 57) work.percentage = (float)((float)work.stats /
(float)work.stats[45] * 100.0);
if (work.percentage > 100.0) {
printf("Day coast error!\n");
goto error;
}
}
if (work.stats[46] > 0) { /*night coast */
if (i == 49) work.percentage = (float)((float)work.stats /
(float)work.stats[46] * 100.0);
if (i == 52) work.percentage = (float)((float)work.stats /
(float)work.stats[46] * 100.0);
if (i == 55) work.percentage = (float)((float)work.stats /
(float)work.stats[46] * 100.0);
if (i == 58) work.percentage = (float)((float)work.stats /
(float)work.stats[46] * 100.0);
if (work.percentage > 100.0) {
printf("Night coast error!\n");
goto error;
}
}
if (work.stats[47] > 0) { /*twilight coast */
if (i == 50) work.percentage = (float)((float)work.stats /
(float)work.stats[47] * 100.0);
if (i == 53) work.percentage = (float)((float)work.stats /
(float)work.stats[47] * 100.0);
if (i == 56) work.percentage = (float)((float)work.stats /
(float)work.stats[47] * 100.0);
if (i == 59) work.percentage = (float)((float)work.stats /
(float)work.stats[47] * 100.0);
if (work.percentage > 100.0) {
printf("Twilight coast error!\n");
goto error;
}
}
/*printf("At end of loop: percent[%d]: %f\t",i,work.percentage);
*/
}
return 1;
error:
exit(EXIT_FAILURE);
}
/* Interpolate array to lower resolution */
/* Reduces the resolution of the 1km array to res km by
taking the average of all valid pixels in the res x res window. */
static int InterpAvgF(float** array, float** array15km) {
int i, j, r, k;
float avg;
float num;
printf("Calling InterpAvgF\n");
/* Slides over the window determined by res */
for (r = 0; r < ResRow; r++) { /* Start loop r */
for (k = 0; k < ResCol; k++) { /* Start loop k */
avg = 0.0;
num = 0.0;
for (i = (r*res); i < r*res+res; i++) {
for (j = (k*res); j < k*res+res; j++) {
if (array[j] != 255.0) {
avg += array[j];
num++;
}
}
}
if (num > 0) {
array15km[r][k] = avg / num;
} else {
array15km[r][k] = 255.0;
}
}
}
return 1;
}
static int InterpAvgI(int** array, int** array15km) {
int i, j, r, k;
int avg;
int num;
printf("Calling InterpAvgI\n");
/* Slides over the window determined by res */
for (r = 0; r < ResRow; r++) { /* Start loop r */
for (k = 0; k < ResCol; k++) { /* Start loop k */
avg = 0;
num = 0;
for (i = r*res; i < r*res+res; i++) {
for (j = k*res; j < k*res+res; j++) {
//printf("Val[%d][%d]: %d", i,j,array[j]);
if (array[j] > 0) {
avg += array[j];
num++;
}
}
}
if (num > 0) {
//printf("avg: %d\tnum: %d\n");
//printf("avg: %d", (int)((float)avg / (float)num));
array15km[r][k] = (int)( ((float)avg / (float)num) +.5);
} else {
array15km[r][k] = 255;
}
}
}
return 1;
}
static int Set2Zero(void) {
int i,j;
printf("Initializing all arrays\n");
for (i = 0; i < work.Row; i++) {
for (j = 0; j < work.Col; j++) {
work.day[j] = 0;
work.twilight[j] = 0;
work.ppsFC[j] = 0;
work.msgFC[j] = 0;
work.msgcloudy[j] = 0;
work.ppscloudy[j] = 0;
work.valconcen[j] = 0;
work.frac[j] = 0;