/[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 1773 - (show annotations) (download) (as text)
Mon May 19 02:52:08 2008 UTC (10 years, 6 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 /* 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 char code[3];
107 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 CONSOLE_DEBUG("Comparing code with '%s'",i->code);
208 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 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 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 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 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
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 */
312
313 #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
320 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
345 char field[10];
346 char code[3];
347
348 unsigned i;
349 assert(N_FIELDS == sizeof(fieldsize) / sizeof(unsigned));
350
351 fgets(code, 3, d->f);
352 //CONSOLE_DEBUG("code = '%s'",code);
353 //assert(strcmp(code,"CA")==0);
354
355 /*
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
365 assert(field!=NULL);
366 assert(strlen(field)!=0);
367
368 data[i] = atoi(field);
369 }
370
371 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 /*
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 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
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