Help | Site Map
Connecting Tech Pros Worldwide
 
 
LinkBack Thread Tools
  #1  
Old December 17th, 2006, 10:05 AM
Sheldon
Guest
 
Posts: n/a
Default Core dump revisited

Hi,

I have a python script that uses a C extention. I keep getting a
recurring problem that causes a core dump a few lines after the C
extention return data back tp python. I tried using pbd and gdb but I
have not succeeded in understanding what went wrong and where. I post
the python script here is the error message and gdb output after
reading the core file:
...... printout from with C extention....
Completed freeing 2D arrays.
Now freeing 1D arrays
freed 1D arrays
freeing 15km arrays
Complete
Returning data back to Python!
In python: assigning tuple data to numeric arrays!
Time to do the coastal effects
Segmentation fault (core dumped)
.......

Now there next step after "Time to do the coastal effects" was never
executed. And even if I put another print statement there the same
thing happens. I am using python 2.3 and I have set my stacksize to
unlimited. When reading the core dump file with gdb I get:

gdb msgppscomp.py core.3203
GNU gdb 6.0-2mdk (Mandrake Linux)
Copyright 2003 Free Software Foundation, Inc.
GDB is free software, covered by the GNU General Public License, and
you are
welcome to change it and/or distribute copies of it under certain
conditions.
Type "show copying" to see the conditions.
There is absolutely no warranty for GDB. Type "show warranty" for
details.
This GDB was configured as
"i586-mandrake-linux-gnu"..."/data/proj_ns1/safworks/sheldon/msgppscomp.py":
not in executable format: File format not recognized

