473,803 Members | 3,943 Online
Bytes | Software Development & Data Engineering Community
+ Post

Home Posts Topics Members FAQ

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,o utfile):

self.outfile = outfile
self.COMPRESS_L VL = 6

def writeAreainfo(s elf):
import _pyhl
import sys
sys.float_outpu t_precision = 2

aList = _pyhl.nodelist( )
aNode = _pyhl.node(_pyh l.GROUP_ID,"/PPS_MSG_COMP")
aList.addNode(a Node)
aNode=_pyhl.nod e(_pyhl.DATASET _ID,"/LAT")
aNode.setArrayV alue(1,lat.shap e,lat,"float",-1)
aList.addNode(a Node)
aNode=_pyhl.nod e(_pyhl.DATASET _ID,"/LON")
aNode.setArrayV alue(1,lon.shap e,lon,"float",-1)
aList.addNode(a Node)
aNode=_pyhl.nod e(_pyhl.DATASET _ID,"/VALCONCEN")
aNode.setArrayV alue(1,valconce n.shape,valconc en,"int",-1)
aList.addNode(a Node)
aNode=_pyhl.nod e(_pyhl.DATASET _ID,"/LIGHTCON")
aNode.setArrayV alue(1,lightcon .shape,lightcon ,"float",-1)
aList.addNode(a Node)
aNode=_pyhl.nod e(_pyhl.DATASET _ID,"/BIAS100")
aNode.setArrayV alue(1,bias100. shape,bias100," float",-1)
aList.addNode(a Node)
aNode=_pyhl.nod e(_pyhl.DATASET _ID,"/BIAS75")
aNode.setArrayV alue(1,bias75.s hape,bias75,"fl oat",-1)
aList.addNode(a Node)
aNode=_pyhl.nod e(_pyhl.DATASET _ID,"/BIAS50")
aNode.setArrayV alue(1,bias50.s hape,bias50,"fl oat",-1)
aList.addNode(a Node)
aNode=_pyhl.nod e(_pyhl.DATASET _ID,"/BIAS25")
aNode.setArrayV alue(1,bias25.s hape,bias25,"fl oat",-1)
aList.addNode(a Node)
aNode=_pyhl.nod e(_pyhl.DATASET _ID,"/COAST")
aNode.setArrayV alue(1,coast.sh ape,coast,"floa t",-1)
aList.addNode(a Node)
aNode=_pyhl.nod e(_pyhl.DATASET _ID,"/STATISTICS")
aNode.setArrayV alue(1,stats.sh ape,stats,"int" ,-1)
aList.addNode(a Node)
aNode=_pyhl.nod e(_pyhl.DATASET _ID,"/PERCENTAGE")
aNode.setArrayV alue(1,percenta ge.shape,percen tage,"float",-1)
aList.addNode(a Node)
aNode=_pyhl.nod e(_pyhl.DATASET _ID,"/VA_vs_BIAS")
aNode.setArrayV alue(1,va_area. shape,va_area," float",-1)
aList.addNode(a Node)
aNode=_pyhl.nod e(_pyhl.DATASET _ID,"/VANUM")
aNode.setArrayV alue(1,vanum.sh ape,vanum,"floa t",-1)
aList.addNode(a Node)
aNode = _pyhl.node(_pyh l.ATTRIBUTE_ID, "/NSCENES")
aNode.setScalar Value(-1,areascenes,"i nt",-1)
aList.addNode(a Node)
aNode = _pyhl.node(_pyh l.ATTRIBUTE_ID, "/SATELLITE")
aNode.setScalar Value(-1,satid,"string ",-1)
aList.addNode(a Node)

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

return self.status

