467,117 Members | 1,012 Online
Bytes | Developer Community
Ask Question

Home New Posts Topics Members FAQ

Post your question to a community of 467,117 developers. It's quick & easy.

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

Dec 17 '06 #1
  • viewed: 2905
Share:
14 Replies
On 17 dic, 07:01, "Sheldon" <shejo...@gmail.comwrote:
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

Dec 18 '06 #2
Sheldon <sh******@gmail.comwrote:
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 <ni**@craig-wood.com-- http://www.craig-wood.com/nick
Dec 18 '06 #3

Nick Craig-Wood skrev:
Sheldon <sh******@gmail.comwrote:
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 <ni**@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

Dec 18 '06 #4

Nick Craig-Wood skrev:
Sheldon <sh******@gmail.comwrote:
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 <ni**@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;
}

Dec 18 '06 #5

Nick Craig-Wood skrev:
Sheldon <sh******@gmail.comwrote:
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 <ni**@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

Dec 18 '06 #6

Sheldon skrev:
Nick Craig-Wood skrev:
Sheldon <sh******@gmail.comwrote:
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 <ni**@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] == 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;
}

Dec 18 '06 #7
Sheldon <sh******@gmail.comwrote:
Sheldon skrev:
Wonderful! Now I know how to used gdb with python.
Good!
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
Ambitious!
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
Typing "up" and "down" to move up and down the call stack is a good
thing. Type "l" to see a list of code at the point that things went
wrong.

You can type "print <variablename>" to see the values of variables.
gdb understands most C syntax so you can use "print
this->member[1].data" etc.

This all assumes you compiled and linked with debugging (-g flag to
compiler and linker). Turning off optimisation may make the debugger
more accurate also.

I can't tell exactly where your program crashed because of line
wrapping it your post, but is it certainly within MemoryFreeOrd().

What I would do next is run the program under gdb, wait for it to
crash then print the values of things in the MemoryFreeOrd() function.
You'll find one of the pointers has been corrupted, ie doesn't point
to valid memory or memory allocated with malloc() I expect. 0 or NULL
is an acceptable input to most free() implementations but not all.

Finding out exactly where that pointer got corrupted is harder. Maybe
it was never initialised, or maybe you initialised it with a PyObject
or maybe you have a memory scribble somewhere.

A memory scribble is the hardest bug to find because it happened
elsewhere in your program. Look at the data that was overwritten.
Maybe it is a string you can identify... You'll also want to start
reading the gdb manual on breakpoints and watchpoints at this moment!

Find memory corruptions can be tricky and time consuming.

Valgrind can help also.

Good luck!
--
Nick Craig-Wood <ni**@craig-wood.com-- http://www.craig-wood.com/nick
Dec 19 '06 #8
Sheldon wrote:
<snip>
>
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
>From my experience, segmentation faults in malloc/free is usually the
result of memory corruption.
My experience comes mostly from solaris and I realize the memory
allocator might work differently in Linux, but that would be my guess.

The cause for memory corruption can be writing outside arrays, using
deallocated memory or something like that.

Under Linux, there are tools like valgrind (free), electric fence
(free), purify (not for free last time I checked) and probably a few
more, that may help you track down memory corruption bugs. They are
often hard to find using just a debugger.

HTH

Dec 19 '06 #9

Nick Craig-Wood skrev:
Sheldon <sh******@gmail.comwrote:
Sheldon skrev:
Wonderful! Now I know how to used gdb with python.

Good!
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

Ambitious!
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

Typing "up" and "down" to move up and down the call stack is a good
thing. Type "l" to see a list of code at the point that things went
wrong.

You can type "print <variablename>" to see the values of variables.
gdb understands most C syntax so you can use "print
this->member[1].data" etc.

This all assumes you compiled and linked with debugging (-g flag to
compiler and linker). Turning off optimisation may make the debugger
more accurate also.

I can't tell exactly where your program crashed because of line
wrapping it your post, but is it certainly within MemoryFreeOrd().

What I would do next is run the program under gdb, wait for it to
crash then print the values of things in the MemoryFreeOrd() function.
You'll find one of the pointers has been corrupted, ie doesn't point
to valid memory or memory allocated with malloc() I expect. 0 or NULL
is an acceptable input to most free() implementations but not all.

Finding out exactly where that pointer got corrupted is harder. Maybe
it was never initialised, or maybe you initialised it with a PyObject
or maybe you have a memory scribble somewhere.

