Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: MSG Native Reader
4 : * Purpose: Basic types implementation.
5 : * Author: Frans van den Bergh, fvdbergh@csir.co.za
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2005, Frans van den Bergh <fvdbergh@csir.co.za>
9 : *
10 : * Permission is hereby granted, free of charge, to any person obtaining a
11 : * copy of this software and associated documentation files (the "Software"),
12 : * to deal in the Software without restriction, including without limitation
13 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
14 : * and/or sell copies of the Software, and to permit persons to whom the
15 : * Software is furnished to do so, subject to the following conditions:
16 : *
17 : * The above copyright notice and this permission notice shall be included
18 : * in all copies or substantial portions of the Software.
19 : *
20 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
21 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
22 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
23 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
24 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
25 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
26 : * DEALINGS IN THE SOFTWARE.
27 : ****************************************************************************/
28 :
29 : #include "cpl_port.h"
30 : #include "cpl_error.h"
31 : #include "msg_basic_types.h"
32 :
33 : #include <stdio.h>
34 :
35 : namespace msg_native_format
36 : {
37 :
38 : #ifndef SQR
39 : #define SQR(x) ((x) * (x))
40 : #endif
41 :
42 : // endian conversion routines
43 0 : void to_native(GP_PK_HEADER &h)
44 : {
45 0 : h.sourceSUId = CPL_MSBWORD32(h.sourceSUId);
46 0 : h.sequenceCount = CPL_MSBWORD16(h.sequenceCount);
47 0 : h.packetLength = CPL_MSBWORD32(h.packetLength);
48 0 : }
49 :
50 0 : void to_native(GP_PK_SH1 &h)
51 : {
52 0 : h.spacecraftId = CPL_MSBWORD16(h.spacecraftId);
53 0 : h.packetTime.day = CPL_MSBWORD16(h.packetTime.day);
54 0 : h.packetTime.ms = CPL_MSBWORD32(h.packetTime.ms);
55 0 : }
56 :
57 0 : void to_native(SUB_VISIRLINE &v)
58 : {
59 0 : v.satelliteId = CPL_MSBWORD16(v.satelliteId);
60 0 : v.lineNumberInVisirGrid = CPL_MSBWORD32(v.lineNumberInVisirGrid);
61 0 : }
62 :
63 0 : static void swap_64_bits(unsigned char *b)
64 : {
65 0 : CPL_SWAP64PTR(b);
66 0 : }
67 :
68 0 : void to_native(RADIOMETRIC_PROCESSING_RECORD &r)
69 : {
70 0 : for (int i = 0; i < 12; i++)
71 : {
72 0 : swap_64_bits((unsigned char *)&r.level1_5ImageCalibration[i].cal_slope);
73 0 : swap_64_bits(
74 0 : (unsigned char *)&r.level1_5ImageCalibration[i].cal_offset);
75 : }
76 0 : }
77 :
78 0 : static void to_native(REFERENCEGRID_VISIR &r)
79 : {
80 0 : r.numberOfLines = CPL_MSBWORD32(r.numberOfLines);
81 0 : r.numberOfColumns = CPL_MSBWORD32(r.numberOfColumns);
82 : // should floats be swapped too?
83 : float f;
84 :
85 : // convert float using CPL_MSBPTR32
86 0 : memcpy(&f, &r.lineDirGridStep, sizeof(f));
87 0 : CPL_MSBPTR32(&f);
88 0 : r.lineDirGridStep = f;
89 :
90 : // convert float using CPL_MSBPTR32
91 0 : memcpy(&f, &r.columnDirGridStep, sizeof(f));
92 0 : CPL_MSBPTR32(&f);
93 0 : r.columnDirGridStep = f;
94 0 : }
95 :
96 0 : static void to_native(PLANNED_COVERAGE_VISIR &r)
97 : {
98 0 : r.southernLinePlanned = CPL_MSBWORD32(r.southernLinePlanned);
99 0 : r.northernLinePlanned = CPL_MSBWORD32(r.northernLinePlanned);
100 0 : r.easternColumnPlanned = CPL_MSBWORD32(r.easternColumnPlanned);
101 0 : r.westernColumnPlanned = CPL_MSBWORD32(r.westernColumnPlanned);
102 0 : }
103 :
104 0 : static void to_native(PLANNED_COVERAGE_HRV &r)
105 : {
106 0 : r.lowerSouthLinePlanned = CPL_MSBWORD32(r.lowerSouthLinePlanned);
107 0 : r.lowerNorthLinePlanned = CPL_MSBWORD32(r.lowerNorthLinePlanned);
108 0 : r.lowerEastColumnPlanned = CPL_MSBWORD32(r.lowerEastColumnPlanned);
109 0 : r.lowerWestColumnPlanned = CPL_MSBWORD32(r.lowerWestColumnPlanned);
110 0 : r.upperSouthLinePlanned = CPL_MSBWORD32(r.upperSouthLinePlanned);
111 0 : r.upperNorthLinePlanned = CPL_MSBWORD32(r.upperNorthLinePlanned);
112 0 : r.upperEastColumnPlanned = CPL_MSBWORD32(r.upperEastColumnPlanned);
113 0 : r.upperWestColumnPlanned = CPL_MSBWORD32(r.upperWestColumnPlanned);
114 0 : }
115 :
116 0 : void to_native(IMAGE_DESCRIPTION_RECORD &r)
117 : {
118 : float f;
119 :
120 : // convert float using CPL_MSBPTR32
121 0 : memcpy(&f, &r.longitudeOfSSP, sizeof(f));
122 0 : CPL_MSBPTR32(&f);
123 0 : r.longitudeOfSSP = f;
124 :
125 0 : to_native(r.referencegrid_visir);
126 0 : to_native(r.referencegrid_hrv);
127 0 : to_native(r.plannedCoverage_visir);
128 0 : to_native(r.plannedCoverage_hrv);
129 0 : }
130 :
131 0 : void to_native(ACTUAL_L15_COVERAGE_VISIR_RECORD &r)
132 : {
133 0 : r.southernLineActual = CPL_MSBWORD32(r.southernLineActual);
134 0 : r.northernLineActual = CPL_MSBWORD32(r.northernLineActual);
135 0 : r.easternColumnActual = CPL_MSBWORD32(r.easternColumnActual);
136 0 : r.westernColumnActual = CPL_MSBWORD32(r.westernColumnActual);
137 0 : }
138 :
139 0 : void to_native(ACTUAL_L15_COVERAGE_HRV_RECORD &r)
140 : {
141 0 : r.lowerSouthLineActual = CPL_MSBWORD32(r.lowerSouthLineActual);
142 0 : r.lowerNorthLineActual = CPL_MSBWORD32(r.lowerNorthLineActual);
143 0 : r.lowerEastColumnActual = CPL_MSBWORD32(r.lowerEastColumnActual);
144 0 : r.lowerWestColumnActual = CPL_MSBWORD32(r.lowerWestColumnActual);
145 0 : r.upperSouthLineActual = CPL_MSBWORD32(r.upperSouthLineActual);
146 0 : r.upperNorthLineActual = CPL_MSBWORD32(r.upperNorthLineActual);
147 0 : r.upperEastColumnActual = CPL_MSBWORD32(r.upperEastColumnActual);
148 0 : r.upperWestColumnActual = CPL_MSBWORD32(r.upperWestColumnActual);
149 0 : }
150 :
151 0 : void to_string(PH_DATA &d)
152 : {
153 0 : d.name[29] = 0;
154 0 : d.value[49] = 0;
155 0 : }
156 :
157 : #ifdef notdef
158 : // unit tests on structures
159 : bool perform_type_size_check(void)
160 : {
161 : bool success = true;
162 : if (sizeof(MAIN_PROD_HEADER) != 3674)
163 : {
164 : fprintf(stderr, "MAIN_PROD_HEADER size not 3674 (%lu)\n", /*ok*/
165 : (unsigned long)sizeof(MAIN_PROD_HEADER));
166 : success = false;
167 : }
168 : if (sizeof(SECONDARY_PROD_HEADER) != 1120)
169 : {
170 : fprintf(stderr, "SECONDARY_PROD_HEADER size not 1120 (%lu)\n", /*ok*/
171 : (unsigned long)sizeof(SECONDARY_PROD_HEADER));
172 : success = false;
173 : }
174 : if (sizeof(SUB_VISIRLINE) != 27)
175 : {
176 : fprintf(stderr, "SUB_VISIRLINE size not 17 (%lu)\n", /*ok*/
177 : (unsigned long)sizeof(SUB_VISIRLINE));
178 : success = false;
179 : }
180 : if (sizeof(GP_PK_HEADER) != 22)
181 : {
182 : fprintf(stderr, "GP_PK_HEADER size not 22 (%lu)\n", /*ok*/
183 : (unsigned long)sizeof(GP_PK_HEADER));
184 : success = false;
185 : }
186 : if (sizeof(GP_PK_SH1) != 16)
187 : {
188 : fprintf(stderr, "GP_PK_SH1 size not 16 (%lu)\n", /*ok*/
189 : (unsigned long)sizeof(GP_PK_SH1));
190 : success = false;
191 : }
192 : return success;
193 : }
194 : #endif
195 :
196 : const double Conversions::altitude = 42164; // km from origin
197 : // the spheroid in CGMS 03 4.4.3.2 is unique - flattening is 1/295.488
198 : // note the req and rpol were revised in issue 2.8 of CGMS/DOC/12/0017 - these
199 : // are the revised values
200 : const double Conversions::req = 6378.1370; // earth equatorial radius
201 : const double Conversions::rpol = 6356.7523; // earth polar radius
202 :
203 : // replace the magic constants in the paper with the calculated values in case
204 : // of further change
205 : const double Conversions::dtp2 =
206 : (SQR(altitude) -
207 : SQR(req)); // square of the distance to the equatorial tangent point
208 : // first/last point sensed on the equator
209 :
210 : // given req and rpol, oblate is already defined. Unused afaik in the gdal code
211 : const double Conversions::oblate = ((req - rpol) / req); // oblateness of earth
212 : const double Conversions::eccentricity2 =
213 : (1.0 - (SQR(rpol) / SQR(req))); // 0.00669438...
214 : const double Conversions::ratio2 =
215 : SQR(rpol / req); // 0.9933056 1/x = 1.006739501
216 : const double Conversions::deg_to_rad = (M_PI / 180.0);
217 : const double Conversions::rad_to_deg = (180.0 / M_PI);
218 : const double Conversions::nlines = 3712; // number of lines in an image
219 : const double Conversions::step =
220 : 17.83 / nlines; // pixel / line step in degrees
221 :
222 : const int Conversions::CFAC = -781648343;
223 : const int Conversions::LFAC = -781648343;
224 : const int Conversions::COFF = 1856;
225 : const int Conversions::LOFF = 1856;
226 : const double Conversions::CFAC_scaled = ((double)CFAC / (1 << 16));
227 : const double Conversions::LFAC_scaled = ((double)LFAC / (1 << 16));
228 :
229 : #define SQR(x) ((x) * (x))
230 :
231 0 : void Conversions::convert_pixel_to_geo(double line, double column,
232 : double &longitude, double &latitude)
233 : {
234 : // x and y are angles in radians
235 0 : double x = (column - COFF - 0.0) / CFAC_scaled;
236 0 : double y = (line - LOFF - 0.0) / LFAC_scaled;
237 :
238 0 : double sd = sqrt(SQR(altitude * cos(x) * cos(y)) -
239 0 : (SQR(cos(y)) + SQR(sin(y)) / ratio2) * dtp2);
240 0 : double sn = (altitude * cos(x) * cos(y) - sd) /
241 0 : (SQR(cos(y)) + SQR(sin(y)) / ratio2);
242 0 : double s1 = altitude - sn * SQR(cos(y));
243 0 : double s2 = sn * sin(x) * cos(y);
244 0 : double s3 = -sn * sin(y);
245 0 : double sxy = sqrt(s1 * s1 + s2 * s2);
246 :
247 0 : longitude = atan(s2 / s1);
248 0 : latitude = atan((s3 / sxy) / ratio2);
249 :
250 0 : longitude = longitude * rad_to_deg;
251 0 : latitude = latitude * rad_to_deg;
252 0 : }
253 :
254 0 : void Conversions::compute_pixel_xyz(double line, double column, double &x,
255 : double &y, double &z)
256 : {
257 0 : double asamp = -(column - (nlines / 2.0 + 0.5)) * step;
258 0 : double aline = (line - (nlines / 2.0 + 0.5)) * step;
259 :
260 0 : asamp *= deg_to_rad;
261 0 : aline *= deg_to_rad;
262 :
263 0 : double tanal = tan(aline);
264 0 : double tanas = tan(asamp);
265 :
266 0 : double p = -1;
267 0 : double q = tanas;
268 0 : double r = tanal * sqrt(1 + q * q);
269 :
270 0 : double a = q * q + SQR((r * req / rpol)) + p * p;
271 0 : double b = 2 * altitude * p;
272 0 : double c = altitude * altitude - req * req;
273 :
274 0 : double det = b * b - 4 * a * c;
275 :
276 0 : if (det > 0)
277 : {
278 0 : double k = (-b - sqrt(det)) / (2 * a);
279 0 : x = altitude + k * p;
280 0 : y = k * q;
281 0 : z = k * r;
282 : }
283 : else
284 : {
285 0 : x = y = z = 0;
286 0 : CPLError(CE_Warning, CPLE_AppDefined, "Warning: pixel not visible");
287 : }
288 0 : }
289 :
290 0 : double Conversions::compute_pixel_area_sqkm(double line, double column)
291 : {
292 : double x1, x2;
293 : double y1, y2;
294 : double z1, z2;
295 :
296 0 : compute_pixel_xyz(line - 0.5, column - 0.5, x1, y1, z1);
297 0 : compute_pixel_xyz(line + 0.5, column - 0.5, x2, y2, z2);
298 :
299 0 : double xlen = sqrt(SQR(x1 - x2) + SQR(y1 - y2) + SQR(z1 - z2));
300 :
301 0 : compute_pixel_xyz(line - 0.5, column + 0.5, x2, y2, z2);
302 :
303 0 : double ylen = sqrt(SQR(x1 - x2) + SQR(y1 - y2) + SQR(z1 - z2));
304 :
305 0 : return xlen * ylen;
306 : }
307 :
308 0 : void Conversions::convert_geo_to_pixel(double longitude, double latitude,
309 : unsigned int &line, unsigned int &column)
310 : {
311 :
312 0 : latitude = latitude * deg_to_rad;
313 0 : longitude = longitude * deg_to_rad;
314 :
315 0 : double c_lat = atan(ratio2 * tan(latitude));
316 0 : double r_l = rpol / sqrt(1 - eccentricity2 * cos(c_lat) * cos(c_lat));
317 0 : double r1 = altitude - r_l * cos(c_lat) * cos(longitude);
318 0 : double r2 = -r_l * cos(c_lat) * sin(longitude);
319 0 : double r3 = r_l * sin(c_lat);
320 0 : double rn = sqrt(r1 * r1 + r2 * r2 + r3 * r3);
321 :
322 0 : double x = atan(-r2 / r1) * CFAC_scaled + COFF;
323 0 : double y = asin(-r3 / rn) * LFAC_scaled + LOFF;
324 :
325 0 : line = (unsigned int)floor(x + 0.5);
326 0 : column = (unsigned int)floor(y + 0.5);
327 0 : }
328 :
329 : } // namespace msg_native_format
|
