Home > atmlab > sensors > atovs > avhrr_gac_read.m

avhrr_gac_read

PURPOSE ^

AVHRR_GAC_READ Read, calibrate, geolocate AVHRR GAC L1B data.

SYNOPSIS ^

function S = avhrr_gac_read(filename, varargin)

DESCRIPTION ^

 AVHRR_GAC_READ Read, calibrate, geolocate AVHRR GAC L1B data.

 Read a AVHRR GAC L1B file, apply calibration and geolocation.  The actual
 reading is done in avhrr_gac_read_raw.  Structures of valid header fields
 can be obtained with avhrr_define_gac_l1b.  By default, it returns
 measurements (y), time, lat, lon, epoch.  All other fields that may be
 requested are returned as-is.

 FORMAT

   S = avhrr_gac_read(filename, header_fields, data_fields)

 IN

       filename    

           String, path to file to be read (AVHRR L1B file).  Required.

       extra_header_fields

           Cell array of strings.  Fields will be copied from the L1B
           header.  Call avhrr_define_gac_l1b for valid fields.
           Optional, defaults to {}.

       extra_data_fields

           Cell array of strings.  Arrays will be build by copying data
           from the L1B line records.  Call avhrr_define_gac_l1b for valid
           fields.
           Optional, defaults to {}.

 OUT

   Structure with at least 'time', 'lat', 'lon', 'y', 'epoch'.  Otherwise
   fields copied from data as described above.

 See also: avhrr_gac_define_l1b, avhrr_gac_read_raw

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

DOWNLOAD ^

avhrr_gac_read.m

SOURCE CODE ^

