johnpye/datareader/tmy.c

/*  ASCEND modelling environment
    Copyright (C) 2006 Carnegie Mellon University

    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2, or (at your option)
    any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place - Suite 330,
    Boston, MA 02111-1307, USA.
*//**
    @file
    Data Reader implementation for the TMY2 format.

    These functions implement a reader interface for meteorological data in the
    TMY2 format as specified at http://rredc.nrel.gov/solar/pubs/tmy2/tab3-2.html
*//*
    by John Pye, Aug 2006
*/

#include <stdio.h>
/* #include <libradtran/sun.h> */
#include "sun.h"

#include <utilities/ascMalloc.h>
#include <utilities/error.h>

#include "tmy.h"


typedef struct{
    double t;
    float G;
    float Gbn;
    float T;
    float v_wind;
    
} Tmy2Point;

#define DATA(D) ((Tmy2Point *)(D->data))[D->i]

/**
    @return 0 on success
*/
int datareader_tmy2_header(DataReader *d){
    char wban[-2 + 6 +2];
    char city[-8 +29 +2];
    char state[-31+32+2];
    int zone;
    char lathemi;
    int latdeg, latmin;
    char longhemi;
    int longdeg, longmin;
    int elev;

    fscanf(d->f,"%s %s %s %d" 
        " %c %d %d"
        " %c %d %d"
        " %d"
            ,wban,city,state,&zone
            ,&lathemi,&latdeg,&latmin
            ,&longhemi,&longdeg,&longmin
            ,&elev
    );

    double lat = latdeg + latmin/60;
    if(lathemi=='S')lat=-lat;
    double lng = longdeg + longmin/60;
    if(longhemi=='E')lng=-lng;
    CONSOLE_DEBUG( "TMY2 data for city '%s' (WBAN %s, time zone %+d) at lat=%.3f, long=%.3f, elev=%d m"
        , city, wban, zone, lat, lng, elev
    );
    ERROR_REPORTER_HERE(ASC_PROG_NOTE,"TMY2 data is for %s, %s",city, state);

    d->i = 0;
    d->ndata=8760;
    d->data = ASC_NEW_ARRAY(Tmy2Point,d->ndata);
    return 0;
}

int datareader_tmy2_eof(DataReader *d){
    if(feof(d->f)){
        CONSOLE_DEBUG("REACHED END OF FILE");
        d->ndata=d->i;
        ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Read %d rows",d->ndata);
        return 1;
    }

    /* set the number of inputs and outputs */
    d->ninputs = 1;
    d->noutputs = 2;
    return 0;
}

#define MEAS(N) int N; char N##_source; int N##_uncert
#define READ(N) &N, &N##_source, &N##_uncert
    
