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