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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

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