#---------------------------------------------------------------------------------------------------------------------------
if __name__ == "__main__":
from Numeric import *
import sys, os, string, math, glob
import msgppsarea,msgp pscoast
import shelve
date = str(200510)
#date = sys.argv[1]
#s = sys.argv[2]
cp = 'cfc'
global
valconcen,bias1 00,bias75,light con,bias50,bias 25,percentage,v a_area,lat,lon
global stats,areascene s,satid,vanum,c oast
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+".sh elve")
msglist = d["Area.msgli st"]
tilescenes = d["Area.tilescene s"]
tiles = d["Area.tiles "]
ppslist = d["Area.ppsli st"]
lllist = d["Area.lllis t"]
areascenes = d["Area.areascene s"]
d.close()
print "Deleting file pointer and calling C function"
RES =
msgppsarea.msgp psarea(tilescen es,tiles,msglis t,ppslist,lllis t,int(date),are ascenes)
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.msg ppscoast(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(ou tfile)
res = chk.writeAreain fo()
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
14 3364
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_ar ea.c:1675
#3 0x405dae0a in msgppsarea (self=0x0, args=0x4042aa7c , kwargs=0x0) at
msgppsmodule_ar ea.c:328
#4 0x40073b16 in PyCFunction_Cal l () 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_W ITH_IMPORT_ARRA Y
#include "acpg_arrayobje ct_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_vhlhd f.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 createPythonObj ect(void);
static int readPythonObjec t(void);
PyMODINIT_FUNC initmsgppsarea( void); /* msgppsarea will be the name of
the SO file */
/* Copied from Sarah */
static int arrShortFromHdf Node(HL_NodeLis t *nodelist, char *nodename);
static int arrIntFromHdfNo de(HL_NodeList *nodelist, char *nodename,
char *satp);
/* *************** ******** */
static int readHDFPPS(void );
static int readHDFMSG(void );
static int readHDFQual(voi d);
static int readHDFPHYS(voi d);
static int evaluateBIT(uns igned short int value, int bitnumber);
static int Statistics(void );
static int FillArea(void);
static int Percent(void);
static int InterpAvgF(floa t** 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(vo id);
static int MemoryFreeOrd(v oid);
static int MemoryFree15km( void);
static int MemoryFreeA(voi d);

/* 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_68n2 9w","0","1215", "0","1215",
"06","ibla_68n0 0e","0","1215", "1215","243 0",
"07","ibla_68n2 9e","0","1215", "2430","364 5",
"0B","ibla_57n2 0w","1215","243 0","0","1215 ",
"0C","ibla_57n0 0e","1215","243 0","1215","2430 ",
"0D","ibla_57n2 0e","1215","243 0","2430","3645 ",
"13","ibla_46n1 6w","2430","364 5","0","1215 ",
"14","ibla_46n0 0e","2430","364 5","1215","2430 ",
"15","ibla_46n1 6e","2430","364 5","2430","3645 ",
"1D","ibla_35n1 3w","3645","486 0","0","1215 ",
"1E","ibla_35n0 0e","3645","486 0","1215","2430 ",
"1F","ibla_35n1 3e","3645","486 0","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=NU LL;

static PyObject* msgppsarea(PyOb ject* self, PyObject* args, PyObject
*kwargs) { /* Start main function */

mtrace ();
mcheck(NULL);

static char *kwlist[] =
{"nscenesobj"," tileobj","msgob j","ppsobj","da teint","sumscen es",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_ParseTu pleAndKeywords( args,kwargs,"OO OOii:nothing",
kwlist,
&work.nscenesob j,
&work.tileob j,
&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.da te,"%d",work.da tein);
printf("Date: %s\n", work.date);
/* Allocate memory to python objects */
work.pps_scenes = malloc(sizeof *work.pps_scene s*work.sumscene s);
work.msg_scenes = malloc(sizeof *work.msg_scene s*work.sumscene s);
if (work.pps_scene s == NULL) {
printf("fail to allocate memory!");
exit(EXIT_FAILU RE);
}
if (work.msg_scene s == NULL) {
printf("fail to allocate memory!");
exit(EXIT_FAILU RE);
}
work.tiles = malloc(sizeof*w ork.tiles*Numbe rOfTiles);
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(i nt)*NumberOfTil es);
if (work.pps_scene s==NULL) {
printf("Failed to allocate memory for pps and msg scenes\n");
goto fail;
}
if (work.msg_scene s==NULL) {
printf("Failed to allocate memory for pps and msg scenes\n");
goto fail;
}
if (work.tiles==NU LL) {
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 (!readPythonObj ect()) {
printf("readPyt honObject_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.til e,work.tiles[i]);
/*printf("%d\n%d \n%s\n",work.Ro w,work.Col,work .Region);*/
/*printf("\nDoin g tile: %s\n", work.tile);*/
for (xp = 0; xp < NumberOfTiles; xp++) {
if (strncmp(work.t ile,tileinfo[xp][0],2) == 0) {
strcpy(work.Reg ion,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.PHY Sfile,PHYS);
strcat(work.PHY Sfile,work.Regi on);
strcat(work.PHY Sfile,".hdf");
printf("PHYSfil e: %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("\nProce ssing %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.MSG file,work.msg_s cenes[q]);
/*printf("%s\n", work.MSGfile);*/
strcpy(work.PPS file,work.pps_s cenes[q]);
/*printf("%s\n", work.PPSfile);*/
/*printf("Copyin g 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_s cenes[i]);
free(work.msg_s cenes[i]);
}
free(work.pps_s cenes);
free(work.msg_s cenes);

for (i = 0; i < NumberOfTiles; i++) {
free(work.tiles[i]);
}
free(work.tiles );
free(work.nscen es);

if(!Percent()) {
printf("Failed calculate percentages\n") ;
goto fail;
}
/* for (i = 0; i < work.Cat; i++) {
printf("percent age[%d]: %f\n",i,work.pe rcentage[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 (!(MemoryFreeOr d())) {
printf("Allocat ed free memory for ordinary arrays\n");
goto fail;
}
work.lightcon15 km = malloc(sizeof(f loat*)*ResRow);
work.bias10015k m = malloc(sizeof(f loat*)*ResRow);
work.bias7515km = malloc(sizeof(f loat*)*ResRow);
work.bias5015km = malloc(sizeof(f loat*)*ResRow);
work.bias2515km = malloc(sizeof(f loat*)*ResRow);
work.valconcen1 5km = malloc(sizeof(i nt*)*ResRow);
for (i = 0; i < ResRow; i++) {
work.lightcon15 km[i] = malloc(sizeof(f loat)*ResCol);
work.bias10015k m[i] = malloc(sizeof(f loat)*ResCol);
work.bias7515km[i] = malloc(sizeof(f loat)*ResCol);
work.bias5015km[i] = malloc(sizeof(f loat)*ResCol);
work.bias2515km[i] = malloc(sizeof(f loat)*ResCol);
work.valconcen1 5km[i] = malloc(sizeof(i nt)*ResCol);
}
if ((work.lightcon 15km)==NULL) {
printf("Failed to allocating lightcon memory\n");
exit(EXIT_FAILU RE);
}
if ((work.bias1001 5km)==NULL) {
printf("Failed to allocating bias100 memory\n");
exit(EXIT_FAILU RE);
}
if ((work.bias7515 km)==NULL) {
printf("Failed to allocating bias75 memory\n");
exit(EXIT_FAILU RE);
}
if ((work.bias5015 km)==NULL) {
printf("Failed to allocating bias50 memory\n");
exit(EXIT_FAILU RE);
}
if ((work.bias2515 km)==NULL) {
printf("Failed to allocating bias25 memory\n");
exit(EXIT_FAILU RE);
}
if ((work.valconce n15km)==NULL) {
printf("Failed to allocating val memory\n");
exit(EXIT_FAILU RE);
}
/* Interpolating data array */
if(!InterpAvgF( work.lightconA, work.lightcon15 km)) {
printf("Failed interpolation lightcon\n");
goto fail;
}
if (!InterpAvgI(wo rk.valconcenA,w ork.valconcen15 km)) {
printf("Failed interpolation val\n");
goto fail;
}
if(!InterpAvgF( work.bias100A,w ork.bias10015km )) {
printf("Failed interpolation b100\n");
goto fail;
}
if (!InterpAvgF(wo rk.bias75A,work .bias7515km)) {
printf("Failed interpolation b75\n");
goto fail;
}
if(!InterpAvgF( work.bias50A,wo rk.bias5015km)) {
printf("Failed interpolation b50\n");
goto fail;
}
if (!InterpAvgF(wo rk.bias25A,work .bias2515km)) {
printf("Failed interpolation b25\n");
goto fail;
}
/* Free Area arrays */
if ((MemoryFreeA() )) {
printf("Allocat ed free memory for area arrays\n");
goto fail;
}
/*printf("VA[%d]: %f\n",i, work.va_area[i]);
printf("biasdat a[%d]: %f\n",i, work.biasdata[i]);
printf("Number[%d]: %f\n",i, work.number[i]);*/

printf("Creatin g python objects\n");

Rlightcon=PyTup le_New(ResRow*R esCol);
Rbias100=PyTupl e_New(ResRow*Re sCol);
Rbias75=PyTuple _New(ResRow*Res Col);
Rbias50=PyTuple _New(ResRow*Res Col);
Rbias25=PyTuple _New(ResRow*Res Col);
Rvalcon=PyTuple _New(ResRow*Res Col);
Rstats=PyTuple_ New(work.Cat);
Rpercentage=PyT uple_New(work.C at);

if (!(createPython Object())) {
printf("Failed to create python arrays objects\n");
goto fail;
}
printf("Returni ng data back to Python!\n");
muntrace ();
return Py_BuildValue(" OOOOOOOO",Rbias 100,
Rbias75,
Rbias50,
Rbias25,
Rlightcon,
Rvalcon,
Rstats,
Rpercentage);

fail:
exit(EXIT_FAILU RE);
}
/* Functions for the C-wrapper */
/* This structure is needed for the python-C-interface */
static struct PyMethodDef msgppsareaMetho ds[] = {
{"msgppsarea ", (PyCFunction)ms gppsarea,
METH_VARARGS|ME TH_KEYWORDS,"HE LP for msgppsarea\n"},
{NULL,NULL,0,NU LL}
};

PyMODINIT_FUNC initmsgppsarea( void) {
(void) Py_InitModule(" msgppsarea",msg ppsareaMethods) ;
import_array(); /* access to Numeric PyArray functions */
}

/*Read an list of strings from a python object*/
static int readPythonObjec t(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_AsStri ng(msgop);
ppsop = PyList_GetItem( work.ppsobj, i);
work.pps_scenes[i] = PyString_AsStri ng(ppsop);
}
for (i = 0; i < NumberOfTiles; i++) {
tileop = PyList_GetItem( work.tileobj, i);
work.tiles[i] = PyString_AsStri ng(tileop);
sceneop = PyList_GetItem( work.nscenesobj , i);
work.nscenes[i] = PyInt_AsLong(sc eneop);
}
return 1;
} /*end readPythonObjec t*/

static int createPythonObj ect(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.stat s[i]);
op2 = PyFloat_FromDou ble((double)wor k.percentage[i]);
if (PyTuple_SetIte m(Rstats,i,op1) !=0) {
printf("Error in creating python object Rstats\n");
exit(EXIT_FAILU RE);
}
if (PyTuple_SetIte m(Rpercentage,i ,op2) !=0) {
printf("Error in creating python object Rpercentage\n") ;
exit(EXIT_FAILU RE);
}
}
printf("Complet ed 1D arrays\n");
k=0;
for (i = 0; i < ResRow; i++) {
for (j = 0; j < ResCol; j++) {
ligop = PyFloat_FromDou ble((double)wor k.lightcon15km[i][j]);
if (PyTuple_SetIte m(Rlightcon,k,l igop) !=0) {
printf("Error in creating python light object\n");
exit(EXIT_FAILU RE);
}
valop = PyInt_FromLong( (long)work.valc oncen15km[i][j]);
if (PyTuple_SetIte m(Rvalcon,k,val op) !=0) {
printf("Error in creating python val object\n");
exit(EXIT_FAILU RE);
}
b100op = PyFloat_FromDou ble((double)wor k.bias10015km[i][j]);
if (PyTuple_SetIte m(Rbias100,k,b1 00op) !=0) {
printf("Error in creating python bias100 object\n");
exit(EXIT_FAILU RE);
}
b75op = PyFloat_FromDou ble((double)wor k.bias7515km[i][j]);
if (PyTuple_SetIte m(Rbias75,k,b75 op) !=0) {
printf("Error in creating python bias75 object\n");
exit(EXIT_FAILU RE);
}
b50op = PyFloat_FromDou ble((double)wor k.bias5015km[i][j]);
if (PyTuple_SetIte m(Rbias50,k,b50 op) !=0) {
printf("Error in creating python bias50 object\n");
exit(EXIT_FAILU RE);
}
b25op = PyFloat_FromDou ble((double)wor k.bias2515km[i][j]);
if (PyTuple_SetIte m(Rbias25,k,b25 op) !=0) {
printf("Error in creating python bias25 object\n");
exit(EXIT_FAILU RE);
}
k++;
}
}
printf("Complet ed 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 (!(MemoryFree15 km())) {
printf("Couldn' t free memory for 15km arrays\n");
goto fail;
}
printf("Complet e\n");
return 1;
fail:
exit(EXIT_FAILU RE);
} /* end createPythonObj ect */

static int readHDFPHYS(voi d) {

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("Finishe d setting inithdf\n");
}

if (!(nodelist=rea dHL_NodeList(wo rk.PHYSfile))) {
printf("Error reading nodelist! from files: %s\n", work.PHYSfile);
goto fail;
}
/*printf("Checki ng dataset's existens\n"); */
if (!selectAllNode s(nodelist)) {
printf("Fail to select items\n");
goto fail;
}
if(!fetchMarked Nodes(nodelist) ) {
printf("Failed to fetch data\n");
goto fail;
}
/* Reading cloudmask data into arrays */
if (!arrIntFromHdf Node(nodelist,f ield,"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_FAILU RE);
}

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("Finishe d setting inithdf\n");
}

if (!(nodelist=rea dHL_NodeList(wo rk.PPSfile))) {
printf("Error reading nodelist! from files: %s\n", work.PPSfile);
goto fail;
}
/*printf("Checki ng dataset's existens\n"); */
if (!selectAllNode s(nodelist)) {
printf("Fail to select items\n");
goto fail;
}
if(!fetchMarked Nodes(nodelist) ) {
printf("Failed to fetch data\n");
goto fail;
}
/* Reading cloudmask data into arrays */
if (!arrIntFromHdf Node(nodelist,f ield,"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_FAILU RE);
}

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("Finishe d setting inithdf\n");
}

if (!(nodelist=rea dHL_NodeList(wo rk.MSGfile))) {
printf("Error reading nodelist! from files: %s\n", work.MSGfile);
goto fail;
}
/*printf("Checki ng dataset's existens\n"); */
if (!selectAllNode s(nodelist)) {
printf("Fail to select items\n");
goto fail;
}
if(!fetchMarked Nodes(nodelist) ) {
printf("Failed to fetch data\n");
goto fail;
}
/* Reading cloudmask data into arrays */
if (!arrIntFromHdf Node(nodelist,f ield,"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_FAILU RE);
}

static int readHDFQual(voi d) {

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("Finishe d setting inithdf\n");
}

if (!(nodelist=rea dHL_NodeList(wo rk.PPSfile))) {
printf("Error reading nodelist from files: %s\n", work.PPSfile);
goto fail;
}
if (!selectAllNode s(nodelist)) {
printf("Fail to select items\n");
goto fail;
}
if(!fetchMarked Nodes(nodelist) ) {
printf("Failed to fetch data\n");
goto fail;
}
if (!(arrShortFrom HdfNode(nodelis t,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_FAILU RE);
}

/* Statistics function */
static int Statistics(void ) { /* Start */
int clear = 1;
int cloudy = 3;
int i, j;
int island, iscoast, isnight, istwilight, isday;
/*printf("Calcul ating 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((un signed short int)work.qual[i][j],0);
/*iscoast = evaluateBIT((un signed short int)work.qual[i][j],1);
coast not set in DWD files */
isnight = evaluateBIT((un signed short int)work.qual[i][j],2);
istwilight = evaluateBIT((un signed 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(uns igned 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 arrIntFromHdfNo de(HL_NodeList *nodelist, char *nodename,
char *satp) {

HL_Node *node;
unsigned char *hdfdata=NULL;
int ii,i,j;

/*Read from node*/
if(!(node=getNo de(nodelist,nod ename))){
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(i nt)*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_FAILU RE);
} /*end arrIntFromHdfNo de*/

/*Read an array of short from node.
Give it as output in form of an array of short */
static int arrShortFromHdf Node(HL_NodeLis t *nodelist, char *nodename) {

HL_Node *node;
unsigned short *hdfdata=NULL;
int ii;
int i,j;

/*Read from node*/
if(!(node=getNo de(nodelist,nod ename))) {
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(u nsigned 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!=NUL L)
free(hdfdata);

return 1;
fail:
exit(EXIT_FAILU RE);
} /*end arrShortFromHdf Node*/

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("Calcula ting 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.sta ts[0] * 100.0);
if (work.percentag e[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.sta ts[5] * 100.0);
if (work.percentag e[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.sta ts[10] * 100.0);
if (work.percentag e[i] 100.0) {
printf("Twiligh t 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.sta ts[15] * 100.0);
if (work.percentag e[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.sta ts[20] * 100.0);
if (work.percentag e[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.sta ts[25] * 100.0);
if (work.percentag e[i] 100.0) {
printf("Twiligh t 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.sta ts[30] * 100.0);
if (work.percentag e[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.sta ts[35] * 100.0);
if (work.percentag e[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.sta ts[40] * 100.0);
if (work.percentag e[i] 100.0) {
printf("Twiligh t 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.sta ts[45] * 100.0);
if (i == 51) work.percentage[i] = (float)((float) work.stats[i] /
(float)work.sta ts[45] * 100.0);
if (i == 54) work.percentage[i] = (float)((float) work.stats[i] /
(float)work.sta ts[45] * 100.0);
if (i == 57) work.percentage[i] = (float)((float) work.stats[i] /
(float)work.sta ts[45] * 100.0);
if (work.percentag e[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.sta ts[46] * 100.0);
if (i == 52) work.percentage[i] = (float)((float) work.stats[i] /
(float)work.sta ts[46] * 100.0);
if (i == 55) work.percentage[i] = (float)((float) work.stats[i] /
(float)work.sta ts[46] * 100.0);
if (i == 58) work.percentage[i] = (float)((float) work.stats[i] /
(float)work.sta ts[46] * 100.0);
if (work.percentag e[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.sta ts[47] * 100.0);
if (i == 53) work.percentage[i] = (float)((float) work.stats[i] /
(float)work.sta ts[47] * 100.0);
if (i == 56) work.percentage[i] = (float)((float) work.stats[i] /
(float)work.sta ts[47] * 100.0);
if (i == 59) work.percentage[i] = (float)((float) work.stats[i] /
(float)work.sta ts[47] * 100.0);
if (work.percentag e[i] 100.0) {
printf("Twiligh t coast error!\n");
goto error;
}
}
/*printf("At end of loop: percent[%d]: %f\t",i,work.pe rcentage[i]);
*/
}
return 1;
error:
exit(EXIT_FAILU RE);
}

/* 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(floa t** 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)av g / (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("Initial izing 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("Resetin g 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("Calcula ting 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.valconcen A[i][j] 0) {
work.bias100A[i][j] =
(float)(((((flo at)work.msgclou dyA[i][j]+(float)work.ms gFCA[i][j])/(float)work.val concenA[i][j])*100.0)-

(float)((((floa t)work.ppscloud yA[i][j]+(float)work.pp sFCA[i][j])/(float)work.val concenA[i][j])*100.0));

work.bias75A[i][j] =
(float)(((((flo at)work.msgclou dyA[i][j]+(float)work.ms gFCA[i][j])*FC75/(float)work.val concenA[i][j])*100.0)-

(float)((((floa t)work.ppscloud yA[i][j]+(float)work.pp sFCA[i][j])*FC75/(float)work.val concenA[i][j])*100.0));

work.bias50A[i][j] =
(float)(((((flo at)work.msgclou dyA[i][j]+(float)work.ms gFCA[i][j]*FC50)/(float)work.val concenA[i][j])*100.0)-

(float)((((floa t)work.ppscloud yA[i][j]+(float)work.pp sFCA[i][j]*FC50)/(float)work.val concenA[i][j])*100.0));

work.bias25A[i][j] =
(float)(((((flo at)work.msgclou dyA[i][j]+(float)work.ms gFCA[i][j]*FC25)/(float)work.val concenA[i][j])*100.0)-

(float)((((floa t)work.ppscloud yA[i][j]+(float)work.pp sFCA[i][j]*FC25)/(float)work.val concenA[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("Complet e\n");
return 1;
}

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

printf("Calcula ting 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.valconc enA[i][j] 0) {
work.lightconA[i][j] = (float)((
((float)work.da yA[i][j]+((float)work.t wilightA[i][j]*twi_weight)) /
(float)work.val concenA[i][j])*100.0);
if (work.lightconA[i][j] 100.0) {
printf("Light condition greater than 100%\n");
exit(EXIT_FAILU RE);
} else if (work.lightconA[i][j] < 0) {
printf("Light condition less than 0%\n");
exit(EXIT_FAILU RE);
}
} else {
work.lightconA[i][j] = 255.0;
}
}
}
printf("Complet e\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("Arraymi n: %d\tArraymax: %d\n",min,max);
return 1;
}

static int PrintFloat(floa t **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("Arraymi n: %f\tArraymax: %f\n",min,max);
return 1;
}

static int Memory(void) {

int i;

/* Allocate memory */
work.day = malloc(sizeof(i nt*)*work.Row);
work.twilight = malloc(sizeof(i nt*)*work.Row);
work.msgdata = malloc(sizeof(i nt*)*work.Row);
work.ppsdata = malloc(sizeof(i nt*)*work.Row);
work.ppsFC = malloc(sizeof(i nt*)*work.Row);
work.msgFC = malloc(sizeof(i nt*)*work.Row);
work.msgcloudy = malloc(sizeof(i nt*)*work.Row);
work.ppscloudy = malloc(sizeof(i nt*)*work.Row);
work.valconcen = malloc(sizeof(i nt*)*work.Row);
work.frac = malloc(sizeof(i nt*)*work.Row);
work.qual = malloc(sizeof(u nsigned short*)*work.Ro w);
work.stats = malloc(sizeof(f loat)*work.Cat) ;
work.percentage = malloc(sizeof(f loat)*work.Cat) ;
for (i = 0; i < work.Row; i++) {
work.day[i] = malloc(sizeof(i nt)*work.Col);
work.twilight[i] = malloc(sizeof(i nt)*work.Col);
work.msgdata[i] = malloc(sizeof(i nt)*work.Col);
work.ppsdata[i] = malloc(sizeof(i nt)*work.Col);
work.ppsFC[i] = malloc(sizeof(i nt)*work.Col);
work.msgFC[i] = malloc(sizeof(i nt)*work.Col);
work.msgcloudy[i] = malloc(sizeof(i nt)*work.Col);
work.ppscloudy[i] = malloc(sizeof(i nt)*work.Col);
work.valconcen[i] = malloc(sizeof(i nt)*work.Col);
work.frac[i] = malloc(sizeof(i nt)*work.Col);
work.qual[i] = malloc(sizeof(u nsigned short)*work.Col );
}
if (work.frac==NUL L) {
printf("Failed to allocate memory: frac\n");
exit(EXIT_FAILU RE);
}
if (work.day==NULL ) {
printf("Failed to allocate memory: day\n");
exit(EXIT_FAILU RE);
}
if (work.twilight= =NULL) {
printf("Failed to allocate memory: twilight\n");
exit(EXIT_FAILU RE);
}
if(work.msgdata ==NULL) {
printf("Failed to allocate memory: day\n");
exit(EXIT_FAILU RE);
}
if (work.ppsdata== NULL) {
printf("Failed to allocating memory ppsdata\n");
exit(EXIT_FAILU RE);
}
if (work.ppsFC==NU LL) {
printf("Failed to allocate memory: ppsFC\n");
exit(EXIT_FAILU RE);
}
if (work.msgFC==NU LL) {
printf("Failed to allocate memory: msgFC\n");
exit(EXIT_FAILU RE);
}
if (work.msgcloudy ==NULL) {
printf("Failed to allocate memory: day\n");
exit(EXIT_FAILU RE);
}
if (work.ppscloudy ==NULL) {
printf("Failed to allocate memory: ppscloudy\n");
exit(EXIT_FAILU RE);
}
if (work.valconcen ==NULL) {
printf("Failed to allocate memory: val\n");
exit(EXIT_FAILU RE);
}
if (work.qual==NUL L) {
printf("Failed to allocate memory: qual\n");
exit(EXIT_FAILU RE);
}
if (work.stats==NU LL) {
printf("Failed to allocate memory: stats\n");
exit(EXIT_FAILU RE);
}
if(work.percent age==NULL) {
printf("Failed to allocate memory: percentage\n");
exit(EXIT_FAILU RE);
}
/* 2D memory */
/* Allocate memory using MACRO */
printf("Allocat ing 2D memory\n");
work.dayA = malloc(sizeof(i nt*)*AreaRow);
work.twilightA = malloc(sizeof(i nt*)*AreaRow);
work.ppsFCA = malloc(sizeof(i nt*)*AreaRow);
work.msgFCA = malloc(sizeof(i nt*)*AreaRow);
work.msgcloudyA = malloc(sizeof(i nt*)*AreaRow);
work.ppscloudyA = malloc(sizeof(i nt*)*AreaRow);
work.lightconA = malloc(sizeof(f loat*)*AreaRow) ;
work.bias100A = malloc(sizeof(f loat*)*AreaRow) ;
work.bias75A = malloc(sizeof(f loat*)*AreaRow) ;
work.bias50A = malloc(sizeof(f loat*)*AreaRow) ;
work.bias25A = malloc(sizeof(f loat*)*AreaRow) ;
work.valconcenA = malloc(sizeof(i nt*)*AreaRow);
for (i = 0; i < AreaRow; i++) {
work.dayA[i] = malloc(sizeof(i nt)*AreaCol);
work.twilightA[i] = malloc(sizeof(i nt)*AreaCol);
work.ppsFCA[i] = malloc(sizeof(i nt)*AreaCol);
work.msgFCA[i] = malloc(sizeof(i nt)*AreaCol);
work.msgcloudyA[i] = malloc(sizeof(i nt)*AreaCol);
work.ppscloudyA[i] = malloc(sizeof(i nt)*AreaCol);
work.lightconA[i] = malloc(sizeof(f loat)*AreaCol);
work.bias100A[i] = malloc(sizeof(f loat)*AreaCol);
work.bias75A[i] = malloc(sizeof(f loat)*AreaCol);
work.bias50A[i] = malloc(sizeof(f loat)*AreaCol);
work.bias25A[i] = malloc(sizeof(f loat)*AreaCol);
work.valconcenA[i] = malloc(sizeof(i nt)*AreaCol);
}
if (work.dayA==NUL L) {
printf("Failed to allocate memory: dayA\n");
exit(EXIT_FAILU RE);
}
if (work.twilightA ==NULL) {
printf("Failed to allocate memory: twilightA\n");
exit(EXIT_FAILU RE);
}
if (work.ppsFCA==N ULL) {
printf("Failed to allocate memory: ppsFCA\n");
exit(EXIT_FAILU RE);
}
if (work.msgFCA==N ULL) {
printf("Failed to allocate memory: msgFCA\n");
exit(EXIT_FAILU RE);
}
if (work.msgcloudy A==NULL) {
printf("Failed to allocate msgcloudyA\n");
exit(EXIT_FAILU RE);
}
if (work.ppscloudy A==NULL) {
printf("Failed to allocate memory: ppscloudyA\n");
exit(EXIT_FAILU RE);
}
if (work.lightconA ==NULL) {
printf("Failed to allocate memory: lightconA\n");
exit(EXIT_FAILU RE);
}
if (work.bias100A= =NULL) {
printf("Failed to allocate memory: bias100A\n");
exit(EXIT_FAILU RE);
}
if (work.bias75A== NULL) {
printf("Failed to allocate memory: bias75A\n");
exit(EXIT_FAILU RE);
}
if (work.bias50A== NULL) {
printf("Failed to allocate memory: bias50A\n");
exit(EXIT_FAILU RE);
}
if (work.bias25A== NULL) {
printf("Failed to allocate memory: bias25A\n");
exit(EXIT_FAILU RE);
}
if (work.valconcen A==NULL) {
printf("Failed to allocate memory valA\n");
exit(EXIT_FAILU RE);
}
printf("Allocat ed memory for all arrays\n");
return 1;
}

static int MemoryFreeA(voi d) {
int i;

for (i = 0; i < AreaCol; i++) {
free(work.dayA[i]);
free(work.twili ghtA[i]);
free(work.msgFC A[i]);
free(work.ppsFC A[i]);
free(work.msgcl oudyA[i]);
free(work.ppscl oudyA[i]);
free(work.light conA[i]);
free(work.bias1 00A[i]);
free(work.bias7 5A[i]);
free(work.bias5 0A[i]);
free(work.bias2 5A[i]);
free(work.valco ncenA[i]);
}
free(work.dayA) ;
free(work.twili ghtA);
free(work.msgFC A);
free(work.ppsFC A);
free(work.msgcl oudyA);
free(work.ppscl oudyA);
free(work.light conA);
free(work.bias1 00A);
free(work.bias7 5A);
free(work.bias5 0A);
free(work.bias2 5A);
free(work.valco ncenA);

return 1;
}

static int MemoryFree15km( void) {

int i;

for (i = 0; i < ResCol; i++) {
free(work.light con15km[i]);
free(work.bias1 0015km[i]);
free(work.bias7 515km[i]);
free(work.bias5 015km[i]);
free(work.bias2 515km[i]);
free(work.valco ncen15km[i]);
}
free(work.light con15km);
free(work.bias1 0015km);
free(work.bias7 515km);
free(work.bias5 015km);
free(work.bias2 515km);
free(work.valco ncen15km);

return 1;
}

static int MemoryFreeOrd(v oid) {

int i;

for (i = 0; i < work.Col; i++) {
free(work.day[i]);
free(work.twili ght[i]);
free(work.msgcl oudy[i]);
free(work.msgFC[i]);
free(work.ppsFC[i]);
free(work.ppscl oudy[i]);
free(work.msgda ta[i]);
free(work.ppsda ta[i]);
free(work.valco ncen[i]);
free(work.frac[i]);
free(work.qual[i]);
}
free(work.day);
free(work.twili ght);
free(work.msgcl oudy);
free(work.msgFC );
free(work.ppsFC );
free(work.ppscl oudy);
free(work.msgda ta);
free(work.ppsda ta);
free(work.valco ncen);
free(work.frac) ;
free(work.qual) ;

free(work.date) ;
free(work.tile) ;
free(work.Regio n);
free(work.Repro dir);
free(work.PPSfi le);
free(work.MSGfi le);
free(work.PHYSf ile);

return 1;
}

static int MemoryFree1D(vo id) {

int i;

free(work.stats );
free(work.nscen es);
free(work.perce ntage);

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_ar ea.c:1675
#3 0x405dae0a in msgppsarea (self=0x0, args=0x4042aa7c , kwargs=0x0) at
msgppsmodule_ar ea.c:328
#4 0x40073b16 in PyCFunction_Cal l () 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_ar ea.c:1675
#3 0x405dae0a in msgppsarea (self=0x0, args=0x4042aa7c , kwargs=0x0) at
msgppsmodule_ar ea.c:328
#4 0x40073b16 in PyCFunction_Cal l () 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_W ITH_IMPORT_ARRA Y
#include "acpg_arrayobje ct_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_vhlhd f.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 createPythonObj ect(void);
static int readPythonObjec t(void);
PyMODINIT_FUNC initmsgppsarea( void); /* msgppsarea will be the name of
the SO file */
/* Copied from Sarah */
static int arrShortFromHdf Node(HL_NodeLis t *nodelist, char *nodename);
static int arrIntFromHdfNo de(HL_NodeList *nodelist, char *nodename,
char *satp);
/* *************** ******** */
static int readHDFPPS(void );
static int readHDFMSG(void );
static int readHDFQual(voi d);
static int readHDFPHYS(voi d);
static int evaluateBIT(uns igned short int value, int bitnumber);
static int Statistics(void );
static int FillArea(void);
static int Percent(void);
static int InterpAvgF(floa t** 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(vo id);
static int MemoryFreeOrd(v oid);
static int MemoryFree15km( void);
static int MemoryFreeA(voi d);

/* 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_68n2 9w","0","1215", "0","1215",
"06","ibla_68n0 0e","0","1215", "1215","243 0",
"07","ibla_68n2 9e","0","1215", "2430","364 5",
"0B","ibla_57n2 0w","1215","243 0","0","1215 ",
"0C","ibla_57n0 0e","1215","243 0","1215","2430 ",
"0D","ibla_57n2 0e","1215","243 0","2430","3645 ",
"13","ibla_46n1 6w","2430","364 5","0","1215 ",
"14","ibla_46n0 0e","2430","364 5","1215","2430 ",
"15","ibla_46n1 6e","2430","364 5","2430","3645 ",
"1D","ibla_35n1 3w","3645","486 0","0","1215 ",
"1E","ibla_35n0 0e","3645","486 0","1215","2430 ",
"1F","ibla_35n1 3e","3645","486 0","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=NU LL;

static PyObject* msgppsarea(PyOb ject* self, PyObject* args, PyObject
*kwargs) { /* Start main function */

mtrace ();
mcheck(NULL);

static char *kwlist[] =
{"nscenesobj"," tileobj","msgob j","ppsobj","da teint","sumscen es",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_ParseTu pleAndKeywords( args,kwargs,"OO OOii:nothing",
kwlist,
&work.nscenesob j,
&work.tileob j,
&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.da te,"%d",work.da tein);
printf("Date: %s\n", work.date);
/* Allocate memory to python objects */
work.pps_scenes = malloc(sizeof *work.pps_scene s*work.sumscene s);
work.msg_scenes = malloc(sizeof *work.msg_scene s*work.sumscene s);
if (work.pps_scene s == NULL) {
printf("fail to allocate memory!");
exit(EXIT_FAILU RE);
}
if (work.msg_scene s == NULL) {
printf("fail to allocate memory!");
exit(EXIT_FAILU RE);
}
work.tiles = malloc(sizeof*w ork.tiles*Numbe rOfTiles);
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(i nt)*NumberOfTil es);
if (work.pps_scene s==NULL) {
printf("Failed to allocate memory for pps and msg scenes\n");
goto fail;
}
if (work.msg_scene s==NULL) {
printf("Failed to allocate memory for pps and msg scenes\n");
goto fail;
}
if (work.tiles==NU LL) {
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 (!readPythonObj ect()) {
printf("readPyt honObject_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.til e,work.tiles[i]);
/*printf("%d\n%d \n%s\n",work.Ro w,work.Col,work .Region);*/
/*printf("\nDoin g tile: %s\n", work.tile);*/
for (xp = 0; xp < NumberOfTiles; xp++) {
if (strncmp(work.t ile,tileinfo[xp][0],2) == 0) {
strcpy(work.Reg ion,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.PHY Sfile,PHYS);
strcat(work.PHY Sfile,work.Regi on);
strcat(work.PHY Sfile,".hdf");
printf("PHYSfil e: %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("\nProce ssing %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.MSG file,work.msg_s cenes[q]);
/*printf("%s\n", work.MSGfile);*/
strcpy(work.PPS file,work.pps_s cenes[q]);
/*printf("%s\n", work.PPSfile);*/
/*printf("Copyin g 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_s cenes[i]);
free(work.msg_s cenes[i]);
}
free(work.pps_s cenes);
free(work.msg_s cenes);

for (i = 0; i < NumberOfTiles; i++) {
free(work.tiles[i]);
}
free(work.tiles );
free(work.nscen es);

if(!Percent()) {
printf("Failed calculate percentages\n") ;
goto fail;
}
/* for (i = 0; i < work.Cat; i++) {
printf("percent age[%d]: %f\n",i,work.pe rcentage[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 (!(MemoryFreeOr d())) {
printf("Allocat ed free memory for ordinary arrays\n");
goto fail;
}
work.lightcon15 km = malloc(sizeof(f loat*)*ResRow);
work.bias10015k m = malloc(sizeof(f loat*)*ResRow);
work.bias7515km = malloc(sizeof(f loat*)*ResRow);
work.bias5015km = malloc(sizeof(f loat*)*ResRow);
work.bias2515km = malloc(sizeof(f loat*)*ResRow);
work.valconcen1 5km = malloc(sizeof(i nt*)*ResRow);
for (i = 0; i < ResRow; i++) {
work.lightcon15 km[i] = malloc(sizeof(f loat)*ResCol);
work.bias10015k m[i] = malloc(sizeof(f loat)*ResCol);
work.bias7515km[i] = malloc(sizeof(f loat)*ResCol);
work.bias5015km[i] = malloc(sizeof(f loat)*ResCol);
work.bias2515km[i] = malloc(sizeof(f loat)*ResCol);
work.valconcen1 5km[i] = malloc(sizeof(i nt)*ResCol);
}
if ((work.lightcon 15km)==NULL) {
printf("Failed to allocating lightcon memory\n");
exit(EXIT_FAILU RE);
}
if ((work.bias1001 5km)==NULL) {
printf("Failed to allocating bias100 memory\n");
exit(EXIT_FAILU RE);
}
if ((work.bias7515 km)==NULL) {
printf("Failed to allocating bias75 memory\n");
exit(EXIT_FAILU RE);
}
if ((work.bias5015 km)==NULL) {
printf("Failed to allocating bias50 memory\n");
exit(EXIT_FAILU RE);
}
if ((work.bias2515 km)==NULL) {
printf("Failed to allocating bias25 memory\n");
exit(EXIT_FAILU RE);
}
if ((work.valconce n15km)==NULL) {
printf("Failed to allocating val memory\n");
exit(EXIT_FAILU RE);
}
/* Interpolating data array */
if(!InterpAvgF( work.lightconA, work.lightcon15 km)) {
printf("Failed interpolation lightcon\n");
goto fail;
}
if (!InterpAvgI(wo rk.valconcenA,w ork.valconcen15 km)) {
printf("Failed interpolation val\n");
goto fail;
}
if(!InterpAvgF( work.bias100A,w ork.bias10015km )) {
printf("Failed interpolation b100\n");
goto fail;
}
if (!InterpAvgF(wo rk.bias75A,work .bias7515km)) {
printf("Failed interpolation b75\n");
goto fail;
}
if(!InterpAvgF( work.bias50A,wo rk.bias5015km)) {
printf("Failed interpolation b50\n");
goto fail;
}
if (!InterpAvgF(wo rk.bias25A,work .bias2515km)) {
printf("Failed interpolation b25\n");
goto fail;
}
/* Free Area arrays */
if ((MemoryFreeA() )) {
printf("Allocat ed free memory for area arrays\n");
goto fail;
}
/*printf("VA[%d]: %f\n",i, work.va_area[i]);
printf("biasdat a[%d]: %f\n",i, work.biasdata[i]);
printf("Number[%d]: %f\n",i, work.number[i]);*/

printf("Creatin g python objects\n");

Rlightcon=PyTup le_New(ResRow*R esCol);
Rbias100=PyTupl e_New(ResRow*Re sCol);
Rbias75=PyTuple _New(ResRow*Res Col);
Rbias50=PyTuple _New(ResRow*Res Col);
Rbias25=PyTuple _New(ResRow*Res Col);
Rvalcon=PyTuple _New(ResRow*Res Col);
Rstats=PyTuple_ New(work.Cat);
Rpercentage=PyT uple_New(work.C at);

if (!(createPython Object())) {
printf("Failed to create python arrays objects\n");
goto fail;
}
printf("Returni ng data back to Python!\n");
muntrace ();
return Py_BuildValue(" OOOOOOOO",Rbias 100,
Rbias75,
Rbias50,
Rbias25,
Rlightcon,
Rvalcon,
Rstats,
Rpercentage);

fail:
exit(EXIT_FAILU RE);
}
/* Functions for the C-wrapper */
/* This structure is needed for the python-C-interface */
static struct PyMethodDef msgppsareaMetho ds[] = {
{"msgppsarea ", (PyCFunction)ms gppsarea,
METH_VARARGS|ME TH_KEYWORDS,"HE LP for msgppsarea\n"},
{NULL,NULL,0,NU LL}
};

PyMODINIT_FUNC initmsgppsarea( void) {
(void) Py_InitModule(" msgppsarea",msg ppsareaMethods) ;
import_array(); /* access to Numeric PyArray functions */
}

/*Read an list of strings from a python object*/
static int readPythonObjec t(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_AsStri ng(msgop);
ppsop = PyList_GetItem( work.ppsobj, i);
work.pps_scenes[i] = PyString_AsStri ng(ppsop);
}
for (i = 0; i < NumberOfTiles; i++) {
tileop = PyList_GetItem( work.tileobj, i);
work.tiles[i] = PyString_AsStri ng(tileop);
sceneop = PyList_GetItem( work.nscenesobj , i);
work.nscenes[i] = PyInt_AsLong(sc eneop);
}
return 1;
} /*end readPythonObjec t*/

static int createPythonObj ect(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.stat s[i]);
op2 = PyFloat_FromDou ble((double)wor k.percentage[i]);
if (PyTuple_SetIte m(Rstats,i,op1) !=0) {
printf("Error in creating python object Rstats\n");
exit(EXIT_FAILU RE);
}
if (PyTuple_SetIte m(Rpercentage,i ,op2) !=0) {
printf("Error in creating python object Rpercentage\n") ;
exit(EXIT_FAILU RE);
}
}
printf("Complet ed 1D arrays\n");
k=0;
for (i = 0; i < ResRow; i++) {
for (j = 0; j < ResCol; j++) {
ligop = PyFloat_FromDou ble((double)wor k.lightcon15km[i][j]);
if (PyTuple_SetIte m(Rlightcon,k,l igop) !=0) {
printf("Error in creating python light object\n");
exit(EXIT_FAILU RE);
}
valop = PyInt_FromLong( (long)work.valc oncen15km[i][j]);
if (PyTuple_SetIte m(Rvalcon,k,val op) !=0) {
printf("Error in creating python val object\n");
exit(EXIT_FAILU RE);
}
b100op = PyFloat_FromDou ble((double)wor k.bias10015km[i][j]);
if (PyTuple_SetIte m(Rbias100,k,b1 00op) !=0) {
printf("Error in creating python bias100 object\n");
exit(EXIT_FAILU RE);
}
b75op = PyFloat_FromDou ble((double)wor k.bias7515km[i][j]);
if (PyTuple_SetIte m(Rbias75,k,b75 op) !=0) {
printf("Error in creating python bias75 object\n");
exit(EXIT_FAILU RE);
}
b50op = PyFloat_FromDou ble((double)wor k.bias5015km[i][j]);
if (PyTuple_SetIte m(Rbias50,k,b50 op) !=0) {
printf("Error in creating python bias50 object\n");
exit(EXIT_FAILU RE);
}
b25op = PyFloat_FromDou ble((double)wor k.bias2515km[i][j]);
if (PyTuple_SetIte m(Rbias25,k,b25 op) !=0) {
printf("Error in creating python bias25 object\n");
exit(EXIT_FAILU RE);
}
k++;
}
}
printf("Complet ed 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 (!(MemoryFree15 km())) {
printf("Couldn' t free memory for 15km arrays\n");
goto fail;
}
printf("Complet e\n");
return 1;
fail:
exit(EXIT_FAILU RE);
} /* end createPythonObj ect */

static int readHDFPHYS(voi d) {

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("Finishe d setting inithdf\n");
}

if (!(nodelist=rea dHL_NodeList(wo rk.PHYSfile))) {
printf("Error reading nodelist! from files: %s\n", work.PHYSfile);
goto fail;
}
/*printf("Checki ng dataset's existens\n"); */
if (!selectAllNode s(nodelist)) {
printf("Fail to select items\n");
goto fail;
}
if(!fetchMarked Nodes(nodelist) ) {
printf("Failed to fetch data\n");
goto fail;
}
/* Reading cloudmask data into arrays */
if (!arrIntFromHdf Node(nodelist,f ield,"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_FAILU RE);
}

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("Finishe d setting inithdf\n");
}

if (!(nodelist=rea dHL_NodeList(wo rk.PPSfile))) {
printf("Error reading nodelist! from files: %s\n", work.PPSfile);
goto fail;
}
/*printf("Checki ng dataset's existens\n"); */
if (!selectAllNode s(nodelist)) {
printf("Fail to select items\n");
goto fail;
}
if(!fetchMarked Nodes(nodelist) ) {
printf("Failed to fetch data\n");
goto fail;
}
/* Reading cloudmask data into arrays */
if (!arrIntFromHdf Node(nodelist,f ield,"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_FAILU RE);
}

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("Finishe d setting inithdf\n");
}

if (!(nodelist=rea dHL_NodeList(wo rk.MSGfile))) {
printf("Error reading nodelist! from files: %s\n", work.MSGfile);
goto fail;
}
/*printf("Checki ng dataset's existens\n"); */
if (!selectAllNode s(nodelist)) {
printf("Fail to select items\n");
goto fail;
}
if(!fetchMarked Nodes(nodelist) ) {
printf("Failed to fetch data\n");
goto fail;
}
/* Reading cloudmask data into arrays */
if (!arrIntFromHdf Node(nodelist,f ield,"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_FAILU RE);
}

static int readHDFQual(voi d) {

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("Finishe d setting inithdf\n");
}

if (!(nodelist=rea dHL_NodeList(wo rk.PPSfile))) {
printf("Error reading nodelist from files: %s\n", work.PPSfile);
goto fail;
}
if (!selectAllNode s(nodelist)) {
printf("Fail to select items\n");
goto fail;
}
if(!fetchMarked Nodes(nodelist) ) {
printf("Failed to fetch data\n");
goto fail;
}
if (!(arrShortFrom HdfNode(nodelis t,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_FAILU RE);
}

/* Statistics function */
static int Statistics(void ) { /* Start */
int clear = 1;
int cloudy = 3;
int i, j;
int island, iscoast, isnight, istwilight, isday;
/*printf("Calcul ating 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((un signed short int)work.qual[i][j],0);
/*iscoast = evaluateBIT((un signed short int)work.qual[i][j],1);
coast not set in DWD files */
isnight = evaluateBIT((un signed short int)work.qual[i][j],2);
istwilight = evaluateBIT((un signed 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(uns igned 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 arrIntFromHdfNo de(HL_NodeList *nodelist, char *nodename,
char *satp) {

HL_Node *node;
unsigned char *hdfdata=NULL;
int ii,i,j;

/*Read from node*/
if(!(node=getNo de(nodelist,nod ename))){
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(i nt)*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_FAILU RE);
} /*end arrIntFromHdfNo de*/

/*Read an array of short from node.
Give it as output in form of an array of short */
static int arrShortFromHdf Node(HL_NodeLis t *nodelist, char *nodename) {

HL_Node *node;
unsigned short *hdfdata=NULL;
int ii;
int i,j;

/*Read from node*/
if(!(node=getNo de(nodelist,nod ename))) {
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(u nsigned 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!=NUL L)
free(hdfdata);

return 1;
fail:
exit(EXIT_FAILU RE);
} /*end arrShortFromHdf Node*/

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("Calcula ting 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.sta ts[0] * 100.0);
if (work.percentag e[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.sta ts[5] * 100.0);
if (work.percentag e[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.sta ts[10] * 100.0);
if (work.percentag e[i] 100.0) {
printf("Twiligh t 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.sta ts[15] * 100.0);
if (work.percentag e[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.sta ts[20] * 100.0);
if (work.percentag e[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.sta ts[25] * 100.0);
if (work.percentag e[i] 100.0) {
printf("Twiligh t 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.sta ts[30] * 100.0);
if (work.percentag e[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.sta ts[35] * 100.0);
if (work.percentag e[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.sta ts[40] * 100.0);
if (work.percentag e[i] 100.0) {
printf("Twiligh t 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.sta ts[45] * 100.0);
if (i == 51) work.percentage[i] = (float)((float) work.stats[i] /
(float)work.sta ts[45] * 100.0);
if (i == 54) work.percentage[i] = (float)((float) work.stats[i] /
(float)work.sta ts[45] * 100.0);
if (i == 57) work.percentage[i] = (float)((float) work.stats[i] /
(float)work.sta ts[45] * 100.0);
if (work.percentag e[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.sta ts[46] * 100.0);
if (i == 52) work.percentage[i] = (float)((float) work.stats[i] /
(float)work.sta ts[46] * 100.0);
if (i == 55) work.percentage[i] = (float)((float) work.stats[i] /
(float)work.sta ts[46] * 100.0);
if (i == 58) work.percentage[i] = (float)((float) work.stats[i] /
(float)work.sta ts[46] * 100.0);
if (work.percentag e[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.sta ts[47] * 100.0);
if (i == 53) work.percentage[i] = (float)((float) work.stats[i] /
(float)work.sta ts[47] * 100.0);
if (i == 56) work.percentage[i] = (float)((float) work.stats[i] /
(float)work.sta ts[47] * 100.0);
if (i == 59) work.percentage[i] = (float)((float) work.stats[i] /
(float)work.sta ts[47] * 100.0);
if (work.percentag e[i] 100.0) {
printf("Twiligh t coast error!\n");
goto error;
}
}
/*printf("At end of loop: percent[%d]: %f\t",i,work.pe rcentage[i]);
*/
}
return 1;
error:
exit(EXIT_FAILU RE);
}

/* 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(floa t** 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)av g / (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("Initial izing 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("Resetin g 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("Calcula ting 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.valconcen A[i][j] 0) {
work.bias100A[i][j] =
(float)(((((flo at)work.msgclou dyA[i][j]+(float)work.ms gFCA[i][j])/(float)work.val concenA[i][j])*100.0)-

(float)((((floa t)work.ppscloud yA[i][j]+(float)work.pp sFCA[i][j])/(float)work.val concenA[i][j])*100.0));

work.bias75A[i][j] =
(float)(((((flo at)work.msgclou dyA[i][j]+(float)work.ms gFCA[i][j])*FC75/(float)work.val concenA[i][j])*100.0)-

(float)((((floa t)work.ppscloud yA[i][j]+(float)work.pp sFCA[i][j])*FC75/(float)work.val concenA[i][j])*100.0));

work.bias50A[i][j] =
(float)(((((flo at)work.msgclou dyA[i][j]+(float)work.ms gFCA[i][j]*FC50)/(float)work.val concenA[i][j])*100.0)-

(float)((((floa t)work.ppscloud yA[i][j]+(float)work.pp sFCA[i][j]*FC50)/(float)work.val concenA[i][j])*100.0));

work.bias25A[i][j] =
(float)(((((flo at)work.msgclou dyA[i][j]+(float)work.ms gFCA[i][j]*FC25)/(float)work.val concenA[i][j])*100.0)-

(float)((((floa t)work.ppscloud yA[i][j]+(float)work.pp sFCA[i][j]*FC25)/(float)work.val concenA[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("Complet e\n");
return 1;
}

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

printf("Calcula ting 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.valconc enA[i][j] 0) {
work.lightconA[i][j] = (float)((
((float)work.da yA[i][j]+((float)work.t wilightA[i][j]*twi_weight)) /
(float)work.val concenA[i][j])*100.0);
if (work.lightconA[i][j] 100.0) {
printf("Light condition greater than 100%\n");
exit(EXIT_FAILU RE);
} else if (work.lightconA[i][j] < 0) {
printf("Light condition less than 0%\n");
exit(EXIT_FAILU RE);
}
} else {
work.lightconA[i][j] = 255.0;
}
}
}
printf("Complet e\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("Arraymi n: %d\tArraymax: %d\n",min,max);
return 1;
}

static int PrintFloat(floa t **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("Arraymi n: %f\tArraymax: %f\n",min,max);
return 1;
}

static int Memory(void) {

int i;

/* Allocate memory */
work.day = malloc(sizeof(i nt*)*work.Row);
work.twilight = malloc(sizeof(i nt*)*work.Row);
work.msgdata = malloc(sizeof(i nt*)*work.Row);
work.ppsdata = malloc(sizeof(i nt*)*work.Row);
work.ppsFC = malloc(sizeof(i nt*)*work.Row);
work.msgFC = malloc(sizeof(i nt*)*work.Row);
work.msgcloudy = malloc(sizeof(i nt*)*work.Row);
work.ppscloudy = malloc(sizeof(i nt*)*work.Row);
work.valconcen = malloc(sizeof(i nt*)*work.Row);
work.frac = malloc(sizeof(i nt*)*work.Row);
work.qual = malloc(sizeof(u nsigned short*)*work.Ro w);
work.stats = malloc(sizeof(f loat)*work.Cat) ;
work.percentage = malloc(sizeof(f loat)*work.Cat) ;
for (i = 0; i < work.Row; i++) {
work.day[i] = malloc(sizeof(i nt)*work.Col);
work.twilight[i] = malloc(sizeof(i nt)*work.Col);
work.msgdata[i] = malloc(sizeof(i nt)*work.Col);
work.ppsdata[i] = malloc(sizeof(i nt)*work.Col);
work.ppsFC[i] = malloc(sizeof(i nt)*work.Col);
work.msgFC[i] = malloc(sizeof(i nt)*work.Col);
work.msgcloudy[i] = malloc(sizeof(i nt)*work.Col);
work.ppscloudy[i] = malloc(sizeof(i nt)*work.Col);
work.valconcen[i] = malloc(sizeof(i nt)*work.Col);
work.frac[i] = malloc(sizeof(i nt)*work.Col);
work.qual[i] = malloc(sizeof(u nsigned short)*work.Col );
}
if (work.frac==NUL L) {
printf("Failed to allocate memory: frac\n");
exit(EXIT_FAILU RE);
}
if (work.day==NULL ) {
printf("Failed to allocate memory: day\n");
exit(EXIT_FAILU RE);
}
if (work.twilight= =NULL) {
printf("Failed to allocate memory: twilight\n");
exit(EXIT_FAILU RE);
}
if(work.msgdata ==NULL) {
printf("Failed to allocate memory: day\n");
exit(EXIT_FAILU RE);
}
if (work.ppsdata== NULL) {
printf("Failed to allocating memory ppsdata\n");
exit(EXIT_FAILU RE);
}
if (work.ppsFC==NU LL) {
printf("Failed to allocate memory: ppsFC\n");
exit(EXIT_FAILU RE);
}
if (work.msgFC==NU LL) {
printf("Failed to allocate memory: msgFC\n");
exit(EXIT_FAILU RE);
}
if (work.msgcloudy ==NULL) {
printf("Failed to allocate memory: day\n");
exit(EXIT_FAILU RE);
}
if (work.ppscloudy ==NULL) {
printf("Failed to allocate memory: ppscloudy\n");
exit(EXIT_FAILU RE);
}
if (work.valconcen ==NULL) {
printf("Failed to allocate memory: val\n");
exit(EXIT_FAILU RE);
}
if (work.qual==NUL L) {
printf("Failed to allocate memory: qual\n");
exit(EXIT_FAILU RE);
}
if (work.stats==NU LL) {
printf("Failed to allocate memory: stats\n");
exit(EXIT_FAILU RE);
}
if(work.percent age==NULL) {
printf("Failed to allocate memory: percentage\n");
exit(EXIT_FAILU RE);
}
/* 2D memory */
/* Allocate memory using MACRO */
printf("Allocat ing 2D memory\n");
work.dayA = malloc(sizeof(i nt*)*AreaRow);
work.twilightA = malloc(sizeof(i nt*)*AreaRow);
work.ppsFCA = malloc(sizeof(i nt*)*AreaRow);
work.msgFCA = malloc(sizeof(i nt*)*AreaRow);
work.msgcloudyA = malloc(sizeof(i nt*)*AreaRow);
work.ppscloudyA = malloc(sizeof(i nt*)*AreaRow);
work.lightconA = malloc(sizeof(f loat*)*AreaRow) ;
work.bias100A = malloc(sizeof(f loat*)*AreaRow) ;
work.bias75A = malloc(sizeof(f loat*)*AreaRow) ;
work.bias50A = malloc(sizeof(f loat*)*AreaRow) ;
work.bias25A = malloc(sizeof(f loat*)*AreaRow) ;
work.valconcenA = malloc(sizeof(i nt*)*AreaRow);
for (i = 0; i < AreaRow; i++) {
work.dayA[i] = malloc(sizeof(i nt)*AreaCol);
work.twilightA[i] = malloc(sizeof(i nt)*AreaCol);
work.ppsFCA[i] = malloc(sizeof(i nt)*AreaCol);
work.msgFCA[i] = malloc(sizeof(i nt)*AreaCol);
work.msgcloudyA[i] = malloc(sizeof(i nt)*AreaCol);
work.ppscloudyA[i] = malloc(sizeof(i nt)*AreaCol);
work.lightconA[i] = malloc(sizeof(f loat)*AreaCol);
work.bias100A[i] = malloc(sizeof(f loat)*AreaCol);
work.bias75A[i] = malloc(sizeof(f loat)*AreaCol);
work.bias50A[i] = malloc(sizeof(f loat)*AreaCol);
work.bias25A[i] = malloc(sizeof(f loat)*AreaCol);
work.valconcenA[i] = malloc(sizeof(i nt)*AreaCol);
}
if (work.dayA==NUL L) {
printf("Failed to allocate memory: dayA\n");
exit(EXIT_FAILU RE);
}
if (work.twilightA ==NULL) {
printf("Failed to allocate memory: twilightA\n");
exit(EXIT_FAILU RE);
}
if (work.ppsFCA==N ULL) {
printf("Failed to allocate memory: ppsFCA\n");
exit(EXIT_FAILU RE);
}
if (work.msgFCA==N ULL) {
printf("Failed to allocate memory: msgFCA\n");
exit(EXIT_FAILU RE);
}
if (work.msgcloudy A==NULL) {
printf("Failed to allocate msgcloudyA\n");
exit(EXIT_FAILU RE);
}
if (work.ppscloudy A==NULL) {
printf("Failed to allocate memory: ppscloudyA\n");
exit(EXIT_FAILU RE);
}
if (work.lightconA ==NULL) {
printf("Failed to allocate memory: lightconA\n");
exit(EXIT_FAILU RE);
}
if (work.bias100A= =NULL) {
printf("Failed to allocate memory: bias100A\n");
exit(EXIT_FAILU RE);
}
if (work.bias75A== NULL) {
printf("Failed to allocate memory: bias75A\n");
exit(EXIT_FAILU RE);
}
if (work.bias50A== NULL) {
printf("Failed to allocate memory: bias50A\n");
exit(EXIT_FAILU RE);
}
if (work.bias25A== NULL) {
printf("Failed to allocate memory: bias25A\n");
exit(EXIT_FAILU RE);
}
if (work.valconcen A==NULL) {
printf("Failed to allocate memory valA\n");
exit(EXIT_FAILU RE);
}
printf("Allocat ed memory for all arrays\n");
return 1;
}

static int MemoryFreeA(voi d) {
int i;

for (i = 0; i < AreaCol; i++) {
free(work.dayA[i]);
free(work.twili ghtA[i]);
free(work.msgFC A[i]);
free(work.ppsFC A[i]);
free(work.msgcl oudyA[i]);
free(work.ppscl oudyA[i]);
free(work.light conA[i]);
free(work.bias1 00A[i]);
free(work.bias7 5A[i]);
free(work.bias5 0A[i]);
free(work.bias2 5A[i]);
free(work.valco ncenA[i]);
}
free(work.dayA) ;
free(work.twili ghtA);
free(work.msgFC A);
free(work.ppsFC A);
free(work.msgcl oudyA);
free(work.ppscl oudyA);
free(work.light conA);
free(work.bias1 00A);
free(work.bias7 5A);
free(work.bias5 0A);
free(work.bias2 5A);
free(work.valco ncenA);

return 1;
}

static int MemoryFree15km( void) {

int i;

for (i = 0; i < ResCol; i++) {
free(work.light con15km[i]);
free(work.bias1 0015km[i]);
free(work.bias7 515km[i]);
free(work.bias5 015km[i]);
free(work.bias2 515km[i]);
free(work.valco ncen15km[i]);
}
free(work.light con15km);
free(work.bias1 0015km);
free(work.bias7 515km);
free(work.bias5 015km);
free(work.bias2 515km);
free(work.valco ncen15km);

return 1;
}

static int MemoryFreeOrd(v oid) {

int i;

for (i = 0; i < work.Col; i++) {
free(work.day[i]);
free(work.twili ght[i]);
free(work.msgcl oudy[i]);
free(work.msgFC[i]);
free(work.ppsFC[i]);
free(work.ppscl oudy[i]);
free(work.msgda ta[i]);
free(work.ppsda ta[i]);
free(work.valco ncen[i]);
free(work.frac[i]);
free(work.qual[i]);
}
free(work.day);
free(work.twili ght);
free(work.msgcl oudy);
free(work.msgFC );
free(work.ppsFC );
free(work.ppscl oudy);
free(work.msgda ta);
free(work.ppsda ta);
free(work.valco ncen);
free(work.frac) ;
free(work.qual) ;

free(work.date) ;
free(work.tile) ;
free(work.Regio n);
free(work.Repro dir);
free(work.PPSfi le);
free(work.MSGfi le);
free(work.PHYSf ile);

return 1;
}

static int MemoryFree1D(vo id) {

int i;

free(work.stats );
free(work.nscen es);
free(work.perce ntage);

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_ar ea.c:1675
#3 0x405dae0a in msgppsarea (self=0x0, args=0x4042aa7c , kwargs=0x0) at
msgppsmodule_ar ea.c:328
#4 0x40073b16 in PyCFunction_Cal l () 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_ar ea.c:1675
#3 0x405dae0a in msgppsarea (self=0x0, args=0x4042aa7c , kwargs=0x0) at
msgppsmodule_ar ea.c:328
#4 0x40073b16 in PyCFunction_Cal l () 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_ar ea.c:1675
#3 0x405dae0a in msgppsarea (self=0x0, args=0x4042aa7c , kwargs=0x0) at
msgppsmodule_ar ea.c:328
#4 0x40073b16 in PyCFunction_Cal l () 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 readPythonObjec t(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_AsStri ng(msgop);
ppsop = PyList_GetItem( work.ppsobj, i);
work.pps_scenes[i] = PyString_AsStri ng(ppsop);
}
for (i = 0; i < NumberOfTiles; i++) {
tileop = PyList_GetItem( work.tileobj, i);
work.tiles[i] = PyString_AsStri ng(tileop);
sceneop = PyList_GetItem( work.nscenesobj , i);
work.nscenes[i] = PyInt_AsLong(sc eneop);
}
return 1;
} /*end readPythonObjec t*/

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

This thread has been closed and replies have been disabled. Please start a new discussion.

Similar topics

3
6209
by: Nick Craig-Wood | last post by:
I've just discovered that my python (Python 2.3.4 from debian package 2.3.4-1 running on debian testing x86 + linux 2.4.26) core dumps when I set recursionlimit very high and do lots of recursion. Eg $ python -c 'import sys; sys.setrecursionlimit(100000) def f(): return f() f()' Segmentation fault (core dumped)
1
3512
by: Martin | last post by:
I use dbx and i got the following error: Reading GL_CliConnMgr core file header read successfully Reading ld.so.1 dbx: core file read error: address 0xff3e6000 not available dbx: core file read error: address 0xff3e70a0 not available dbx: warning: could not initialize librtld_db.so.1 -- trying libDP_rtld_db.so Make sure this is the same version of Solaris where the core dump
10
4026
by: ken | last post by:
hello, i'm writing a c program on a linux system. i'm debugging a segmentation fault but i don't want it to dump a core file because the memory footprint of the program is over 300Mb and i don't need it to generate a 300Mb file every time I add a new printf statement to debug the code. can i do something to prevent it from dumping the core file even when it seg faults? (is this a unix/linux thing, or a c thing?) thanks!
1
1693
by: invincible | last post by:
hi, assume my program is running , is there a way I can dump its core in solaris, in linux i think u can use gcore Mohan
10
584
by: John Liu | last post by:
We upgraded from 7.2 to 7.4, it looks like everything working, but when I issue a query such as select * from tab (tab has about 2-3 million records), it causes core dump. I tuned some the parameters, it still produce the core. Thanks for any hints.
3
2459
by: John Liu | last post by:
AIX pg version 7.4 Select * from document2 core dump. Did a few more experiments with select * from document2 limit... I limit to 500000 it works, 600000 it exits but says "calloc: There is not
4
1944
by: madhusudan.hv | last post by:
hi, Can you tell me when an application might take long time to dump core? The linux o.s is indicating that a process is dumping core, and only after 40 mins i am seeing the core file. What might be the problem? regards, Madhusudan
10
7447
by: wong_powah | last post by:
I want to find out where (which line) my C program core dump. How to do that? Is there a web site describing the procedure? One approach is to use stack trace of the mdb debugger, but I does not understand its output completely. e.g. How to interpret the stack trace to find out which line inside the coredump_func function of a test program caused the
5
6096
by: johnericaturnbull | last post by:
Hi - I am very new to python. I get this random core dump and am looking for a good way to catch the error. I know the function my core dump occurs. Is there any error catching/handling that I could use in python?
0
9564
by: Hystou | last post by:
Most computers default to English, but sometimes we require a different language, especially when relocating. Forgot to request a specific language before your computer shipped? No problem! You can effortlessly switch the default language on Windows 10 without reinstalling. I'll walk you through it. First, let's disable language synchronization. With a Microsoft account, language settings sync across devices. To prevent any complications,...
0
10548
Oralloy
by: Oralloy | last post by:
Hello folks, I am unable to find appropriate documentation on the type promotion of bit-fields when using the generalised comparison operator "<=>". The problem is that using the GNU compilers, it seems that the internal comparison operator "<=>" tries to promote arguments from unsigned to signed. This is as boiled down as I can make it. Here is my compilation command: g++-12 -std=c++20 -Wnarrowing bit_field.cpp Here is the code in...
1
10295
by: Hystou | last post by:
Overview: Windows 11 and 10 have less user interface control over operating system update behaviour than previous versions of Windows. In Windows 11 and 10, there is no way to turn off the Windows Update option using the Control Panel or Settings app; it automatically checks for updates and installs any it finds, whether you like it or not. For most users, this new feature is actually very convenient. If you want to control the update process,...
0
10069
tracyyun
by: tracyyun | last post by:
Dear forum friends, With the development of smart home technology, a variety of wireless communication protocols have appeared on the market, such as Zigbee, Z-Wave, Wi-Fi, Bluetooth, etc. Each protocol has its own unique characteristics and advantages, but as a user who is planning to build a smart home system, I am a bit confused by the choice of these technologies. I'm particularly interested in Zigbee because I've heard it does some...
1
7604
isladogs
by: isladogs | last post by:
The next Access Europe User Group meeting will be on Wednesday 1 May 2024 starting at 18:00 UK time (6PM UTC+1) and finishing by 19:30 (7.30PM). In this session, we are pleased to welcome a new presenter, Adolph Dupré who will be discussing some powerful techniques for using class modules. He will explain when you may want to use classes instead of User Defined Types (UDT). For example, to manage the data in unbound forms. Adolph will...
0
6842
by: conductexam | last post by:
I have .net C# application in which I am extracting data from word file and save it in database particularly. To store word all data as it is I am converting the whole word file firstly in HTML and then checking html paragraph one by one. At the time of converting from word file to html my equations which are in the word document file was convert into image. Globals.ThisAddIn.Application.ActiveDocument.Select();...
0
5500
by: TSSRALBI | last post by:
Hello I'm a network technician in training and I need your help. I am currently learning how to create and manage the different types of VPNs and I have a question about LAN-to-LAN VPNs. The last exercise I practiced was to create a LAN-to-LAN VPN between two Pfsense firewalls, by using IPSEC protocols. I succeeded, with both firewalls in the same network. But I'm wondering if it's possible to do the same thing, with 2 Pfsense firewalls...
1
4275
by: 6302768590 | last post by:
Hai team i want code for transfer the data from one system to another through IP address by using C# our system has to for every 5mins then we have to update the data what the data is updated we have to send another system
2
3798
muto222
by: muto222 | last post by:
How can i add a mobile payment intergratation into php mysql website.

By using Bytes.com and it's services, you agree to our Privacy Policy and Terms of Use.

To disable or enable advertisements and analytics tracking please visit the manage ads & tracking page.