A memory scribble is the hardest bug to find because it happened
elsewhere in your program. Look at the data that was overwritten.
Maybe it is a string you can identify... You'll also want to start
reading the gdb manual on breakpoints and watchpoints at this moment!

Find memory corruptions can be tricky and time consuming.

Valgrind can help also.

Good luck!
--
Nick Craig-Wood <ni**@craig-wood.com-- http://www.craig-wood.com/nick

Man. You are good. This is most insight I have had from anyone. I did
initialize the arrays with PyObjects and today, after hours of
debugging and now with your insight, I think the problem lies here:

/* here a python list of strings is read into a C string array */
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*/

I am new to this and copied this code from a colleague. So, it corrupts
the pointer. How do I do this properly?

Your help is truly appreciated!
/S

Dec 19 '06 #10
"Sheldon" <sh******@gmail.comwrote:
I am new to this and copied this code from a colleague. So, it
corrupts the pointer. How do I do this properly?
Here is at least part of your problem:

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);
....
free(work.pps_scenes[i]);
free(work.msg_scenes[i]);

You initialised msg_scenes and pps_scenes with a malloc'ed block but you
then just overwrote the pointer with the result of PyString_AsString. You
don't own the memory for the string returned from PyString_AsString, so
freeing it will cause a corruption. You should copy the string data into
the malloc'ed block (with appropriate length checks).
Dec 19 '06 #11

Duncan Booth skrev:
"Sheldon" <sh******@gmail.comwrote:
I am new to this and copied this code from a colleague. So, it
corrupts the pointer. How do I do this properly?
Here is at least part of your problem:

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);
...
free(work.pps_scenes[i]);
free(work.msg_scenes[i]);

You initialised msg_scenes and pps_scenes with a malloc'ed block but you
then just overwrote the pointer with the result of PyString_AsString. You
don't own the memory for the string returned from PyString_AsString, so
freeing it will cause a corruption. You should copy the string data into
the malloc'ed block (with appropriate length checks).
Do you mean with: PyString_FromStringAndSize() and
PyString_Size(PyObject *string)

/S

Dec 19 '06 #12
"Sheldon" <sh******@gmail.comwrote:
>
Duncan Booth skrev:
>"Sheldon" <sh******@gmail.comwrote:
I am new to this and copied this code from a colleague. So, it
corrupts the pointer. How do I do this properly?
Here is at least part of your problem:

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);
...
free(work.pps_scenes[i]);
free(work.msg_scenes[i]);

You initialised msg_scenes and pps_scenes with a malloc'ed block but
you then just overwrote the pointer with the result of
PyString_AsString. You don't own the memory for the string returned
from PyString_AsString, so freeing it will cause a corruption. You
should copy the string data into the malloc'ed block (with
appropriate length checks).

Do you mean with: PyString_FromStringAndSize() and
PyString_Size(PyObject *string)
If you wish, or even just strlen() if you aren't concerned about embedded
nulls.

Dec 19 '06 #13
Sheldon <sh******@gmail.comwrote:
Man. You are good. This is most insight I have had from anyone.
:-)
I did initialize the arrays with PyObjects and today, after hours of
debugging and now with your insight, I think the problem lies here:
Good!

You need to release some python references otherwise you'll have a
memory leak and copy some strings you don't own

It is worth reading the definitions for those functions in the Python
API docs

http://docs.python.org/api/api.html

In particular read about reference counts.

With the Python C API you have to know at all time (by reading the
doc) whether you own the reference or not, and whether you own the
memory or not.
/* here a python list of strings is read into a C string array */
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);
work.msg_scenes[i] = strdup(PyString_AsString(msgop));
Py_DECREF(msgop);
ppsop = PyList_GetItem(work.ppsobj, i);
work.pps_scenes[i] = PyString_AsString(ppsop);
work.pps_scenes[i] = strdup(PyString_AsString(ppsop));
Py_DECREF(ppsop);
}
for (i = 0; i < NumberOfTiles; i++) {
tileop = PyList_GetItem(work.tileobj, i);
work.tiles[i] = PyString_AsString(tileop);
work.tiles[i] = strdup(PyString_AsString(tileop));
Py_DECREF(tileop);
sceneop = PyList_GetItem(work.nscenesobj, i);
work.nscenes[i] = PyInt_AsLong(sceneop);
Py_DECREF(sceneop);
}
return 1;
} /*end readPythonObject*/
You free() the strings later which is fine.

The above ignores errors which PyList_GetItem may return and strdup()
returning 0, but it should get you out of trouble hopefully.

....