/**
    Read a line of data and store in d.
    @return 0 on success
*/
int datareader_tmy2_data(DataReader *d){
    /* static int lastmonth=-1;
    static int lastday=-1; */
    int res = 0;

    Tmy2Point *tmy;
    int year,month,day,hour;
    int Iegh,Iedn; // Irradiation
    MEAS(Igh); /* Global horizontal irradiation in the specified interval / (Wh/m2) */
    MEAS(Idn); /* Direct normal irradiation in the specified interval / (Wh/m2) */
    MEAS(Idh); /* Diffuse horizontal irradiation in the specified interval / (Wh/m2) */

    MEAS(Lgh); /* Global horiz illuminance / (100 lux) */
    MEAS(Ldn); /* Direct normal illuminance / (100 lux) */
    MEAS(Ldh); /* Diffuse horiz illuminance / (100 lux) */
    MEAS(Lz); /* Zenith illuminance / (10 Cd/m2) */

    MEAS(covtot); /* Total sky cover / tenths */
    MEAS(covopq); /* Opaque sky cover / tenths */

    MEAS(T); /* temperature / (0.1degC) */
    MEAS(Tdew); /* dew point temperature / (0.1degC) */
    MEAS(p); /* pressure / mbar */
    MEAS(rh); /* rel humidity / % */
    MEAS(wdir); /* wind dir, N=0, E=90,... */
    MEAS(wvel); /* wind speed / (m/s) */
    MEAS(vis); /* visibility / (100m) */
    MEAS(ch); /* ceiling height / m (or special value) */
    
    MEAS(rain); /* preciptable water / mm */
    MEAS(aer); /* aerosol optical depth in thousandths (???) */
    MEAS(snow); /* snow depth on the specified day / cm (999=missing data) */
    MEAS(dsno); /* days since last snow (or special value) */

    /* weather observations from Appendix B */
    int obs, storm, precip, drizz, snowtype, snowshower, sleet, fog, smog, hail;

    /* brace yourself for this one... */

    res = fscanf(d->f, 
        /* 1 */ "%2d%2d%2d%2d" "%4d%4d" "%4d%1s%1d" "%4d%1s%1d" "%4d%1s%1d" /* =15 */
        /* 2 */ "%4d%1s%1d" "%4d%1s%1d" "%4d%1s%1d" "%4d%1s%1d" /* +12=27 */
        /* 3 */ "%2d%1s%1d" "%2d%1s%1d" "%4d%1s%1d" "%4d%1s%1d" "%3d%1s%1d" "%4d%1s%1d" /* +18=45 */
        /* 4 */ "%3d%1s%1d" "%3d%1s%1d" "%4d%1s%1d" "%5d%1s%1d" /* +12=57 */
        /* 5 */ "%1d%1d%1d%1d%1d%1d%1d%1d%1d%1d" /* +10=67 */
        /* 6 */ "%3d%1s%1d" "%3d%1s%1d" "%3d%1s%1d" "%2d%1s%1d" /* +12=79 */
        " " /* to ensure that we move to the start of the next line, else end of file */

    /* 1 */
        ,&year, &month, &day, &hour
        ,&Iegh, &Iedn
        ,READ(Igh) /* I values in Wh/m2 */
        ,READ(Idn)
        ,READ(Idh)
    /* 2 */
        ,READ(Lgh) /* L values in kCd/m2 */
        ,READ(Ldn)
        ,READ(Ldh)
        ,READ(Lz)
    /* 3 */
        ,READ(covtot)
        ,READ(covopq)
        ,READ(T)
        ,READ(Tdew)
        ,READ(rh)
        ,READ(p)
    /* 4 */
        ,READ(wdir)
        ,READ(wvel)
        ,READ(vis)
        ,READ(ch)
    /* 5 */     
        ,&obs, &storm, &precip, &drizz, &snowtype, &snowshower, &sleet, &fog, &smog, &hail
    /* 6 */
        ,READ(rain)
        ,READ(aer)
        ,READ(snow)
        ,READ(dsno)
    );

    if(res!=79){
        CONSOLE_DEBUG("Bad input data in data row %d (read %d items OK) (%d/%d/%d %2d:00",d->i,res,day,month,year,hour);
        return 1;
    }

    /*
    if(month!=lastmonth || day!=lastday){
        CONSOLE_DEBUG("Reading data for %d/%d",day,month);
        lastmonth=month;
        lastday=day;
    }
    */

    /* 
        for the moment, we only record global horizontal, direct normal,
        ambient temperature, wind speed.
    */

    tmy = &DATA(d);
    tmy->t = ((day_of_year_specific(day,month,year) - 1)*24 + hour)*3600;
    tmy->G = (float)Igh; /* average W/m2 for the hour in question */
    tmy->Gbn = (float)Idn; /* normal beam radiation */
    tmy->T = 0.1*(float)T + 273.15; /* temperature */
    tmy->v_wind = (float)wvel;
    d->i++;

    return 0;
}

int datareader_tmy2_time(DataReader *d, double *t){
    *t = DATA(d).t;
    return 0;
}