0001 function S = avhrr_gac_read(filename, varargin)
0002 
0003 % AVHRR_GAC_READ Read, calibrate, geolocate AVHRR GAC L1B data.
0004 %
0005 % Read a AVHRR GAC L1B file, apply calibration and geolocation.  The actual
0006 % reading is done in avhrr_gac_read_raw.  Structures of valid header fields
0007 % can be obtained with avhrr_define_gac_l1b.  By default, it returns
0008 % measurements (y), time, lat, lon, epoch.  All other fields that may be
0009 % requested are returned as-is.
0010 %
0011 % FORMAT
0012 %
0013 %   S = avhrr_gac_read(filename, header_fields, data_fields)
0014 %
0015 % IN
0016 %
0017 %       filename
0018 %
0019 %           String, path to file to be read (AVHRR L1B file).  Required.
0020 %
0021 %       extra_header_fields
0022 %
0023 %           Cell array of strings.  Fields will be copied from the L1B
0024 %           header.  Call avhrr_define_gac_l1b for valid fields.
0025 %           Optional, defaults to {}.
0026 %
0027 %       extra_data_fields
0028 %
0029 %           Cell array of strings.  Arrays will be build by copying data
0030 %           from the L1B line records.  Call avhrr_define_gac_l1b for valid
0031 %           fields.
0032 %           Optional, defaults to {}.
0033 %
0034 % OUT
0035 %
0036 %   Structure with at least 'time', 'lat', 'lon', 'y', 'epoch'.  Otherwise
0037 %   fields copied from data as described above.
0038 %
0039 % See also: avhrr_gac_define_l1b, avhrr_gac_read_raw
0040 
0041 % TODO:
0042 %
0043 % - read arbitrary fields
0044 % - write documentation
0045 % - verify correctness
0046 % - finalise
0047 
0048 % Questions:
0049 %
0050 % - AVHRR 3A/3B, FRAME TELEMETRY avh_telem_id or SCAN LINE INFORMATION
0051 % avh_scnlinbin?
0052 % - Should I care about a TIP parity error?
0053 
0054 % $Id: avhrr_gac_read.m 8345 2013-04-17 18:16:40Z gerrit $
0055 
0056 [extra_header_fields, extra_data_fields] = optargs(varargin, {{}, {}});
0057 
0058 core_header_fields = {'avh_h_satid', 'avh_h_instid', 'avh_h_siteid', 'avh_h_hdrcnt', ...
0059     'avh_h_radtempcnv3b1', 'avh_h_radtempcnv3b2', 'avh_h_radtempcnv3bc', ...
0060     'avh_h_radtempcnv41', 'avh_h_radtempcnv42', 'avh_h_radtempcnv4c', ...
0061     'avh_h_radtempcnv51', 'avh_h_radtempcnv52', 'avh_h_radtempcnv5c', ...
0062     'avh_h_scnlin'};
0063 
0064 core_data_fields = {'avh_calvis_os11', 'avh_calvis_oi11', ...
0065      'avh_calvis_os12', 'avh_calvis_oi12', ...
0066      'avh_calvis_oi1', ...
0067      'avh_calvis_os21', 'avh_calvis_oi21', ...
0068      'avh_calvis_os22', 'avh_calvis_oi22', ...
0069      'avh_calvis_oi2', ...
0070      'avh_calvis_os3a1', 'avh_calvis_oi3a1', ...
0071      'avh_calvis_os3a2', 'avh_calvis_oi3a2', ...
0072      'avh_calvis_oi3a', ...
0073      'avh_calir_o3b1', 'avh_calir_o3b2', 'avh_calir_o3b3', ...
0074      'avh_calir_o41', 'avh_calir_o42', 'avh_calir_o43', ...
0075      'avh_calir_o51', 'avh_calir_o52', 'avh_calir_o53', ...
0076      'avh_video', 'avh_scnlin', 'avh_scnlinyr', ...
0077      'avh_telem_id', ...
0078      'avh_scnlintime', 'avh_pos', 'avh_scnlindy', ...
0079      'avh_scnlinbit', 'avh_calqual', 'avh_qualind', 'avh_scnlinqual_c', 'avh_scnlinqual_t', 'avh_scnlinqual_e', 'avh_navstat'};
0080 
0081  
0082 all_header_fields = vec2row(union(core_header_fields, extra_header_fields));
0083 all_data_fields = vec2row(union(core_data_fields, extra_data_fields));
0084  
0085 %limit = 500; % while still in progress
0086 [data_head, data_line] = avhrr_gac_read_raw(filename, ...
0087     all_header_fields, all_data_fields);
0088 
0089 S.y = permute(calibrate_avhrr(data_head, data_line), [2, 1, 3]);
0090 S.time = compensate_wraparound(single(data_line.avh_scnlintime)/1000).';
0091 [S.lat, S.lon] = navigate_avhrr(data_head, data_line);
0092 S.lat = S.lat.';
0093 S.lon = S.lon.';
0094 doys = dayofyear_inverse(data_line.avh_scnlinyr(1), data_line.avh_scnlindy(1));
0095 S.epoch = date2unixsecs(doys.year, doys.month, doys.day);
0096 
0097 % copy over remaining fields
0098 
0099 for hfield = extra_header_fields
0100     S.(hfield{1}) = data_head.(hfield{1});
0101 end
0102 
0103 for dfield = extra_data_fields
0104     S.(dfield{1}) = data_line.(dfield{1}).';
0105 end
0106 
0107 end
0108 
0109 function y = calibrate_avhrr(data_head, data_line)
0110 
0111 % Calibrate AVHRR GAC solar and thermal channels
0112 %
0113 % FORMAT
0114 %
0115 %   calib = calibrated_avhrr(data_head, data_line)
0116 
0117 N_lines = length(data_line.avh_scnlin);
0118 
0119 % Field avh_telem_id contains bit flags (Scan Line Bit Field).
0120 % The first 2 bits specify "channel 3 select (0=3B; 1=3A; 2=transition)".
0121 
0122 chan3_is_3a_alt = logical(bitand(par(data_line.avh_scnlinbit, 1, ':'), 1));
0123 chan3_is_3a = logical(bitand(data_line.avh_scnlinbit, bitshift(1, 0)));
0124 
0125 if ~isequal(chan3_is_3a, chan3_is_3a_alt)
0126     error(['atmlab:' mfilename ':ambiguity'], ...
0127         'Different filds give different into as to 3A/3B status.  Investigate!');
0128 end
0129 
0130 % look at various flags.  See NOAA KLM User's Guide, Table 8.3.1.4.3.2-1
0131 do_not_use = logical(bitand(data_line.avh_qualind, bitshift(1, 31)));
0132 % FIXME: should I care about a tip_parity_error?
0133 %tip_parity_error = logical(bitand(data_line.avh_qualind, bitshift(1, 8)));
0134 cannot_calibrate = logical(bitand(data_line.avh_qualind, bitshift(1, 28)));
0135 bad_calibration = logical(data_line.avh_calqual); % any bit on means bad
0136 bad_calibration2 = logical(data_line.avh_scnlinqual_c);
0137 bad_geolocation = logical(data_line.avh_scnlinqual_e); % any bit on means bad
0138 bad_time = logical(data_line.avh_scnlinqual_t);
0139 
0140 % Unpack 10-bit data into 16-bits for easier handling
0141 %
0142 % TODO: make this call vectorised; consumes 70% of time for
0143 % calibrate_avhrr
0144 % (but calibrate_avhrr consumes only 17% of avhrr_gac_read)
0145 %ce = par(unpack_bip(data_line(1).avh_video), ':', 4);
0146 X = arrayfun(@(i) unpack_bip(data_line.avh_video(:, i)), 1:N_lines, 'UniformOutput', false);
0147 ce = reshape(vertcat(X{:}), [409, N_lines, 5]);
0148 
0149 % pre-allocate calibrated radiance array
0150 Ne = zeros(size(ce), 'single');
0151 
0152 %% calibrate visible channels
0153 %
0154 % According to NOAA KLM User's Guide, Section 7.1.1, page 7-2 (page 248)
0155 % and onward
0156 % http://www.ncdc.noaa.gov/oa/pod-guide/ncdc/docs/klm/html/c7/sec7-1.htm#a041801aa
0157 %
0158 % A = SC + I
0159 %
0160 % with one pair of (S, I) for C<avh_calvis_oi1 and one pair (S, I) for
0161 % C>avh_calvis_oi1
0162 
0163 intercept = bsxfun(@lt, ...
0164                    single(ce(:, :, 1:3)), ...
0165                    shiftdim([data_line.avh_calvis_oi1; ...
0166                              data_line.avh_calvis_oi2; ...
0167                              data_line.avh_calvis_oi3a].', ...
0168                             -1));
0169 
0170 below = bsxfun(@plus, ...
0171                bsxfun(@times, ...
0172                       shiftdim([data_line.avh_calvis_os11; ...
0173                                 data_line.avh_calvis_os21; ...
0174                                 data_line.avh_calvis_os3a1].', ...
0175                                -1), ...
0176                       single(ce(:, :, 1:3))), ...
0177                shiftdim([data_line.avh_calvis_oi11; ...
0178                          data_line.avh_calvis_oi21; ...
0179                          data_line.avh_calvis_oi3a1].', ...
0180                         -1));
0181 
0182 above = bsxfun(@plus, ...
0183                bsxfun(@times, ...
0184                       shiftdim([data_line.avh_calvis_os12; ...
0185                                 data_line.avh_calvis_os22; ...
0186                                 data_line.avh_calvis_os3a2].', ...
0187                                -1), ...
0188                       single(ce(:, :, 1:3))), ...
0189                shiftdim([data_line.avh_calvis_oi12; ...
0190                          data_line.avh_calvis_oi22; ...
0191                          data_line.avh_calvis_oi3a2].', ...
0192                         -1));
0193            
0194 albedo = zeros(size(below), 'single');
0195 albedo(intercept) = below(intercept);
0196 albedo(~intercept) = above(~intercept);
0197 
0198 %% calibrate thermal channels
0199 %
0200 % According to NOAA KLM User's Guide, Section 7.1.2.3, page 7-7 (page 253),
0201 % http://www.ncdc.noaa.gov/oa/pod-guide/ncdc/docs/klm/html/c7/sec7-1.htm#a091405a
0202 %
0203 % Equation (7.1.2.3-1):
0204 %
0205 % Ne = a0 + a1*ce + a2*ce^2;
0206 
0207 % extract calibration coefficients
0208 
0209 a0 = [data_line.avh_calir_o3b1; ...
0210       data_line.avh_calir_o41; ...
0211       data_line.avh_calir_o51];
0212 a1 = [data_line.avh_calir_o3b2; ...
0213        data_line.avh_calir_o42; ...
0214        data_line.avh_calir_o52];
0215 a2 = [data_line.avh_calir_o3b3; ...
0216        data_line.avh_calir_o43; ...
0217        data_line.avh_calir_o53];
0218 
0219 % calibrate channel 3B where appropriate
0220 
0221 Ne(:, ~chan3_is_3a, 3) = ...
0222     bsxfun(@plus, ...
0223            a0(1, :), ...
0224            bsxfun(@plus, ...
0225                   bsxfun(@times, a1(1, :), single(ce(:, :, 3))), ...
0226                   bsxfun(@times, a2(1, :), single(ce(:, :, 3)).^2)));
0227 
0228 % calibrate channels 4 and 5
0229 
0230 term2 = bsxfun(@times, shiftdim(a1(2:3, :).', -1), single(ce(:, :, 4:5)));
0231 term3 = bsxfun(@times, shiftdim(a2(2:3, :).', -1), single(ce(:, :, 4:5)).^2);
0232 Ne(:, :, 4:5) = bsxfun(@plus, ...
0233                        shiftdim(a0(2:3, :).', -1), ...
0234                        bsxfun(@plus, term2, term3));
0235 
0236 if any(any(any(Ne<0)))                  
0237 %    logtext(atmlab('ERR'), 'Setting %d negative radiances to 0!\n', ...
0238 %        length(find(Ne<0)));
0239     Ne(Ne<0) = 0;
0240 end
0241 % now Ne contains radiance in mW/(m2-sr-cm-1)
0242 
0243 %% Convert radiances to brightness temperatures
0244 %
0245 % According to NOAA KLM User's Guide, equation (7.1.2.4-8, -9).
0246 % Page 7-10 (page 256),
0247 % http://www.ncdc.noaa.gov/oa/pod-guide/ncdc/docs/klm/html/c7/sec7-1.htm#a070102d
0248 %
0249 % Te* = c2*vc/log(1 + c1*v^3/Ne)
0250 % Te = (Te*-A)/B
0251 %
0252 % with (c1, c2) from Nigel Atkinson and (A, B) from headers
0253 
0254 vc = [data_head.avh_h_radtempcnv3bc, data_head.avh_h_radtempcnv4c, data_head.avh_h_radtempcnv5c];
0255 b1 = [data_head.avh_h_radtempcnv3b1, data_head.avh_h_radtempcnv41, data_head.avh_h_radtempcnv51];
0256 b2 = [data_head.avh_h_radtempcnv3b2, data_head.avh_h_radtempcnv42, data_head.avh_h_radtempcnv42];
0257 %[vc_, A_, B_] = avhrr_coef_therm2rad(data_head.avh_h_satid);
0258 
0259 % Coefficients from word-document attached by Nigel Atkinson
0260 c1 = 1.191e-5; % mW m^-2 sr^-1 cm^4
0261 c2 = 1.439; % K cm
0262 
0263 Tep = zeros(size(Ne), 'single');
0264 
0265 Tep(:, ~chan3_is_3a, 3) = ...
0266     c2 * vc(1) ./ log(1 + c1 * vc(1).^3 ./ Ne(:, :, 3));
0267 
0268 Tep(:, :, 4:5) = ...
0269     bsxfun(@rdivide, ...
0270            c2 * reshape(vc(2:3), [1, 1, 2]), ...
0271            log(1 + c1 * bsxfun(@rdivide, ...
0272                               reshape(vc(2:3).^3, [1, 1, 2]), ...
0273                               Ne(:, :, 4:5))));
0274 % The following two equations are equivalent:
0275 % Tb = (Tep - A) / B;
0276 % Tb = b1 + b2 * Tep;
0277 % ...with A, B from KLM User's Guide or b1, b2 from headers
0278 
0279 Tb = zeros(size(Tep), 'single');
0280 
0281 Tb(:, ~chan3_is_3a, 3) = b1(1) + b2(1) * Tep(:, ~chan3_is_3a, 3);
0282 
0283 Tb(:, :, 4:5) = bsxfun(@plus, ...
0284                        reshape(b1(2:3), [1, 1, 2]), ...
0285                        bsxfun(@times, ...
0286                               reshape(b2(2:3), [1, 1, 2]), ...
0287                               Tep(:, :, 4:5)));
0288 
0289 %% Put it all together
0290 
0291 y = zeros(size(Ne), 'single');
0292 
0293 y(:, :, 1:2) = albedo(:, :, 1:2);
0294 
0295 y(:, chan3_is_3a, 3) = albedo(:, chan3_is_3a, 3);
0296 
0297 y(:, ~chan3_is_3a, 3) = Tb(:, ~chan3_is_3a, 3);
0298 
0299 y(:, :, 4:5) = Tb(:, :, 4:5);
0300 
0301 wrong = do_not_use | bad_geolocation | bad_time | bad_calibration2 | any(bad_calibration, 1) | cannot_calibrate;
0302 y(:, wrong, :) = -9999.99;
0303 
0304 end
0305 
0306 function [lat, lon] = navigate_avhrr(~, data_line)
0307 % Calculate lan/lon for each observation
0308 %
0309 % From NOAA KLM User's Guide, Table 8.3.1.4.3.2 (page 8-79 / 363)
0310 %
0311 % lat/lon 5, 13, ..., 405
0312 % for viewpos 1, 2, ..., 409
0313 %
0314 % Distance between anchorpoints is 30–150 km
0315 %
0316 % Implementation inspired by IDL function by Nigel Atkinson
0317 
0318 nlines = size(data_line.avh_pos, 2);
0319 loc = reshape(data_line.avh_pos, [2, 51, nlines]);
0320 lat = zeros(409, nlines, 'single');
0321 lon = zeros(409, nlines, 'single');
0322 
0323 pos_given = 5:8:405;
0324 lat_in = squeeze(loc(1, :, :));
0325 lon_in = squeeze(loc(2, :, :));
0326 
0327 x_in = cosd(lon_in).*cosd(lat_in);
0328 y_in = sind(lon_in).*cosd(lat_in);
0329 z_in = sind(lat_in);
0330 
0331 xf = spline(pos_given, x_in.', 1:409).';
0332 yf = spline(pos_given, y_in.', 1:409).';
0333 zf = spline(pos_given, z_in.', 1:409).';
0334 lon(:, :) = atan2(yf, xf) .* constants('RAD2DEG');
0335 lat(:, :) = atan2(zf, sqrt(xf.^2 + yf.^2)) .* constants('RAD2DEG');
0336 
0337 %{
0338 IDL code
0339 
0340 spots_nav = lindgen(51)*8 + 5
0341 d2r = !pi/180
0342 
0343 nspots = 409
0344 spots_req = lindgen(nspots) + 1
0345 
0346 nlines = n_elements(pos(0,0,*))
0347 lat = fltarr(nspots,nlines)
0348 lon = fltarr(nspots,nlines)
0349 
0350 for iline=0,nlines-1 do begin
0351   lat_in = pos(0,*,iline)*1.0E-4*d2r
0352   lon_in = pos(1,*,iline)*1.0E-4*d2r
0353 
0354   x_in = cos(lon_in)*cos(lat_in);
0355   y_in = sin(lon_in)*cos(lat_in);
0356   z_in = sin(lat_in);
0357 
0358   x = spline(spots_nav,x_in,spots_req)
0359   y = spline(spots_nav,y_in,spots_req)
0360   z = spline(spots_nav,z_in,spots_req)
0361 
0362   lon(*,iline) = atan(y,x)/d2r
0363   lat(*,iline) = atan(z,sqrt(x*x + y*y))/d2r
0364 endfor
0365 
0366 %}
0367 
0368 end
0369 %{
0370 function [vc, A, B] = avhrr_coef_therm2rad(satid)
0371 % get coefficients for thermal channel temperature-to-radiance
0372 %
0373 % satid as per field avh_h_satid 'Spacecraft Identification Code'.
0374 % From NOAA KLM User's Guide, Table 8.3.1.4.2.1-1:
0375 %
0376 % 2=NOAA-L
0377 % 4=NOAA-K
0378 % 6=NOAA-M
0379 % 7=NOAA-N
0380 % 8=NOAA-P
0381 % 11=MetOp-1 ( = MetOp-B)
0382 % 12=MetOp-A
0383 %
0384 % Not in use, could be used as backup if header unavailable?
0385 
0386 % NB! MetOp-A = MetOp-2; MetOp-B = MetOp-1!
0387 % http://database.eohandbook.com/database/missionsummary.aspx?missionID=261
0388 % NOAA KLM User's Guide Table D-2 (page 906)
0389 
0390 switch satid
0391     case 2 % NOAA-16
0392         % Table D.2-12 (page D-91 / page 995)
0393         vc = [2700.1148, 917.2289, 838.1255];
0394         A = [1.592459, 0.332380, 0.674623];
0395         B = [0.998147, 0.998522, 0.998363];
0396     case 4 % NOAA-15
0397         % Table D.1-11 (page D-19 / page 923)
0398         vc = [2695.9743, 925.4075, 839.8979];
0399         A = [1.621256, 0.337810, 0.304558];
0400         B = [0.998015, 0.998719, 0.999024];
0401     case 6 % NOAA-17
0402         % Table D.3-7 (page D-117 / page 1021)
0403         vc = [2669.3554, 926.2947, 839.8246];
0404         A = [1.702380, 0.271682, 0.309180];
0405         B = [0.997378, 0.998794, 0.999012];
0406     case 7 % NOAA-18
0407         % Table D.4-7 (page D-183 / page 1087)
0408         vc = [2659.7952, 928.1460, 833.232];
0409         A = [1.698704, 0.436645, 0.253179];
0410         B = [0.996960, 0.998607, 0.999057];
0411     case 8 % NOAA-19
0412         % Table D.6-7 (page D-785 / page 1689)
0413         vc = [2670.0, 928.9, 831.9];
0414         A = [1.67396, 0.53959, 0.36064];
0415         B = [0.997364, 0.998534, 0.998913];
0416 %    case 11 % MetOp-1 ??
0417     case 12 % MetOp-A
0418         % Table D.5-7 (page D-492 / page 1396)
0419         vc = [2687.0, 927.2, 837.7];
0420         A = [2.06699, 0.55126, 0.34716];
0421         B = [0.996577, 0.998509, 0.998947];
0422     otherwise
0423         error(['atmlab:', mfilename, ':unknownid'], ...
0424             'Unknown satellite id: %d', id);
0425 end
0426 end
0427 %}

Generated on Mon 15-Sep-2014 13:31:28 by m2html © 2005