/[ascend]/trunk/models/johnpye/datareader/acdb.c
ViewVC logotype

Annotation of /trunk/models/johnpye/datareader/acdb.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1773 - (hide annotations) (download) (as text)
Mon May 19 02:52:08 2008 UTC (16 years, 2 months ago) by jpye
File MIME type: text/x-csrc
File size: 13629 byte(s)
More on ACDB data reader. Reads data OK now, need to check the values are coming through correctly.
1 jpye 1772 /* ASCEND modelling environment
2     Copyright (C) 2006 Carnegie Mellon University
3    
4     This program is free software; you can redistribute it and/or modify
5     it under the terms of the GNU General Public License as published by
6     the Free Software Foundation; either version 2, or (at your option)
7     any later version.
8    
9     This program is distributed in the hope that it will be useful,
10     but WITHOUT ANY WARRANTY; without even the implied warranty of
11     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12     GNU General Public License for more details.
13    
14     You should have received a copy of the GNU General Public License
15     along with this program; if not, write to the Free Software
16     Foundation, Inc., 59 Temple Place - Suite 330,
17     Boston, MA 02111-1307, USA.
18     *//**
19     @file
20     Data Reader implementation for the Australian Climate Data Bank .DAT format.
21    
22     These functions implement a reader interface for the Australian Climate
23     Data Bank meteorological data as given in .DAT files from the Australian
24     bureau of Meteorology.
25    
26     http://www.bom.gov.au/climate/averages/tables/supply.shtml
27    
28     The following comes from the README.txt file supplied by the BoM, to explain
29     this data format:
30    
31     Characters Item
32    
33     1 - 2 location identification (e.g.ME represents Melbourne)
34     3 - 4 year (e.g. 67)
35     5 - 6 month (i.e. 1 - 12)
36     7 - 8 day (i.e. 1 - 31)
37     9 - 10 hour standard (i.e. 0-23, 0 = midnight)
38     11 - 14 dry bulb temperature (10-1 ��C)
39     15 - 17 absolute moisture content (10-1 g/kg)
40     18 - 21 atmospheric pressure (10-1 kPa)
41     22 - 24 wind speed (10-1 m/s)
42     25 - 26 wind direction (0-16; 0 = CALM. 1 = NNE ,...,16 = N
43     27 total cloud cover (oktas, 0 - 8)
44     28 flag relating to dry bulb temperature
45     29 flag relating to absolute moisture content
46     30 flag relating to atmospheric pressure
47     31 flag relating to wind speed and direction
48     32 flag relating to total cloud cover
49     33 blank
50     34 - 37 global solar irradiance on a horizontal plane (W/m2)
51     38 - 40 diffuse solar irradiance on a horizontal plane (W/m2)
52     41 - 44 direct solar irradiance on a plane normal to the beam ((W/m2)
53     45 - 46 solar altitude (degrees, 0-90)
54     47 - 49 solar azimuth (degrees, 0-360)
55     50 flag relating to global and diffuse solar irradiance
56     51 flag }
57     52 - 56 Australian Met Station Number } Some locations only
58     57 - 61 wet bulb temperature (10-1 ��C) }
59     62 - 81 Station name (first line only) }
60    
61    
62     Values for flags relating to standard surface meteorological data
63     columns 28 - 32)
64    
65     0 means that the value is measured value
66     1 means that the value is estimated to replace a missing measurement
67     2 means that the value is an interpolating between three-hourly measurements
68     3 missing value
69    
70     Values for flag relating to solar radiation data (column 50)
71    
72     0 means that both global and diffuse irradiance values are based on measurements
73     1 means that both global and diffuse irradiance values are estimated to
74     replace a missing measurement
75     2 means that the global irradiance value is based on measurement but the
76     diffuse irradiance value is estimated to replace a missing measurement
77     3 missing value or estimated value from cloud cover data
78     4 interpolated value from three hourly data
79     *//*
80     by John Pye, Aug 2006
81     */
82    
83     #include <stdio.h>
84     /* #include <libradtran/sun.h> */
85     #include "sun.h"
86    
87     #include <utilities/ascMalloc.h>
88     #include <utilities/error.h>
89    
90     #include "tmy.h"
91    
92    
93     typedef struct AcdbPoint_struct{
94     double t;
95     float I;
96     float Ibn;
97     float Id;
98     float T;
99     float v_wind;
100    
101     } AcdbPoint;
102    
103     #define DATA(D) ((AcdbPoint *)(D->data))[D->i]
104    
105     struct AcdbCity{
106 jpye 1773 char code[3];
107 jpye 1772 char name[40];
108     char state[40];
109     char dir[9];
110     };
111    
112     static const struct AcdbCity acdb_city_info[] = {
113     {"CN","Cairns AMO","QLD","CAIRNS"}
114     ,{"CV","Charleville AMO","QLD","CHARLEVI"}
115     ,{"CL","Cloncurry AMO","QLD","CLONCURR"}
116     ,{"GL","Gladstone MO","QLD","GLADSTON"}
117     ,{"LO","Longreach","QLD","LONGREAC"}
118     ,{"MK","Mackay MO","QLD","MACKAY"}
119     ,{"IS","Mt Isa AMO","QLD","MTISA"}
120     ,{"OA","Oakey Army Aviation","QLD","OAKEY"}
121     ,{"RO","Rockhampton AMO","QLD","ROCKHPTN"}
122     ,{"TO","Townsville AMO","QLD","TOWNSVL"}
123     ,{"WS","Wllis Island MO","QLD","WILLISIS"}
124     ,{"CO","Cobar AMO","NSW","COBAR"}
125     ,{"CH","Coffs Harbour MO","NSW","COFFSHAR"}
126     ,{"MA","Mascot (Sydney AMO)","NSW","MASCOT"}
127     ,{"MO","Moree MO","NSW","MOREE"}
128     ,{"NO","Nowra RAN AIR STN","NSW","NOWRA"}
129     ,{"OR","Orange AP","NSW","ORANGE"}
130     ,{"RI","Richmond","NSW","RICHMOND"}
131     ,{"SY","Sydney RO","NSW","SYDNEYRO"}
132     ,{"TM","Tamworth AMO","NSW","TAMWORTH"}
133     ,{"WA","Wagga Wagga AMO","NSW","WAGGA"}
134     ,{"WE","Williamtown AMO","NSW","WILLIAM"}
135     ,{"SE","East Sale AMO","VIC","EASTSALE"}
136     ,{"ES","Essendon AMO","VIC","ESSENDON"}
137     ,{"LA","Laverton AMO","VIC","LAVERTON"}
138     ,{"ME","Melbourne RO","VIC","MELBOURN"}
139     ,{"MI","Mildura AMO","VIC","MILDURA"}
140     ,{"NH","Nhill COMPOSITE","VIC","NHILL"}
141     ,{"TU","Tullamarine","VIC","TULLAMAR"}
142     ,{"HB","Hobart AMO","TAS","HOBRTAMO"}
143     ,{"HO","Hobart RO","TAS","HOBARTRO"}
144     ,{"LU","Launceston AP","TAS","LAUNCESN"}
145     ,{"AD","Adelaide AP","SA","ADELAP"}
146     ,{"AR","Adelaide RO","SA","ADELRO"}
147     ,{"AW","Adelaide West Terrace","SA","ADELWEST"}
148     ,{"CE","Ceduna","SA","CEDUNA"}
149     ,{"MG","Mt Gambier AMO","SA","MTGAMB"}
150     ,{"OO","Oodnadatta AMO","SA","OODNADAT"}
151     ,{"WO","Woomera AMO","SA","WOOMERA"}
152     ,{"AL","Alice Springs","NT","ALICESPR"}
153     ,{"DA","Darwin AP","NT","DARWINAP"}
154     ,{"DR","Darwin RO","NT","DARWINRO"}
155     ,{"TE","Tennant Crk MO","NT","TENNANT"}
156     ,{"AY","Albany (Eclipse Is)","WA","ALBANY"}
157     ,{"AB","Albany AMO","WA","ALBANAMO"}
158     ,{"BM","Broome AMO","WA","BROOME"}
159     ,{"CR","Carnarvon AMO","WA","CARNARVO"}
160     ,{"DE","Derby","WA","DERBY"}
161     ,{"EP","Esperance","WA","ESPERANC"}
162     ,{"FO","Forrest AMO","WA","FORREST"}
163     ,{"GE","Geraldton AMO","WA","GERALDTO"}
164     ,{"GI","Giles","WA","GILES"}
165     ,{"HA","Halls Creek","WA","HALLSCRK"}
166     ,{"KA","Kalgoorlie","WA","KALGOORL"}
167     ,{"LM","Learmonth MO","WA","LEARMTH"}
168     ,{"MT","Meekatharra","WA","MEEKATH"}
169     ,{"OW","Onslow","WA","ONSLOW"}
170     ,{"ON","Onslow PO","WA","ONSLOWPO"}
171     ,{"PR","Perth AMO","WA","PERTHAMO"}
172     ,{"PE","Perth RO","WA","PERTHRO"}
173     ,{"HE","Pt Hedland","WA","PORTHEDL"}
174     ,{"CB","Canberra AMO","ACT","CANBAMO"}
175     ,{"CA","Canberra CITY","ACT","CANBERRA"}
176     ,{"LE","Lae","New Guinea","LAE"}
177     ,{"PM","Pt Moresby Jackson","New Guinea","PTMORESB"}
178     ,{"RB","Rabual","New Guinea","RABAUL"}
179     ,{"LH","Lord Howe Island","Aust. Islands","LORDHOWE"}
180     ,{"NI","Norfolk Island AMO","Aust. Islands","NORFOLK"}
181     ,{"CI","Cocos Island COMP","Aust. Islands","COCOSIS"}
182     ,{"HH","Honiara Henderson AP","Solomon Islands","HONIARAP"}
183     ,{"HN","Honiara MO","Solomon Islands","HONIARMO"}
184     ,{"AN","Aneityum AERO","Vanuatu","ANEITYUM"}
185     ,{"VL","Vanua Lava","Vanuatu","VANUA"}
186     ,{"VI","Vila","Vanuatu","VILA"}
187     ,{"BT","Butterworth RAAF","Malaysia","BUTTRAAF"}
188     ,{"","","",""}
189     };
190    
191     /**
192     @return 0 on success
193     */
194     int datareader_acdb_header(DataReader *d){
195    
196     char code[3];
197     unsigned yr;
198     unsigned data_rows;
199    
200     fscanf(d->f,"%2c%2ud",code,&yr);
201    
202     code[2] = '\0';
203    
204     const struct AcdbCity *i;
205     unsigned found=0;
206     for(i=acdb_city_info; i->code[0] != '\0'; ++i){
207 jpye 1773 CONSOLE_DEBUG("Comparing code with '%s'",i->code);
208 jpye 1772 if(strcmp(i->code,code)==0){
209     found=1;
210     break;
211     }
212     }
213     if(!found){
214     CONSOLE_DEBUG("Unknown city '%s' in ACDB data file",code);
215     ERROR_REPORTER_HERE(ASC_PROG_WARNING,"ACDB data file contains unrecognised city code '%s'",code);
216     }else{
217     CONSOLE_DEBUG("ACDB data file contains data for %s, %s.",i->name, i->state);
218     ERROR_REPORTER_HERE(ASC_PROG_NOTE,"ACDB data file contains data for %s, %s.",i->name, i->state);
219     }
220    
221     if(yr < 50)yr+=2000;
222     else yr+=1900;
223    
224     if(is_leap_year(yr)){
225     data_rows = 366 * 24;
226     }else{
227     data_rows = 365 * 24;
228     }
229 jpye 1773 CONSOLE_DEBUG("ACDB data file is for year %u, expect %u data rows.",yr,data_rows);
230     ERROR_REPORTER_HERE(ASC_PROG_NOTE,"ACDB data file is for year %u, expect %u data rows.",yr,data_rows);
231    
232 jpye 1772 d->ndata = data_rows;
233     d->i = 0;
234     d->ndata=8760;
235     d->data = ASC_NEW_ARRAY(AcdbPoint,d->ndata);
236    
237     /* every line contains the city name, so rewind the file now */
238     rewind(d->f);
239    
240     return 0;
241     }
242    
243     int datareader_acdb_eof(DataReader *d){
244     if(feof(d->f)){
245     CONSOLE_DEBUG("REACHED END OF FILE");
246 jpye 1773 if(d->i < d->ndata){
247     ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Incomplete data set found (%d rows < %d expected",d->i, d->ndata);
248     }
249 jpye 1772 d->ndata=d->i;
250     ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Read %d rows",d->ndata);
251     return 1;
252     }
253    
254     /* set the number of inputs and outputs */
255     d->ninputs = 1;
256     d->noutputs = 5;
257     return 0;
258     }
259    
260     #define MEAS(N) int N; char N##_source; int N##_uncert
261     #define READ(N) &N, &N##_source, &N##_uncert
262    
263     /**
264     Read a line of data and store in d.
265     @return 0 on success
266     */
267     int datareader_acdb_data(DataReader *d){
268     /* static int lastmonth=-1;
269     static int lastday=-1; */
270    
271     AcdbPoint *dat;
272    
273     /* still need to implement...
274    
275     51 flag }
276     52 - 56 Australian Met Station Number } Some locations only
277     57 - 61 wet bulb temperature (10-1 ��C) }
278     62 - 81 Station name (first line only) }
279    
280 jpye 1773
281    
282     1 - 2 location identification (e.g.ME represents Melbourne)
283     3 - 4 year (e.g. 67)
284     5 - 6 month (i.e. 1 - 12)
285     7 - 8 day (i.e. 1 - 31)
286     9 - 10 hour standard (i.e. 0-23, 0 = midnight)
287    
288     11 - 14 dry bulb temperature (10-1 ��C)
289     15 - 17 absolute moisture content (10-1 g/kg)
290     18 - 21 atmospheric pressure (10-1 kPa)
291     22 - 24 wind speed (10-1 m/s)
292     25 - 26 wind direction (0-16; 0 = CALM. 1 = NNE ,...,16 = N
293     27 total cloud cover (oktas, 0 - 8)
294     28 flag relating to dry bulb temperature
295     29 flag relating to absolute moisture content
296     30 flag relating to atmospheric pressure
297     31 flag relating to wind speed and direction
298     32 flag relating to total cloud cover
299     33 blank
300     34 - 37 global solar irradiance on a horizontal plane (W/m2)
301     38 - 40 diffuse solar irradiance on a horizontal plane (W/m2)
302     41 - 44 direct solar irradiance on a plane normal to the beam ((W/m2)
303     45 - 46 solar altitude (degrees, 0-90)
304     47 - 49 solar azimuth (degrees, 0-360)
305     50 flag relating to global and diffuse solar irradiance
306     51 flag }
307     52 - 56 Australian Met Station Number } Some locations only
308     57 - 61 wet bulb temperature (10-1 ��C) }
309     62 - 81 Station name (first line only) }
310    
311 jpye 1772 */
312    
313 jpye 1773 #define N_FIELDS 22
314     const unsigned fieldsize[N_FIELDS] = {
315     2,2,2,2, 4,3,4,3,2,1, 1,1,1,1,1, 1, 4,3,4,3,3,1
316     };
317    
318     int data[N_FIELDS];
319 jpye 1772
320 jpye 1773 enum {
321     ACDB_YEAR
322     ,ACDB_MONTH
323     ,ACDB_DAY
324     ,ACDB_HOUR
325     ,ACDB_DRY_BULB_TEMP/* [0.1��C] */
326     ,ACDB_W /* ABSOLUTE MOISTURE CONTENT [0.1 G/KG] */
327     ,ACDB_P_ATM /* ATMOSPHERIC PRESSURE [0.1 KPA] */
328     ,ACDB_V_WIND /* WIND SPEED [0.1 M/S] */
329     ,ACDB_DIR_WIND /* DIR OF WIND: 0-16; 0 = CALM. 1 = NNE ,...,16 = N */
330     ,ACDB_CLOUD /* TOTAL CLOUD COVER (OKTAS, 0 - 8) */
331     ,ACDB_DRY_BULB_TEMP_FLAG
332     ,ACDB_W_FLAG
333     ,ACDB_P_ATM_FLAG
334     ,ACDB_WIND_FLAG
335     ,ACDB_CLOUD_FLAG
336     ,ACDB_WHITESPACE1
337     ,ACDB_GHI
338     ,ACDB_IHI /* indirect (diffuse) horizontal radiation */
339     ,ACDB_DNI
340     ,ACDB_ALTITUDE_DEG
341     ,ACDB_AZIMUTH_DEG
342     ,ACDB_SOLAR_FLAG
343     };
344 jpye 1772
345 jpye 1773 char field[10];
346     char code[3];
347    
348     unsigned i;
349     assert(N_FIELDS == sizeof(fieldsize) / sizeof(unsigned));
350 jpye 1772
351 jpye 1773 fgets(code, 3, d->f);
352     //CONSOLE_DEBUG("code = '%s'",code);
353     //assert(strcmp(code,"CA")==0);
354 jpye 1772
355 jpye 1773 /*
356     this format includes spaces in the data file, so we can't use fscanf
357     because of the way it soaks up leading spaces without counting them
358     in the field width.
359     */
360     for(i=0; i < N_FIELDS; ++i){
361     fgets(field, fieldsize[i] + 1, d->f);
362     //CONSOLE_DEBUG("field %d: size = %d, str = '%s'", i, fieldsize[i], field);
363     if(i==ACDB_WHITESPACE1)continue;
364 jpye 1772
365 jpye 1773 assert(field!=NULL);
366     assert(strlen(field)!=0);
367    
368     data[i] = atoi(field);
369 jpye 1772 }
370    
371 jpye 1773 fscanf(d->f," ");
372    
373     CONSOLE_DEBUG("Time: %d/%d/%d %2d:00",data[ACDB_DAY],data[ACDB_MONTH],data[ACDB_YEAR],data[ACDB_HOUR]);
374     CONSOLE_DEBUG("code = %s, dry_bulb_temp = %f, w = %f, p_atm = %f, v_wind = %f, dir_wind = %d"
375     ,code, data[ACDB_DRY_BULB_TEMP]/10., data[ACDB_W]/10., data[ACDB_P_ATM]/10., data[ACDB_V_WIND]/10., data[ACDB_DIR_WIND]
376     );
377     assert(strcmp(code,"CA")==0);
378     assert(data[ACDB_DRY_BULB_TEMP]/10. <= 70.);
379     assert(data[ACDB_DIR_WIND] >= 0 && data[ACDB_DIR_WIND] <= 16);
380    
381 jpye 1772 /*
382     if(month!=lastmonth || day!=lastday){
383     CONSOLE_DEBUG("Reading data for %d/%d",day,month);
384     lastmonth=month;
385     lastday=day;
386     }
387     */
388    
389     /*
390     for the moment, we only record global horizontal, direct normal,
391     ambient temperature, wind speed.
392     */
393    
394     dat = &DATA(d);
395 jpye 1773 dat->t = ((day_of_year_specific(data[ACDB_DAY],data[ACDB_MONTH],data[ACDB_YEAR]) - 1)*24.0 + data[ACDB_HOUR])*3600.0;
396     dat->I = data[ACDB_GHI] * 0.1; /* average W/m2 for the hour in question */
397     dat->Ibn = data[ACDB_DNI] * 0.1; /* normal beam radiation, W/m2 */
398     dat->Id = data[ACDB_IHI] * 0.1;
399     dat->T = 0.1*data[ACDB_DRY_BULB_TEMP] + 273.15; /* temperature */
400     dat->v_wind = data[ACDB_V_WIND] * 0.1;
401 jpye 1772
402     #if 0
403     if(d->i < 20){
404     CONSOLE_DEBUG("ROW %d: year %d, month %d, day %d, hour %d (t = %.0f), Iegh %d, Iedn %d, Igh (%c) = %d --> I = %f"
405     , d->i, year, month, day, hour, tmy->t, Iegh, Iedn, Igh_source, Igh, tmy->I
406     );
407     }
408     #endif
409    
410     d->i++;
411     return 0;
412     }
413    
414     int datareader_acdb_time(DataReader *d, double *t){
415     *t = DATA(d).t;
416     return 0;
417     }
418    
419     int datareader_acdb_vals(DataReader *d, double *v){
420     CONSOLE_DEBUG("At t=%f h, T = %lf, I = %f J/m2"
421     ,(DATA(d).t/3600.0),DATA(d).T, DATA(d).I
422     );
423     v[0]=DATA(d).I;
424     v[1]=DATA(d).Ibn;
425     v[2]=DATA(d).Id;
426     v[3]=DATA(d).T;
427     v[4]=DATA(d).v_wind;
428     return 0;
429     }

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22