int datareader_tmy2_vals(DataReader *d, double *v){
    v[0]=DATA(d).G;
    v[1]=DATA(d).Gbn;
    return 0;
}
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 TMY2 format.
21
22	These functions implement a reader interface for meteorological data in the
23	TMY2 format as specified at http://rredc.nrel.gov/solar/pubs/tmy2/tab3-2.html
24	//
25	by John Pye, Aug 2006
26	*/
27
28	#include <stdio.h>
29	/* #include <libradtran/sun.h> */
30	#include "sun.h"
31
32	#include <utilities/ascMalloc.h>
33	#include <utilities/error.h>
34
35	#include "tmy.h"
36
37
38	typedef struct{
39	double t;
40	float G;
41	float Gbn;
42	float T;
43	float v_wind;
44
45	} Tmy2Point;
46
47	#define DATA(D) ((Tmy2Point *)(D->data))[D->i]
48
49	/**
50	@return 0 on success
51	*/
52	int datareader_tmy2_header(DataReader *d){
53	char wban[-2 + 6 +2];
54	char city[-8 +29 +2];
55	char state[-31+32+2];
56	int zone;
57	char lathemi;
58	int latdeg, latmin;
59	char longhemi;
60	int longdeg, longmin;
61	int elev;
62
63	fscanf(d->f,"%s %s %s %d"
64	" %c %d %d"
65	" %c %d %d"
66	" %d"
67	,wban,city,state,&zone
68	,&lathemi,&latdeg,&latmin
69	,&longhemi,&longdeg,&longmin
70	,&elev
71	);
72
73	double lat = latdeg + latmin/60;
74	if(lathemi=='S')lat=-lat;
75	double lng = longdeg + longmin/60;
76	if(longhemi=='E')lng=-lng;
77	CONSOLE_DEBUG( "TMY2 data for city '%s' (WBAN %s, time zone %+d) at lat=%.3f, long=%.3f, elev=%d m"
78	, city, wban, zone, lat, lng, elev
79	);
80	ERROR_REPORTER_HERE(ASC_PROG_NOTE,"TMY2 data is for %s, %s",city, state);
81
82	d->i = 0;
83	d->ndata=8760;
84	d->data = ASC_NEW_ARRAY(Tmy2Point,d->ndata);
85	return 0;
86	}
87
88	int datareader_tmy2_eof(DataReader *d){
89	if(feof(d->f)){
90	CONSOLE_DEBUG("REACHED END OF FILE");
91	d->ndata=d->i;
92	ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Read %d rows",d->ndata);
93	return 1;
94	}
95
96	/* set the number of inputs and outputs */
97	d->ninputs = 1;
98	d->noutputs = 2;
99	return 0;
100	}
101
102	#define MEAS(N) int N; char N##_source; int N##_uncert
103	#define READ(N) &N, &N##_source, &N##_uncert
104
105	/**
106	Read a line of data and store in d.
107	@return 0 on success
108	*/
109	int datareader_tmy2_data(DataReader *d){
110	/* static int lastmonth=-1;
111	static int lastday=-1; */
112	int res = 0;
113
114	Tmy2Point *tmy;
115	int year,month,day,hour;
116	int Iegh,Iedn; // Irradiation
117	MEAS(Igh); /* Global horizontal irradiation in the specified interval / (Wh/m2) */
118	MEAS(Idn); /* Direct normal irradiation in the specified interval / (Wh/m2) */
119	MEAS(Idh); /* Diffuse horizontal irradiation in the specified interval / (Wh/m2) */
120
121	MEAS(Lgh); /* Global horiz illuminance / (100 lux) */
122	MEAS(Ldn); /* Direct normal illuminance / (100 lux) */
123	MEAS(Ldh); /* Diffuse horiz illuminance / (100 lux) */
124	MEAS(Lz); /* Zenith illuminance / (10 Cd/m2) */
125
126	MEAS(covtot); /* Total sky cover / tenths */
127	MEAS(covopq); /* Opaque sky cover / tenths */
128
129	MEAS(T); /* temperature / (0.1degC) */
130	MEAS(Tdew); /* dew point temperature / (0.1degC) */
131	MEAS(p); /* pressure / mbar */
132	MEAS(rh); /* rel humidity / % */
133	MEAS(wdir); /* wind dir, N=0, E=90,... */
134	MEAS(wvel); /* wind speed / (m/s) */
135	MEAS(vis); /* visibility / (100m) */
136	MEAS(ch); /* ceiling height / m (or special value) */
137
138	MEAS(rain); /* preciptable water / mm */
139	MEAS(aer); /* aerosol optical depth in thousandths (???) */
140	MEAS(snow); /* snow depth on the specified day / cm (999=missing data) */
141	MEAS(dsno); /* days since last snow (or special value) */
142
143	/* weather observations from Appendix B */
144	int obs, storm, precip, drizz, snowtype, snowshower, sleet, fog, smog, hail;
145
146	/* brace yourself for this one... */
147
148	res = fscanf(d->f,
149	/* 1 / "%2d%2d%2d%2d" "%4d%4d" "%4d%1s%1d" "%4d%1s%1d" "%4d%1s%1d" / =15 */
150	/* 2 / "%4d%1s%1d" "%4d%1s%1d" "%4d%1s%1d" "%4d%1s%1d" / +12=27 */
151	/* 3 / "%2d%1s%1d" "%2d%1s%1d" "%4d%1s%1d" "%4d%1s%1d" "%3d%1s%1d" "%4d%1s%1d" / +18=45 */
152	/* 4 / "%3d%1s%1d" "%3d%1s%1d" "%4d%1s%1d" "%5d%1s%1d" / +12=57 */
153	/* 5 / "%1d%1d%1d%1d%1d%1d%1d%1d%1d%1d" / +10=67 */
154	/* 6 / "%3d%1s%1d" "%3d%1s%1d" "%3d%1s%1d" "%2d%1s%1d" / +12=79 */
155	" " /* to ensure that we move to the start of the next line, else end of file */
156
157	/* 1 */
158	,&year, &month, &day, &hour
159	,&Iegh, &Iedn
160	,READ(Igh) /* I values in Wh/m2 */
161	,READ(Idn)
162	,READ(Idh)
163	/* 2 */
164	,READ(Lgh) /* L values in kCd/m2 */
165	,READ(Ldn)
166	,READ(Ldh)
167	,READ(Lz)
168	/* 3 */
169	,READ(covtot)
170	,READ(covopq)
171	,READ(T)
172	,READ(Tdew)
173	,READ(rh)
174	,READ(p)
175	/* 4 */
176	,READ(wdir)
177	,READ(wvel)
178	,READ(vis)
179	,READ(ch)
180	/* 5 */
181	,&obs, &storm, &precip, &drizz, &snowtype, &snowshower, &sleet, &fog, &smog, &hail
182	/* 6 */
183	,READ(rain)
184	,READ(aer)
185	,READ(snow)
186	,READ(dsno)
187	);
188
189	if(res!=79){
190	CONSOLE_DEBUG("Bad input data in data row %d (read %d items OK) (%d/%d/%d %2d:00",d->i,res,day,month,year,hour);
191	return 1;
192	}
193
194	/*
195	if(month!=lastmonth \|\| day!=lastday){
196	CONSOLE_DEBUG("Reading data for %d/%d",day,month);
197	lastmonth=month;
198	lastday=day;
199	}
200	*/
201
202	/*
203	for the moment, we only record global horizontal, direct normal,
204	ambient temperature, wind speed.
205	*/
206
207	tmy = &DATA(d);
208	tmy->t = ((day_of_year_specific(day,month,year) - 1)24 + hour)3600;
209	tmy->G = (float)Igh; /* average W/m2 for the hour in question */
210	tmy->Gbn = (float)Idn; /* normal beam radiation */
211	tmy->T = 0.1(float)T + 273.15; / temperature */
212	tmy->v_wind = (float)wvel;
213	d->i++;
214
215	return 0;
216	}
217
218	int datareader_tmy2_time(DataReader d, double t){
219	*t = DATA(d).t;
220	return 0;
221	}
222
223	int datareader_tmy2_vals(DataReader d, double v){
224	v[0]=DATA(d).G;
225	v[1]=DATA(d).Gbn;
226	return 0;
227	}