Core was generated by `python msgppscomp.py'.
Program terminated with signal 11, Segmentation fault.
#0 0x40079dfa in ?? ()
(gdb)

..................

I am not very familiar with using gdb with python and C extentions. I
would like to get some ideas about where the problem might be. I have
made sure that all the arrays in the C extention are freed before
exiting so I doubt that this might be the problem. Here is the python
script:

#!/usr/bin/env python
#!-*-coding: UTF-8 -*-

class WriteHdfFile:

def __init__(self,outfile):

self.outfile = outfile
self.COMPRESS_LVL = 6

def writeAreainfo(self):
import _pyhl
import sys
sys.float_output_precision = 2

aList = _pyhl.nodelist()
aNode = _pyhl.node(_pyhl.GROUP_ID,"/PPS_MSG_COMP")
aList.addNode(aNode)
aNode=_pyhl.node(_pyhl.DATASET_ID,"/LAT")
aNode.setArrayValue(1,lat.shape,lat,"float",-1)
aList.addNode(aNode)
aNode=_pyhl.node(_pyhl.DATASET_ID,"/LON")
aNode.setArrayValue(1,lon.shape,lon,"float",-1)
aList.addNode(aNode)
aNode=_pyhl.node(_pyhl.DATASET_ID,"/VALCONCEN")
aNode.setArrayValue(1,valconcen.shape,valconcen,"i nt",-1)
aList.addNode(aNode)
aNode=_pyhl.node(_pyhl.DATASET_ID,"/LIGHTCON")
aNode.setArrayValue(1,lightcon.shape,lightcon,"flo at",-1)
aList.addNode(aNode)
aNode=_pyhl.node(_pyhl.DATASET_ID,"/BIAS100")
aNode.setArrayValue(1,bias100.shape,bias100,"float ",-1)
aList.addNode(aNode)
aNode=_pyhl.node(_pyhl.DATASET_ID,"/BIAS75")
aNode.setArrayValue(1,bias75.shape,bias75,"float",-1)
aList.addNode(aNode)
aNode=_pyhl.node(_pyhl.DATASET_ID,"/BIAS50")
aNode.setArrayValue(1,bias50.shape,bias50,"float",-1)
aList.addNode(aNode)
aNode=_pyhl.node(_pyhl.DATASET_ID,"/BIAS25")
aNode.setArrayValue(1,bias25.shape,bias25,"float",-1)
aList.addNode(aNode)
aNode=_pyhl.node(_pyhl.DATASET_ID,"/COAST")
aNode.setArrayValue(1,coast.shape,coast,"float",-1)
aList.addNode(aNode)
aNode=_pyhl.node(_pyhl.DATASET_ID,"/STATISTICS")
aNode.setArrayValue(1,stats.shape,stats,"int",-1)
aList.addNode(aNode)
aNode=_pyhl.node(_pyhl.DATASET_ID,"/PERCENTAGE")
aNode.setArrayValue(1,percentage.shape,percentage, "float",-1)
aList.addNode(aNode)
aNode=_pyhl.node(_pyhl.DATASET_ID,"/VA_vs_BIAS")
aNode.setArrayValue(1,va_area.shape,va_area,"float ",-1)
aList.addNode(aNode)
aNode=_pyhl.node(_pyhl.DATASET_ID,"/VANUM")
aNode.setArrayValue(1,vanum.shape,vanum,"float",-1)
aList.addNode(aNode)
aNode = _pyhl.node(_pyhl.ATTRIBUTE_ID,"/NSCENES")
aNode.setScalarValue(-1,areascenes,"int",-1)
aList.addNode(aNode)
aNode = _pyhl.node(_pyhl.ATTRIBUTE_ID,"/SATELLITE")
aNode.setScalarValue(-1,satid,"string",-1)
aList.addNode(aNode)

self.status = aList.write(self.outfile,self.COMPRESS_LVL)

return self.status

#---------------------------------------------------------------------------------------------------------------------------
if __name__ == "__main__":
from Numeric import *
import sys, os, string, math, glob
import msgppsarea,msgppscoast
import shelve
date = str(200510)
#date = sys.argv[1]
#s = sys.argv[2]
cp = 'cfc'
global
valconcen,bias100,bias75,lightcon,bias50,bias25,pe rcentage,va_area,lat,lon
global stats,areascenes,satid,vanum,coast
valconcen = zeros((324,243),'i')
bias100 = zeros((324,243),'f')
bias75 = zeros((324,243),'f')
lightcon = zeros((324,243),'f')
bias50 = zeros((324,243),'f')
bias25 = zeros((324,243),'f')
coast = zeros((324,243),'f')
percentage = zeros((60,1),'f')
va_area = zeros((90,1),'f')
vanum = zeros((90,1),'f')
lat = zeros((324,243),'f')
lon = zeros((324,243),'f')
stats = zeros((60,1), 'i')
d =
shelve.open("/data/proj/safworks/sheldon/MSG_PPS_COMP/c_codes/"+"_"+date+".shelve")
msglist = d["Area.msglist"]
tilescenes = d["Area.tilescenes"]
tiles = d["Area.tiles"]
ppslist = d["Area.ppslist"]
lllist = d["Area.lllist"]
areascenes = d["Area.areascenes"]
d.close()
print "Deleting file pointer and calling C function"
RES =
msgppsarea.msgppsarea(tilescenes,tiles,msglist,pps list,lllist,int(date),areascenes)
print "In python: assigning tuple data to numeric arrays!"
for x in range(324):
bias100[x] = (RES[0])[x*243:x*243+243]
bias75[x] = (RES[1])[x*243:x*243+243]
bias50[x] = (RES[2])[x*243:x*243+243]
bias25[x] = (RES[3])[x*243:x*243+243]
lat[x] = (RES[4])[x*243:x*243+243]
lon[x] = (RES[5])[x*243:x*243+243]
lightcon[x] = (RES[6])[x*243:x*243+243]
valconcen[x] = (RES[7])[x*243:x*243+243]
for x in range(90):
va_area[x] = (RES[10])[x]
vanum[x] = (RES[11])[x]
for x in range(60):
stats[x] = (RES[8])[x]
percentage[x] = (RES[9])[x]
print "Time to do the coastal effects"
RES = msgppscoast.msgppscoast(ravel(bias100),tiles)
for x in range(324):
coast[x] = (RES[0])[x*243:x*243+243]
outfile =
"/data/proj/safworks/sheldon/MSG_PPS_COMP/results/"+'CFC'+'_'+ '_' +
date + '.h5'
print outfile
chk = WriteHdfFile(outfile)
res = chk.writeAreainfo()
print "File written!"
print "Complete!"
********************************

The C extention is very large so I will include it if anyone wants to
see it.

Thanks for any advice,

Sheldon

  #2  
Old December 18th, 2006, 11:25 AM
Gabriel Genellina
Guest
 
Posts: n/a
Default Re: Core dump revisited

On 17 dic, 07:01, "Sheldon" <shejo...@gmail.comwrote:
Quote:
I have a python script that uses a C extention. I keep getting a
recurring problem that causes a core dump a few lines after the C
extention return data back tp python. I tried using pbd and gdb but I
Most probably the error is inside the C extension, not in Python code.
Contact the author.

--
Gabriel Genellina

  #3  
Old December 18th, 2006, 03:35 PM
Nick Craig-Wood
Guest
 
Posts: n/a
Default Re: Core dump revisited

Sheldon <shejo284@gmail.comwrote:
Quote:
gdb msgppscomp.py core.3203
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.

--
Nick Craig-Wood <nick@craig-wood.com-- http://www.craig-wood.com/nick
  #4  
Old December 18th, 2006, 08:45 PM
Sheldon
Guest
 
Posts: n/a
Default Re: Core dump revisited


Nick Craig-Wood skrev:
Quote:
Sheldon <shejo284@gmail.comwrote:
Quote:
gdb msgppscomp.py core.3203
>
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.
>
--
Nick Craig-Wood <nick@craig-wood.com-- http://www.craig-wood.com/nick
Wonderful! Now I know how to used gdb with python. The are results area
posted below. Since I am new at this I could used some help in
interpreting the problem. What I do know is this: my allocation of the
array is good but when freeing the array, a problem arises. The problem
is given below but since I am new as well to C I sure could use some
help.

Program received signal SIGSEGV, Segmentation fault.
[Switching to Thread 1077321856 (LWP 32710)]
0x40297b83 in mallopt () from /lib/tls/libc.so.6
(gdb) bt
#0 0x40297b83 in mallopt () from /lib/tls/libc.so.6
#1 0x402958ba in free () from /lib/tls/libc.so.6
#2 0x405e070a in MemoryFreeOrd () at msgppsmodule_area.c:1675
#3 0x405dae0a in msgppsarea (self=0x0, args=0x4042aa7c, kwargs=0x0) at
msgppsmodule_area.c:328
#4 0x40073b16 in PyCFunction_Call () from /usr/lib/libpython2.3.so.1.0


I appreciate your help so far and could use some more.

/S

  #5  
Old December 18th, 2006, 08:55 PM
Sheldon
Guest
 
Posts: n/a
Default Re: Core dump revisited


Nick Craig-Wood skrev:
Quote:
Sheldon <shejo284@gmail.comwrote:
Quote:
gdb msgppscomp.py core.3203
>
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.
>
--
Nick Craig-Wood <nick@craig-wood.com-- http://www.craig-wood.com/nick
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[i] = malloc(sizeof 300);
work.msg_scenes[i] = malloc(sizeof 300);
}
for (i = 0; i < NumberOfTiles; i++) {
work.tiles[i] = 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[i];
strcpy(work.tile,work.tiles[i]);
/*printf("%d\n%d\n%s\n",work.Row,work.Col,work.Regi on);*/
/*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[i],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[i]);
free(work.msg_scenes[i]);
}
free(work.pps_scenes);
free(work.msg_scenes);

for (i = 0; i < NumberOfTiles; i++) {
free(work.tiles[i]);
}
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[i]);
} */
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[i] = malloc(sizeof(float)*ResCol);
work.bias10015km[i] = malloc(sizeof(float)*ResCol);
work.bias7515km[i] = malloc(sizeof(float)*ResCol);
work.bias5015km[i] = malloc(sizeof(float)*ResCol);
work.bias2515km[i] = malloc(sizeof(float)*ResCol);
work.valconcen15km[i] = 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[i]);
printf("biasdata[%d]: %f\n",i, work.biasdata[i]);
printf("Number[%d]: %f\n",i, work.number[i]);*/

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[i] = PyString_AsString(msgop);
ppsop = PyList_GetItem(work.ppsobj, i);
work.pps_scenes[i] = PyString_AsString(ppsop);
}
for (i = 0; i < NumberOfTiles; i++) {
tileop = PyList_GetItem(work.tileobj, i);
work.tiles[i] = PyString_AsString(tileop);
sceneop = PyList_GetItem(work.nscenesobj, i);
work.nscenes[i] = 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[i]);
op2 = PyFloat_FromDouble((double)work.percentage[i]);
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[i][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[i][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[i][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[i][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[i][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[i][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[i][j] == 0) work.ppsdata[i][j] = 255;
if (work.msgdata[i][j] == 0) work.msgdata[i][j] = 255;
if (work.ppsdata[i][j] == 5) work.ppsdata[i][j] = 255;
if (work.msgdata[i][j] == 5) work.msgdata[i][j] = 255;
if (work.ppsdata[i][j] == 4) work.ppsdata[i][j] = 1;
if (work.msgdata[i][j] == 4) work.msgdata[i][j] = 1;
if (work.msgdata[i][j] != 255 && work.ppsdata[i][j] != 255) {

work.valconcen[i][j]++;

if (work.ppsdata[i][j] == 2) {
work.ppsFC[i][j]++;
}
if (work.ppsdata[i][j] == 3) {
work.ppscloudy[i][j]++;
}
if (work.msgdata[i][j] == 2) {
work.msgFC[i][j]++;
}
if (work.msgdata[i][j] == 3) {
work.msgcloudy[i][j]++;
}
if (work.frac[i][j] 0 && work.frac[i][j] < 255) {
iscoast = 1;
} else {
iscoast = 0;
}

island = evaluateBIT((unsigned short int)work.qual[i][j],0);
/*iscoast = evaluateBIT((unsigned short int)work.qual[i][j],1);
coast not set in DWD files */
isnight = evaluateBIT((unsigned short int)work.qual[i][j],2);
istwilight = evaluateBIT((unsigned short int)work.qual[i][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[i][j]++;
work.stats[10]++; /*11*/
} else {
isday = 1;
work.day[i][j]++;
work.stats[0]++; /*1*/
}

/* Clear conditions */
if (work.msgdata[i][j] == clear && work.ppsdata[i][j] == clear &&
isday) {
work.stats[1]++; /*2*/
} else if (work.msgdata[i][j] == clear && work.ppsdata[i][j] == clear
&& isnight) {
work.stats[6]++; /*7*/
} else if (work.msgdata[i][j] == clear && work.ppsdata[i][j] == clear
&& istwilight) {
work.stats[11]++; /*12*/
}

/* Cloudy conditions */
if (work.msgdata[i][j] == cloudy && work.ppsdata[i][j] == cloudy &&
isday) {
work.stats[2]++; /*3*/
} else if (work.msgdata[i][j] == cloudy && work.ppsdata[i][j] ==
cloudy && isnight) {
work.stats[7]++; /*8*/
} else if (work.msgdata[i][j] == cloudy && work.ppsdata[i][j] ==
cloudy && istwilight) {
work.stats[12]++; /*13*/
}

/* Clear work.msgdata & cloudy work.ppsdata */
if (work.msgdata[i][j] == clear && work.ppsdata[i][j] == cloudy &&
isday) {
work.stats[3]++; /*4*/
} else if (work.msgdata[i][j] == clear && work.ppsdata[i][j] == cloudy
&& isnight) {
work.stats[8]++; /*9*/
} else if (work.msgdata[i][j] == clear && work.ppsdata[i][j] == cloudy
&& istwilight) {
work.stats[13]++; /*14*/
}

/* Cloudy work.msgdata & clear pps */
if (work.msgdata[i][j] == cloudy && work.ppsdata[i][j] == clear &&
isday) {
work.stats[4]++; /*5*/
} else if (work.msgdata[i][j] == cloudy && work.ppsdata[i][j] == clear
&& isnight) {
work.stats[9]++; /*10*/
} else if (work.msgdata[i][j] == cloudy && work.ppsdata[i][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[i][j] == clear && work.ppsdata[i][j] == clear &&
isday && island) {
work.stats[16]++; /*17*/
} else if (work.msgdata[i][j] == clear && work.ppsdata[i][j] == clear
&& isnight && island) {
work.stats[21]++; /*22*/
} else if (work.msgdata[i][j] == clear && work.ppsdata[i][j] == clear
&& istwilight && island) {
work.stats[26]++; /*27*/
}

/* Land cloudy */
if (work.msgdata[i][j] == cloudy && work.ppsdata[i][j] == cloudy &&
isday && island) {
work.stats[17]++; /*18*/
} else if (work.msgdata[i][j] == cloudy && work.ppsdata[i][j] ==
cloudy && isnight && island) {
work.stats[22]++; /*23*/
} else if (work.msgdata[i][j] == cloudy && work.ppsdata[i][j] ==
cloudy && istwilight && island) {
work.stats[27]++; /*28*/
}

/* work.msgdata clear pps cloudy */
if (work.msgdata[i][j] == clear && work.ppsdata[i][j] == cloudy &&
isday && island) {
work.stats[18]++; /*19*/
} else if (work.msgdata[i][j] == clear && work.ppsdata[i][j] == cloudy
&& isnight && island) {
work.stats[23]++; /*24*/
} else if (work.msgdata[i][j] == clear && work.ppsdata[i][j] == cloudy
&& istwilight && island) {
work.stats[28]++; /*29*/
}

/* work.msgdata cloudy pps clear */
if (work.msgdata[i][j] == cloudy && work.ppsdata[i][j] == clear &&
isday && island) {
work.stats[19]++; /*20*/
} else if (work.msgdata[i][j] == cloudy && work.ppsdata[i][j] == clear
&& isnight && island) {
work.stats[24]++; /*25*/
} else if (work.msgdata[i][j] == cloudy && work.ppsdata[i][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[i][j] == clear && work.ppsdata[i][j] == clear &&
isday && !island) {
work.stats[31]++; /*32*/
} else if (work.msgdata[i][j] == clear && work.ppsdata[i][j] == clear
&& isnight && !island) {
work.stats[36]++; /*37*/
} else if (work.msgdata[i][j] == clear && work.ppsdata[i][j] == clear
&& istwilight && !island) {
work.stats[41]++; /*42*/
}

/* both cloudy */
if (work.msgdata[i][j] == cloudy && work.ppsdata[i][j] == cloudy &&
isday && !island) {
work.stats[32]++; /*33*/
} else if (work.msgdata[i][j] == cloudy && work.ppsdata[i][j] ==
cloudy && isnight && !island) {
work.stats[37]++; /*38*/
} else if (work.msgdata[i][j] == cloudy && work.ppsdata[i][j] ==
cloudy && istwilight && !island) {
work.stats[42]++; /*43*/
}

/*work.msgdata clear pps cloudy */
if (work.msgdata[i][j] == clear && work.ppsdata[i][j] == cloudy &&
isday && !island) {
work.stats[33]++; /*34*/
} else if (work.msgdata[i][j] == clear && work.ppsdata[i][j] == cloudy
&& isnight && !island) {
work.stats[38]++; /*39*/
} else if (work.msgdata[i][j] == clear && work.ppsdata[i][j] == cloudy
&& istwilight && !island) {
work.stats[43]++; /*44*/
}

/* work.msgdata cloudy pps clear */
if (work.msgdata[i][j] == cloudy && work.ppsdata[i][j] == clear &&
isday && !island) {
work.stats[34]++; /*35*/
} else if (work.msgdata[i][j] == cloudy && work.ppsdata[i][j] == clear
&& isnight && !island) {
work.stats[39]++; /*40*/
} else if (work.msgdata[i][j] == cloudy && work.ppsdata[i][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[i][j] == clear && work.ppsdata[i][j] == clear &&
isday && iscoast) {
work.stats[48]++; /*49*/
} else if (work.msgdata[i][j] == clear && work.ppsdata[i][j] == clear
&& isnight && iscoast) {
work.stats[49]++; /* 50*/
} else if (work.msgdata[i][j] == clear && work.ppsdata[i][j] == clear
&& istwilight && iscoast) {
work.stats[50]++; /*51*/
}

/* Cloudy conditions */
if (work.msgdata[i][j] == cloudy && work.ppsdata[i][j] == cloudy &&
isday && iscoast) {
work.stats[51]++; /*52*/
} else if (work.msgdata[i][j] == cloudy && work.ppsdata[i][j] ==
cloudy && isnight && iscoast) {
work.stats[52]++; /*53*/
} else if (work.msgdata[i][j] == cloudy && work.ppsdata[i][j] ==
cloudy && istwilight && iscoast) {
work.stats[53]++; /*54*/
}

/* Mixed conditions */
if (work.msgdata[i][j] == cloudy && work.ppsdata[i][j] == clear &&
isday && iscoast) {
work.stats[54]++; /*55*/
} else if (work.msgdata[i][j] == cloudy && work.ppsdata[i][j] == clear
&& isnight && iscoast) {
work.stats[55]++; /*56*/
} else if (work.msgdata[i][j] == cloudy && work.ppsdata[i][j] == clear
&& istwilight && iscoast) {
work.stats[56]++; /*57*/
}

if (work.msgdata[i][j] == clear && work.ppsdata[i][j] == clear &&
isday && iscoast) {
work.stats[57]++; /*58*/
} else if (work.msgdata[i][j] == clear && work.ppsdata[i][j] == clear
&& isnight && iscoast) {
work.stats[58]++; /*59*/
} else if (work.msgdata[i][j] == clear && work.ppsdata[i][j] == clear
&& istwilight && iscoast) {
work.stats[59]++; /* 60 */
}
}
}
}
/*for (i = 0; i < work.Cat; i++) {
printf("Stats[%d]: %d\t", i, work.stats[i]);
}*/
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[i][j] = (int)hdfdata[ii];
} else if ((strcmp(satp,"p\0")) == 0) {
work.ppsdata[i][j] = (int)hdfdata[ii];
} else if ((strcmp(satp,"f\0")) == 0) {
work.frac[i][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[i][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[i][j] = work.valconcen[x][y];
work.msgFCA[i][j] = work.msgFC[x][y];
work.ppsFCA[i][j] = work.ppsFC[x][y];
work.ppscloudyA[i][j] = work.ppscloudy[x][y];
work.msgcloudyA[i][j] = work.msgcloudy[x][y];
work.dayA[i][j] = work.day[x][y];
work.twilightA[i][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[i] = (float)((float)work.stats[i] /
(float)work.stats[0] * 100.0);
if (work.percentage[i] 100.0) {
printf("Day error!\n");
goto error;
}
}
}
if (i 5 && i < 10) { /* Night only */
if (work.stats[5] 0) {
work.percentage[i] = (float)((float)work.stats[i] /
(float)work.stats[5] * 100.0);
if (work.percentage[i] 100.0) {
printf("Night error!\n");
goto error;
}
}
}
if (i 10 && i < 15) { /* Twilight only */
if (work.stats[10] 0) {
work.percentage[i] = (float)((float)work.stats[i] /
(float)work.stats[10] * 100.0);
if (work.percentage[i] 100.0) {
printf("Twilight error!\n");
goto error;
}
}
}
if (i 15 && i < 20) { /* Day land */
if (work.stats[15] 0) {
work.percentage[i] = (float)((float)work.stats[i] /
(float)work.stats[15] * 100.0);
if (work.percentage[i] 100.0) {
printf("Day land error!\n");
goto error;
}
}
}
if (i 20 && i < 25) { /* Night land */
if (work.stats[20] 0) {
work.percentage[i] = (float)((float)work.stats[i] /
(float)work.stats[20] * 100.0);
if (work.percentage[i] 100.0) {
printf("Night land error!\n");
goto error;
}
}
}
if (i 25 && i < 30) { /* Twilight land */
if (work.stats[25] 0) {
work.percentage[i] = (float)((float)work.stats[i] /
(float)work.stats[25] * 100.0);
if (work.percentage[i] 100.0) {
printf("Twilight land error!\n");
goto error;
}
}
}
if (i 30 && i < 35) { /* Water day */
if (work.stats[30] 0) {
work.percentage[i] = (float)((float)work.stats[i] /
(float)work.stats[30] * 100.0);
if (work.percentage[i] 100.0) {
printf("Water day error!\n");
goto error;
}
}
}
if (i 35 && i < 40) { /* Water night */
if (work.stats[35] != 0) {
work.percentage[i] = (float)((float)work.stats[i] /
(float)work.stats[35] * 100.0);
if (work.percentage[i] 100.0) {
printf("Water night error!\n");
goto error;
}
}
}
if (i 40 && i < 45) {
if (work.stats[40] != 0) { /* Twilight water */
work.percentage[i] = (float)((float)work.stats[i] /
(float)work.stats[40] * 100.0);
if (work.percentage[i] 100.0) {
printf("Twilight water error!\n");
goto error;
}
}
}
if (work.stats[45] 0) { /*day coast*/
if (i == 48) work.percentage[i] = (float)((float)work.stats[i] /
(float)work.stats[45] * 100.0);
if (i == 51) work.percentage[i] = (float)((float)work.stats[i] /
(float)work.stats[45] * 100.0);
if (i == 54) work.percentage[i] = (float)((float)work.stats[i] /
(float)work.stats[45] * 100.0);
if (i == 57) work.percentage[i] = (float)((float)work.stats[i] /
(float)work.stats[45] * 100.0);
if (work.percentage[i] 100.0) {
printf("Day coast error!\n");
goto error;
}
}
if (work.stats[46] 0) { /*night coast */
if (i == 49) work.percentage[i] = (float)((float)work.stats[i] /
(float)work.stats[46] * 100.0);
if (i == 52) work.percentage[i] = (float)((float)work.stats[i] /
(float)work.stats[46] * 100.0);
if (i == 55) work.percentage[i] = (float)((float)work.stats[i] /
(float)work.stats[46] * 100.0);
if (i == 58) work.percentage[i] = (float)((float)work.stats[i] /
(float)work.stats[46] * 100.0);
if (work.percentage[i] 100.0) {
printf("Night coast error!\n");
goto error;
}
}
if (work.stats[47] 0) { /*twilight coast */
if (i == 50) work.percentage[i] = (float)((float)work.stats[i] /
(float)work.stats[47] * 100.0);
if (i == 53) work.percentage[i] = (float)((float)work.stats[i] /
(float)work.stats[47] * 100.0);
if (i == 56) work.percentage[i] = (float)((float)work.stats[i] /
(float)work.stats[47] * 100.0);
if (i == 59) work.percentage[i] = (float)((float)work.stats[i] /
(float)work.stats[47] * 100.0);
if (work.percentage[i] 100.0) {
printf("Twilight coast error!\n");
goto error;
}
}
/*printf("At end of loop: percent[%d]: %f\t",i,work.percentage[i]);
*/
}
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[i][j] != 255.0) {
avg += array[i][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[i][j]);
if (array[i][j] 0) {
avg += array[i][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[i][j] = 0;
work.twilight[i][j] = 0;
work.ppsFC[i][j] = 0;
work.msgFC[i][j] = 0;
work.msgcloudy[i][j] = 0;
work.ppscloudy[i][j] = 0;
work.valconcen[i][j] = 0;
work.frac[i][j] = 0;
}
}
for (i = 0; i < work.Cat; i++) {
work.stats[i] = 0;
work.percentage[i] = 0;
}
for (i = 0; i < AreaRow; i++) {
for (j = 0; j < AreaCol; j++) {
work.dayA[i][j] = 0;
work.twilightA[i][j] = 0;
work.msgFCA[i][j] = 0;
work.ppsFCA[i][j] = 0;
work.msgcloudyA[i][j] = 0;
work.ppscloudyA[i][j] = 0;
work.lightconA[i][j] = 0.0;
work.bias100A[i][j] = 0.0;
work.bias75A[i][j] = 0.0;
work.bias50A[i][j] = 0.0;
work.bias25A[i][j] = 0.0;
work.valconcenA[i][j] = 0;
}
}
return 1;
}

static int Reset(void) {

int i,j;

printf("Reseting the tile arrays to zero\n");
for (i = 0; i < work.Row; i++) {
for (j = 0; j < work.Col; j++) {
work.day[i][j] = 0;
work.twilight[i][j] = 0;
work.ppsFC[i][j] = 0;
work.msgFC[i][j] = 0;
work.msgcloudy[i][j] = 0;
work.ppscloudy[i][j] = 0;
work.valconcen[i][j] = 0;
work.ppsdata[i][j] = 0;
work.msgdata[i][j] = 0;
work.frac[i][j] = 0;
}
}
return 1;
}

static int Bias(void) {

printf("Calculating the cloudiness per pixel\n");

int i,j;
float FC75 = 0.75;
float FC50 = 0.5;
float FC25 = 0.25;

for (i = 0; i < AreaRow; i++) {
for (j = 0; j < AreaCol; j++) {
if (work.valconcenA[i][j] 0) {
work.bias100A[i][j] =
(float)(((((float)work.msgcloudyA[i][j]+(float)work.msgFCA[i][j])/(float)work.valconcenA[i][j])*100.0)-

(float)((((float)work.ppscloudyA[i][j]+(float)work.ppsFCA[i][j])/(float)work.valconcenA[i][j])*100.0));

work.bias75A[i][j] =
(float)(((((float)work.msgcloudyA[i][j]+(float)work.msgFCA[i][j])*FC75/(float)work.valconcenA[i][j])*100.0)-

(float)((((float)work.ppscloudyA[i][j]+(float)work.ppsFCA[i][j])*FC75/(float)work.valconcenA[i][j])*100.0));

work.bias50A[i][j] =
(float)(((((float)work.msgcloudyA[i][j]+(float)work.msgFCA[i][j]*FC50)/(float)work.valconcenA[i][j])*100.0)-

(float)((((float)work.ppscloudyA[i][j]+(float)work.ppsFCA[i][j]*FC50)/(float)work.valconcenA[i][j])*100.0));

work.bias25A[i][j] =
(float)(((((float)work.msgcloudyA[i][j]+(float)work.msgFCA[i][j]*FC25)/(float)work.valconcenA[i][j])*100.0)-

(float)((((float)work.ppscloudyA[i][j]+(float)work.ppsFCA[i][j]*FC25)/(float)work.valconcenA[i][j])*100.0));
} else {
work.bias100A[i][j] = 255.0;
work.bias75A[i][j] = 255.0;
work.bias50A[i][j] = 255.0;
work.bias25A[i][j] = 255.0;
}
}
}
printf("Complete\n");
return 1;
}

static int LightCon(void) { /* Start */

printf("Calculating the light condition array\n");

int i,j;
float twi_weight = 0.3;

for (i = 0; i < AreaRow; i++) {
for (j = 0; j < AreaCol; j++) {
if(work.valconcenA[i][j] 0) {
work.lightconA[i][j] = (float)((
((float)work.dayA[i][j]+((float)work.twilightA[i][j]*twi_weight)) /
(float)work.valconcenA[i][j])*100.0);
if (work.lightconA[i][j] 100.0) {
printf("Light condition greater than 100%\n");
exit(EXIT_FAILURE);
} else if (work.lightconA[i][j] < 0) {
printf("Light condition less than 0%\n");
exit(EXIT_FAILURE);
}
} else {
work.lightconA[i][j] = 255.0;
}
}
}
printf("Complete\n");
return 1;
}

static int PrintInt(int **array, int R, int C) {

int i,j;
int min, max;

min = 0;
max = 0;
for (i = 0; i < R; i++) {
for (j = 0; j< C; j++) {
if ( i == 0 && j == 0) {
min = array[i][j];
max = array[i][j];
}
if (array[i][j] < min) {
min = array[i][j];
}
if (array[i][j] max) {
max = array[i][j];
}
}
}
printf("Arraymin: %d\tArraymax: %d\n",min,max);
return 1;
}

static int PrintFloat(float **array, int R, int C) {

int i,j;
float min,max;

for (i = 0; i < R; i++) {
for (j = 0; j< C; j++) {
if ( i == 0 && j == 0) {
min = array[i][j];
max = array[i][j];
}
if (array[i][j] < min) {
min = array[i][j];
}
if (array[i][j] max) {
max = array[i][j];
}
}
}
printf("Arraymin: %f\tArraymax: %f\n",min,max);
return 1;
}

static int Memory(void) {

int i;

/* Allocate memory */
work.day = malloc(sizeof(int*)*work.Row);
work.twilight = malloc(sizeof(int*)*work.Row);
work.msgdata = malloc(sizeof(int*)*work.Row);
work.ppsdata = malloc(sizeof(int*)*work.Row);
work.ppsFC = malloc(sizeof(int*)*work.Row);
work.msgFC = malloc(sizeof(int*)*work.Row);
work.msgcloudy = malloc(sizeof(int*)*work.Row);
work.ppscloudy = malloc(sizeof(int*)*work.Row);
work.valconcen = malloc(sizeof(int*)*work.Row);
work.frac = malloc(sizeof(int*)*work.Row);
work.qual = malloc(sizeof(unsigned short*)*work.Row);
work.stats = malloc(sizeof(float)*work.Cat);
work.percentage = malloc(sizeof(float)*work.Cat);
for (i = 0; i < work.Row; i++) {
work.day[i] = malloc(sizeof(int)*work.Col);
work.twilight[i] = malloc(sizeof(int)*work.Col);
work.msgdata[i] = malloc(sizeof(int)*work.Col);
work.ppsdata[i] = malloc(sizeof(int)*work.Col);
work.ppsFC[i] = malloc(sizeof(int)*work.Col);
work.msgFC[i] = malloc(sizeof(int)*work.Col);
work.msgcloudy[i] = malloc(sizeof(int)*work.Col);
work.ppscloudy[i] = malloc(sizeof(int)*work.Col);
work.valconcen[i] = malloc(sizeof(int)*work.Col);
work.frac[i] = malloc(sizeof(int)*work.Col);
work.qual[i] = malloc(sizeof(unsigned short)*work.Col);
}
if (work.frac==NULL) {
printf("Failed to allocate memory: frac\n");
exit(EXIT_FAILURE);
}
if (work.day==NULL) {
printf("Failed to allocate memory: day\n");
exit(EXIT_FAILURE);
}
if (work.twilight==NULL) {
printf("Failed to allocate memory: twilight\n");
exit(EXIT_FAILURE);
}
if(work.msgdata==NULL) {
printf("Failed to allocate memory: day\n");
exit(EXIT_FAILURE);
}
if (work.ppsdata==NULL) {
printf("Failed to allocating memory ppsdata\n");
exit(EXIT_FAILURE);
}
if (work.ppsFC==NULL) {
printf("Failed to allocate memory: ppsFC\n");
exit(EXIT_FAILURE);
}
if (work.msgFC==NULL) {
printf("Failed to allocate memory: msgFC\n");
exit(EXIT_FAILURE);
}
if (work.msgcloudy==NULL) {
printf("Failed to allocate memory: day\n");
exit(EXIT_FAILURE);
}
if (work.ppscloudy==NULL) {
printf("Failed to allocate memory: ppscloudy\n");
exit(EXIT_FAILURE);
}
if (work.valconcen==NULL) {
printf("Failed to allocate memory: val\n");
exit(EXIT_FAILURE);
}
if (work.qual==NULL) {
printf("Failed to allocate memory: qual\n");
exit(EXIT_FAILURE);
}
if (work.stats==NULL) {
printf("Failed to allocate memory: stats\n");
exit(EXIT_FAILURE);
}
if(work.percentage==NULL) {
printf("Failed to allocate memory: percentage\n");
exit(EXIT_FAILURE);
}
/* 2D memory */
/* Allocate memory using MACRO */
printf("Allocating 2D memory\n");
work.dayA = malloc(sizeof(int*)*AreaRow);
work.twilightA = malloc(sizeof(int*)*AreaRow);
work.ppsFCA = malloc(sizeof(int*)*AreaRow);
work.msgFCA = malloc(sizeof(int*)*AreaRow);
work.msgcloudyA = malloc(sizeof(int*)*AreaRow);
work.ppscloudyA = malloc(sizeof(int*)*AreaRow);
work.lightconA = malloc(sizeof(float*)*AreaRow);
work.bias100A = malloc(sizeof(float*)*AreaRow);
work.bias75A = malloc(sizeof(float*)*AreaRow);
work.bias50A = malloc(sizeof(float*)*AreaRow);
work.bias25A = malloc(sizeof(float*)*AreaRow);
work.valconcenA = malloc(sizeof(int*)*AreaRow);
for (i = 0; i < AreaRow; i++) {
work.dayA[i] = malloc(sizeof(int)*AreaCol);
work.twilightA[i] = malloc(sizeof(int)*AreaCol);
work.ppsFCA[i] = malloc(sizeof(int)*AreaCol);
work.msgFCA[i] = malloc(sizeof(int)*AreaCol);
work.msgcloudyA[i] = malloc(sizeof(int)*AreaCol);
work.ppscloudyA[i] = malloc(sizeof(int)*AreaCol);
work.lightconA[i] = malloc(sizeof(float)*AreaCol);
work.bias100A[i] = malloc(sizeof(float)*AreaCol);
work.bias75A[i] = malloc(sizeof(float)*AreaCol);
work.bias50A[i] = malloc(sizeof(float)*AreaCol);
work.bias25A[i] = malloc(sizeof(float)*AreaCol);
work.valconcenA[i] = malloc(sizeof(int)*AreaCol);
}
if (work.dayA==NULL) {
printf("Failed to allocate memory: dayA\n");
exit(EXIT_FAILURE);
}
if (work.twilightA==NULL) {
printf("Failed to allocate memory: twilightA\n");
exit(EXIT_FAILURE);
}
if (work.ppsFCA==NULL) {
printf("Failed to allocate memory: ppsFCA\n");
exit(EXIT_FAILURE);
}
if (work.msgFCA==NULL) {
printf("Failed to allocate memory: msgFCA\n");
exit(EXIT_FAILURE);
}
if (work.msgcloudyA==NULL) {
printf("Failed to allocate msgcloudyA\n");
exit(EXIT_FAILURE);
}
if (work.ppscloudyA==NULL) {
printf("Failed to allocate memory: ppscloudyA\n");
exit(EXIT_FAILURE);
}
if (work.lightconA==NULL) {
printf("Failed to allocate memory: lightconA\n");
exit(EXIT_FAILURE);
}
if (work.bias100A==NULL) {
printf("Failed to allocate memory: bias100A\n");
exit(EXIT_FAILURE);
}
if (work.bias75A==NULL) {
printf("Failed to allocate memory: bias75A\n");
exit(EXIT_FAILURE);
}
if (work.bias50A==NULL) {
printf("Failed to allocate memory: bias50A\n");
exit(EXIT_FAILURE);
}
if (work.bias25A==NULL) {
printf("Failed to allocate memory: bias25A\n");
exit(EXIT_FAILURE);
}
if (work.valconcenA==NULL) {
printf("Failed to allocate memory valA\n");
exit(EXIT_FAILURE);
}
printf("Allocated memory for all arrays\n");
return 1;
}

static int MemoryFreeA(void) {
int i;

for (i = 0; i < AreaCol; i++) {
free(work.dayA[i]);
free(work.twilightA[i]);
free(work.msgFCA[i]);
free(work.ppsFCA[i]);
free(work.msgcloudyA[i]);
free(work.ppscloudyA[i]);
free(work.lightconA[i]);
free(work.bias100A[i]);
free(work.bias75A[i]);
free(work.bias50A[i]);
free(work.bias25A[i]);
free(work.valconcenA[i]);
}
free(work.dayA);
free(work.twilightA);
free(work.msgFCA);
free(work.ppsFCA);
free(work.msgcloudyA);
free(work.ppscloudyA);
free(work.lightconA);
free(work.bias100A);
free(work.bias75A);
free(work.bias50A);
free(work.bias25A);
free(work.valconcenA);

return 1;
}

static int MemoryFree15km(void) {

int i;

for (i = 0; i < ResCol; i++) {
free(work.lightcon15km[i]);
free(work.bias10015km[i]);
free(work.bias7515km[i]);
free(work.bias5015km[i]);
free(work.bias2515km[i]);
free(work.valconcen15km[i]);
}
free(work.lightcon15km);
free(work.bias10015km);
free(work.bias7515km);
free(work.bias5015km);
free(work.bias2515km);
free(work.valconcen15km);

return 1;
}

static int MemoryFreeOrd(void) {

int i;

for (i = 0; i < work.Col; i++) {
free(work.day[i]);
free(work.twilight[i]);
free(work.msgcloudy[i]);
free(work.msgFC[i]);
free(work.ppsFC[i]);
free(work.ppscloudy[i]);
free(work.msgdata[i]);
free(work.ppsdata[i]);
free(work.valconcen[i]);
free(work.frac[i]);
free(work.qual[i]);
}
free(work.day);
free(work.twilight);
free(work.msgcloudy);
free(work.msgFC);
free(work.ppsFC);
free(work.ppscloudy);
free(work.msgdata);
free(work.ppsdata);
free(work.valconcen);
free(work.frac);
free(work.qual);

free(work.date);
free(work.tile);
free(work.Region);
free(work.Reprodir);
free(work.PPSfile);
free(work.MSGfile);
free(work.PHYSfile);

return 1;
}

static int MemoryFree1D(void) {

int i;

free(work.stats);
free(work.nscenes);
free(work.percentage);

return 1;
}

  #6  
Old December 18th, 2006, 09:25 PM
Sheldon
Guest
 
Posts: n/a
Default Re: Core dump revisited


Nick Craig-Wood skrev:
Quote:
Sheldon <shejo284@gmail.comwrote:
Quote:
gdb msgppscomp.py core.3203
>
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.
>
--
Nick Craig-Wood <nick@craig-wood.com-- http://www.craig-wood.com/nick
Wonderful! Now I know how to used gdb with python. The are results area
posted below. Since I am new at this I could used some help in
interpreting the problem. What I do know is this: my allocation of the
array is good but when freeing the array, a problem arises. The problem
is given below but since I am new as well to C I sure could use some
help.

Program received signal SIGSEGV, Segmentation fault.
[Switching to Thread 1077321856 (LWP 32710)]
0x40297b83 in mallopt () from /lib/tls/libc.so.6
(gdb) bt
#0 0x40297b83 in mallopt () from /lib/tls/libc.so.6
#1 0x402958ba in free () from /lib/tls/libc.so.6
#2 0x405e070a in MemoryFreeOrd () at msgppsmodule_area.c:1675
#3 0x405dae0a in msgppsarea (self=0x0, args=0x4042aa7c, kwargs=0x0) at
msgppsmodule_area.c:328
#4 0x40073b16 in PyCFunction_Call () from /usr/lib/libpython2.3.so.1.0


I appreciate your help so far and could use some more.

/S

  #7  
Old December 18th, 2006, 09:25 PM
Sheldon
Guest
 
Posts: n/a
Default Re: Core dump revisited


Sheldon skrev:
Quote:
Nick Craig-Wood skrev:
>
Quote:
Sheldon <shejo284@gmail.comwrote:
Quote:
gdb msgppscomp.py core.3203
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.

--
Nick Craig-Wood <nick@craig-wood.com-- http://www.craig-wood.com/nick
>
Wonderful! Now I know how to used gdb with python. The are results area
posted below. Since I am new at this I could used some help in
interpreting the problem. What I do know is this: my allocation of the
array is good but when freeing the array, a problem arises. The problem
is given below but since I am new as well to C I sure could use some
help.
>
Program received signal SIGSEGV, Segmentation fault.
[Switching to Thread 1077321856 (LWP 32710)]
0x40297b83 in mallopt () from /lib/tls/libc.so.6
(gdb) bt
#0 0x40297b83 in mallopt () from /lib/tls/libc.so.6
#1 0x402958ba in free () from /lib/tls/libc.so.6
#2 0x405e070a in MemoryFreeOrd () at msgppsmodule_area.c:1675
#3 0x405dae0a in msgppsarea (self=0x0, args=0x4042aa7c, kwargs=0x0) at
msgppsmodule_area.c:328
#4 0x40073b16 in PyCFunction_Call () from /usr/lib/libpython2.3.so.1.0
>
>
I appreciate your help so far and could use some more.
>
/S
Here is the 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."

#define DEBUGGING 1

#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[i] = malloc(sizeof 300);
work.msg_scenes[i] = malloc(sizeof 300);
}
for (i = 0; i < NumberOfTiles; i++) {
work.tiles[i] = 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[i];
strcpy(work.tile,work.tiles[i]);
/*printf("%d\n%d\n%s\n",work.Row,work.Col,work.Regi on);*/
/*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[i],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[i]);
free(work.msg_scenes[i]);
}
free(work.pps_scenes);
free(work.msg_scenes);

for (i = 0; i < NumberOfTiles; i++) {
free(work.tiles[i]);
}
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[i]);
} */
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[i] = malloc(sizeof(float)*ResCol);
work.bias10015km[i] = malloc(sizeof(float)*ResCol);
work.bias7515km[i] = malloc(sizeof(float)*ResCol);
work.bias5015km[i] = malloc(sizeof(float)*ResCol);
work.bias2515km[i] = malloc(sizeof(float)*ResCol);
work.valconcen15km[i] = 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[i]);
printf("biasdata[%d]: %f\n",i, work.biasdata[i]);
printf("Number[%d]: %f\n",i, work.number[i]);*/

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[i] = PyString_AsString(msgop);
ppsop = PyList_GetItem(work.ppsobj, i);
work.pps_scenes[i] = PyString_AsString(ppsop);
}
for (i = 0; i < NumberOfTiles; i++) {
tileop = PyList_GetItem(work.tileobj, i);
work.tiles[i] = PyString_AsString(tileop);
sceneop = PyList_GetItem(work.nscenesobj, i);
work.nscenes[i] = 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[i]);
op2 = PyFloat_FromDouble((double)work.percentage[i]);
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[i][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[i][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[i][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[i][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[i][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[i][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[i][j] == 0) work.ppsdata[i][j] = 255;
if (work.msgdata[i][j] == 0) work.msgdata[i][j] = 255;
if (work.ppsdata[i][j] == 5) work.ppsdata[i][j] = 255;
if (work.msgdata[i][j] == 5) work.msgdata[i][j] = 255;
if (work.ppsdata[i][j] == 4) work.ppsdata[i][j] = 1;
if (work.msgdata[i][j] == 4) work.msgdata[i][j] = 1;
if (work.msgdata[i][j] != 255 && work.ppsdata[i][j] != 255) {

work.valconcen[i][j]++;

if (work.ppsdata[i][j] == 2) {
work.ppsFC[i][j]++;
}
if (work.ppsdata[i][j] == 3) {
work.ppscloudy[i][j]++;
}
if (work.msgdata[i][j] == 2) {
work.msgFC[i][j]++;
}
if (work.msgdata[i][j] == 3) {
work.msgcloudy[i][j]++;
}
if (work.frac[i][j] 0 && work.frac[i][j] < 255) {
iscoast = 1;
} else {
iscoast = 0;
}

island = evaluateBIT((unsigned short int)work.qual[i][j],0);
/*iscoast = evaluateBIT((unsigned short int)work.qual[i][j],1);
coast not set in DWD files */
isnight = evaluateBIT((unsigned short int)work.qual[i][j],2);
istwilight = evaluateBIT((unsigned short int)work.qual[i][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[i][j]++;
work.stats[10]++; /*11*/
} else {
isday = 1;
work.day[i][j]++;
work.stats[0]++; /*1*/
}

/* Clear conditions */
if (work.msgdata[i][j] == clear && work.ppsdata[i][j] == clear &&
isday) {
work.stats[1]++; /*2*/
} else if (work.msgdata[i][j] == clear && work.ppsdata[i][j] == clear
&& isnight) {
work.stats[6]++; /*7*/
} else if (work.msgdata[i][j] == clear && work.ppsdata[i][j] == clear
&& istwilight) {
work.stats[11]++; /*12*/
}

/* Cloudy conditions */
if (work.msgdata[i][j] == cloudy && work.ppsdata[i][j] == cloudy &&
isday) {
work.stats[2]++; /*3*/
} else if (work.msgdata[i][j] == cloudy && work.ppsdata[i][j] ==
cloudy && isnight) {
work.stats[7]++; /*8*/
} else if (work.msgdata[i][j] == cloudy && work.ppsdata[i][j] ==
cloudy && istwilight) {
work.stats[12]++; /*13*/
}

/* Clear work.msgdata & cloudy work.ppsdata */
if (work.msgdata[i][j] == clear && work.ppsdata[i][j] == cloudy &&
isday) {
work.stats[3]++; /*4*/
} else if (work.msgdata[i][j] == clear && work.ppsdata[i][j] == cloudy
&& isnight) {
work.stats[8]++; /*9*/
} else if (work.msgdata[i][j] == clear && work.ppsdata[i][j] == cloudy
&& istwilight) {
work.stats[13]++; /*14*/
}

/* Cloudy work.msgdata & clear pps */
if (work.msgdata[i][j] == cloudy && work.ppsdata[i][j] == clear &&
isday) {
work.stats[4]++; /*5*/
} else if (work.msgdata[i][j] == cloudy && work.ppsdata[i][j] == clear
&& isnight) {
work.stats[9]++; /*10*/
} else if (work.msgdata[i][j] == cloudy && work.ppsdata[i][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[i][j] == clear && work.ppsdata[i][j] == clear &&
isday && island) {
work.stats[16]++; /*17*/
} else if (work.msgdata[i][j] == clear && work.ppsdata[i][j] == clear
&& isnight && island) {
work.stats[21]++; /*22*/
} else if (work.msgdata[i][j] == clear && work.ppsdata[i][j] == clear
&& istwilight && island) {
work.stats[26]++; /*27*/
}

/* Land cloudy */
if (work.msgdata[i][j] == cloudy && work.ppsdata[i][j] == cloudy &&
isday && island) {
work.stats[17]++; /*18*/
} else if (work.msgdata[i][j] == cloudy && work.ppsdata[i][j] ==
cloudy && isnight && island) {
work.stats[22]++; /*23*/
} else if (work.msgdata[i][j] == cloudy && work.ppsdata[i][j] ==
cloudy && istwilight && island) {
work.stats[27]++; /*28*/
}

/* work.msgdata clear pps cloudy */
if (work.msgdata[i][j] == clear && work.ppsdata[i][j] == cloudy &&
isday && island) {
work.stats[18]++; /*19*/
} else if (work.msgdata[i][j] == clear && work.ppsdata[i][j] == cloudy
&& isnight && island) {
work.stats[23]++; /*24*/
} else if (work.msgdata[i][j] == clear && work.ppsdata[i][j] == cloudy
&& istwilight && island) {
work.stats[28]++; /*29*/
}

/* work.msgdata cloudy pps clear */
if (work.msgdata[i][j] == cloudy && work.ppsdata[i][j] == clear &&
isday && island) {
work.stats[19]++; /*20*/
} else if (work.msgdata[i][j] == cloudy && work.ppsdata[i][j] == clear
&& isnight && island) {
work.stats[24]++; /*25*/
} else if (work.msgdata[i][j] == cloudy && work.ppsdata[i][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[i][j] == clear && work.ppsdata[i][j] == clear &&
isday && !island) {
work.stats[31]++; /*32*/
} else if (work.msgdata[i][j] == clear && work.ppsdata[i][j] == clear
&& isnight && !island) {
work.stats[36]++; /*37*/
} else if (work.msgdata[i][j] == clear && work.ppsdata[i][j] == clear
&& istwilight && !island) {
work.stats[41]++; /*42*/
}

/* both cloudy */
if (work.msgdata[i][j] == cloudy && work.ppsdata[i][j] == cloudy &&
isday && !island) {
work.stats[32]++; /*33*/
} else if (work.msgdata[i][j] == cloudy && work.ppsdata[i][j] ==
cloudy && isnight && !island) {
work.stats[37]++; /*38*/
} else if (work.msgdata[i][j] == cloudy && work.ppsdata[i][j] ==
cloudy && istwilight && !island) {
work.stats[42]++; /*43*/
}

/*work.msgdata clear pps cloudy */
if (work.msgdata[i][j] == clear && work.ppsdata[i][j] =