I've written lots of quite similar code, here is a snippet. Note the
comments about who owns the reference, and the error checking. Also
note xstrdup() which is a strdup() which blows up if no memory is
available.

[snip]
PyObject *item = PySequence_GetItem(value, i); /* real ref */
if (item == 0)
{
fprintf(stderr, "Failed to read '%s[%d]' attribute\n", name, i);
goto err;
}
item_cstr = PyString_AsString(item); /* borrowed */
if (item_cstr == 0)
{
fprintf(stderr, "Failed to read '%s[%d]' as string\n", name, i);
goto err;
}
label[i] = xstrdup(item_cstr);
Py_DECREF(item);
item = 0;
[snip]
err:;
PyErr_Print();
out:;
if (value)
Py_DECREF(value);
if (item)
Py_DECREF(item);
return rc;
--
Nick Craig-Wood <ni**@craig-wood.com-- http://www.craig-wood.com/nick
Dec 19 '06 #14

Nick Craig-Wood skrev:
Sheldon <sh******@gmail.comwrote:
Man. You are good. This is most insight I have had from anyone.

:-)
I did initialize the arrays with PyObjects and today, after hours of
debugging and now with your insight, I think the problem lies here:

Good!

You need to release some python references otherwise you'll have a
memory leak and copy some strings you don't own

It is worth reading the definitions for those functions in the Python
API docs

http://docs.python.org/api/api.html

In particular read about reference counts.

With the Python C API you have to know at all time (by reading the
doc) whether you own the reference or not, and whether you own the
memory or not.
/* here a python list of strings is read into a C string array */
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);

work.msg_scenes[i] = strdup(PyString_AsString(msgop));
Py_DECREF(msgop);
ppsop = PyList_GetItem(work.ppsobj, i);
work.pps_scenes[i] = PyString_AsString(ppsop);

work.pps_scenes[i] = strdup(PyString_AsString(ppsop));
Py_DECREF(ppsop);
}
for (i = 0; i < NumberOfTiles; i++) {
tileop = PyList_GetItem(work.tileobj, i);
work.tiles[i] = PyString_AsString(tileop);

work.tiles[i] = strdup(PyString_AsString(tileop));
Py_DECREF(tileop);
sceneop = PyList_GetItem(work.nscenesobj, i);
work.nscenes[i] = PyInt_AsLong(sceneop);

Py_DECREF(sceneop);
}
return 1;
} /*end readPythonObject*/

You free() the strings later which is fine.

The above ignores errors which PyList_GetItem may return and strdup()
returning 0, but it should get you out of trouble hopefully.

...

I've written lots of quite similar code, here is a snippet. Note the
comments about who owns the reference, and the error checking. Also
note xstrdup() which is a strdup() which blows up if no memory is
available.

[snip]
PyObject *item = PySequence_GetItem(value, i); /* real ref */
if (item == 0)
{
fprintf(stderr, "Failed to read '%s[%d]' attribute\n", name, i);
goto err;
}
item_cstr = PyString_AsString(item); /* borrowed */
if (item_cstr == 0)
{
fprintf(stderr, "Failed to read '%s[%d]' as string\n", name, i);
goto err;
}
label[i] = xstrdup(item_cstr);
Py_DECREF(item);
item = 0;
[snip]
err:;
PyErr_Print();
out:;
if (value)
Py_DECREF(value);
if (item)
Py_DECREF(item);
return rc;
--
Nick Craig-Wood <ni**@craig-wood.com-- http://www.craig-wood.com/nick
Thanks Nick! Man, I really appreciate this. I did dereference before
but the program crashed and I did understand why so I kept on learning.
Now I have your snip. Much obliged.
I will be adding Numpy to my new PC in Jan and there I can pass arrays
instead of list. Still, I would like to save your email address i case
I have some mre questions later - if that ok with you?
I will make the changes and try to run the code.

/S

Dec 19 '06 #15

This discussion thread is closed

Replies have been disabled for this discussion.

Similar topics

3 posts views Thread by Nick Craig-Wood | last post: by
1 post views Thread by Martin | last post: by
10 posts views Thread by ken | last post: by
1 post views Thread by invincible | last post: by
10 posts views Thread by John Liu | last post: by
3 posts views Thread by John Liu | last post: by
4 posts views Thread by madhusudan.hv@gmail.com | last post: by
10 posts views Thread by wong_powah@yahoo.ca | last post: by
5 posts views Thread by johnericaturnbull@yahoo.com | last post: by
By using this site, you agree to our Privacy Policy and Terms of Use.