1 : /******************************************************************************
2 : * $Id: msg_basic_types.cpp 17686 2009-09-25 13:35:11Z dron $
3 : *
4 : * Project: MSG Native Reader
5 : * Purpose: Basic types implementation.
6 : * Author: Frans van den Bergh, fvdbergh@csir.co.za
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2005, Frans van den Bergh <fvdbergh@csir.co.za>
10 : *
11 : * Permission is hereby granted, free of charge, to any person obtaining a
12 : * copy of this software and associated documentation files (the "Software"),
13 : * to deal in the Software without restriction, including without limitation
14 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
15 : * and/or sell copies of the Software, and to permit persons to whom the
16 : * Software is furnished to do so, subject to the following conditions:
17 : *
18 : * The above copyright notice and this permission notice shall be included
19 : * in all copies or substantial portions of the Software.
20 : *
21 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
22 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
23 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
24 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
25 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
26 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
27 : * DEALINGS IN THE SOFTWARE.
28 : ****************************************************************************/
29 :
30 : #include "msg_basic_types.h"
31 : #include "cpl_port.h"
32 :
33 : CPL_CVSID("$Id: msg_basic_types.cpp 17686 2009-09-25 13:35:11Z dron $");
34 :
35 : #include <stdio.h>
36 :
37 : namespace msg_native_format {
38 :
39 : #ifndef SQR
40 : #define SQR(x) ((x)*(x))
41 : #endif
42 :
43 : // endian conversion routines
44 0 : void to_native(GP_PK_HEADER& h) {
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 0 : h.spacecraftId = CPL_MSBWORD16(h.spacecraftId);
52 0 : }
53 :
54 0 : void to_native(SUB_VISIRLINE& v) {
55 0 : v.satelliteId = CPL_MSBWORD16(v.satelliteId);
56 0 : v.lineNumberInVisirGrid = CPL_MSBWORD32(v.lineNumberInVisirGrid);
57 0 : }
58 :
59 0 : static void swap_64_bits(unsigned char* b) {
60 0 : for (int i=0; i < 4; i++) {
61 0 : unsigned char t = b[i];
62 0 : b[i] = b[7-i];
63 0 : b[7-i] = t;
64 : }
65 0 : }
66 :
67 0 : void to_native(RADIOMETRIC_PROCCESSING_RECORD& r) {
68 0 : for (int i=0; i < 12; i++) {
69 0 : swap_64_bits((unsigned char*)&r.level1_5ImageCalibration[i].cal_slope);
70 0 : swap_64_bits((unsigned char*)&r.level1_5ImageCalibration[i].cal_offset);
71 : }
72 0 : }
73 :
74 0 : void to_native(IMAGE_DESCRIPTION_RECORD& r) {
75 0 : r.referencegrid_visir.numberOfLines = CPL_MSBWORD32(r.referencegrid_visir.numberOfLines);
76 0 : r.referencegrid_visir.numberOfColumns = CPL_MSBWORD32(r.referencegrid_visir.numberOfColumns);
77 : // should floats be swapped too?
78 : unsigned int t;
79 :
80 : // convert float using CPL_MSBWORD32
81 0 : t = *(unsigned int *)&r.referencegrid_visir.lineDirGridStep;
82 0 : t = CPL_MSBWORD32(t);
83 0 : r.referencegrid_visir.lineDirGridStep = *(float *)&t;
84 :
85 : // convert float using CPL_MSBWORD32
86 0 : t = *(unsigned int *)&r.referencegrid_visir.columnDirGridStep;
87 0 : t = CPL_MSBWORD32(t);
88 0 : r.referencegrid_visir.columnDirGridStep = *(float *)&t;
89 0 : }
90 :
91 0 : void to_string(PH_DATA& d) {
92 0 : d.name[29] = 0;
93 0 : d.value[49] = 0;
94 0 : }
95 :
96 : // unit tests on structures
97 0 : bool perform_type_size_check(void) {
98 0 : bool success = true;
99 : if (sizeof(MAIN_PROD_HEADER) != 3674) {
100 : fprintf(stderr, "MAIN_PROD_HEADER size not 3674 (%lu)\n", (unsigned long)sizeof(MAIN_PROD_HEADER));
101 : success = false;
102 : }
103 : if (sizeof(SECONDARY_PROD_HEADER) != 1120) {
104 : fprintf(stderr, "SECONDARY_PROD_HEADER size not 1120 (%lu)\n", (unsigned long)sizeof(SECONDARY_PROD_HEADER));
105 : success = false;
106 : }
107 : if (sizeof(SUB_VISIRLINE) != 27) {
108 : fprintf(stderr, "SUB_VISIRLINE size not 17 (%lu)\n", (unsigned long)sizeof(SUB_VISIRLINE));
109 : success = false;
110 : }
111 : if (sizeof(GP_PK_HEADER) != 22) {
112 : fprintf(stderr, "GP_PK_HEADER size not 22 (%lu)\n", (unsigned long)sizeof(GP_PK_HEADER));
113 : success = false;
114 : }
115 : if (sizeof(GP_PK_SH1) != 16) {
116 : fprintf(stderr, "GP_PK_SH1 size not 16 (%lu)\n", (unsigned long)sizeof(GP_PK_SH1));
117 : success = false;
118 : }
119 0 : return success;
120 : }
121 :
122 : const double Conversions::altitude = 42164; // from origin
123 : const double Conversions::req = 6378.1690; // earthequatorial radius
124 : const double Conversions::rpol = 6356.5838; // earth polar radius
125 : const double Conversions::oblate = 1.0/298.257; // oblateness of earth
126 : const double Conversions::deg_to_rad = M_PI/180.0;
127 : const double Conversions::rad_to_deg = 180.0/M_PI;
128 : const double Conversions::nlines = 3712; // number of lines in an image
129 : const double Conversions::step = 17.83/nlines; // pixel / line step in degrees
130 :
131 : const int Conversions::CFAC = -781648343;
132 : const int Conversions::LFAC = -781648343;
133 : const int Conversions::COFF = 1856;
134 : const int Conversions::LOFF = 1856;
135 :
136 : #define SQR(x) ((x)*(x))
137 :
138 0 : void Conversions::convert_pixel_to_geo(double line, double column, double&longitude, double& latitude) {
139 0 : double x = (column - COFF - 0.0) / double(CFAC >> 16);
140 0 : double y = (line - LOFF - 0.0) / double(LFAC >> 16);
141 :
142 0 : double sd = sqrt(SQR(altitude*cos(x)*cos(y)) - (SQR(cos(y)) + 1.006803*SQR(sin(y)))*1737121856);
143 0 : double sn = (altitude*cos(x)*cos(y) - sd)/(SQR(cos(y)) + 1.006803*SQR(sin(y)));
144 0 : double s1 = altitude - sn*cos(x)*cos(y);
145 0 : double s2 = sn*sin(x)*cos(y);
146 0 : double s3 = -sn*sin(y);
147 0 : double sxy = sqrt(s1*s1 + s2*s2);
148 :
149 0 : longitude = atan(s2/s1);
150 0 : latitude = atan(1.006803*s3/sxy);
151 :
152 0 : longitude = longitude / M_PI * 180.0;
153 0 : latitude = latitude / M_PI * 180.0;
154 0 : }
155 :
156 0 : void Conversions::compute_pixel_xyz(double line, double column, double& x,double& y, double& z) {
157 0 : double asamp = -(column - (nlines/2.0 + 0.5)) * step;
158 0 : double aline = (line - (nlines/2.0 + 0.5)) * step;
159 :
160 0 : asamp *= deg_to_rad;
161 0 : aline *= deg_to_rad;
162 :
163 0 : double tanal = tan(aline);
164 0 : double tanas = tan(asamp);
165 :
166 0 : double p = -1;
167 0 : double q = tanas;
168 0 : double r = tanal * sqrt(1 + q*q);
169 :
170 0 : double a = q*q + (r*req/rpol)*(r*req/rpol) + p*p;
171 0 : double b = 2 * altitude * p;
172 0 : double c = altitude * altitude - req*req;
173 :
174 0 : double det = b*b - 4*a*c;
175 :
176 0 : if (det > 0) {
177 0 : double k = (-b - sqrt(det))/(2*a);
178 0 : x = altitude + k*p;
179 0 : y = k * q;
180 0 : z = k * r;
181 :
182 : } else {
183 0 : fprintf(stderr, "Warning: pixel not visible\n");
184 : }
185 0 : }
186 :
187 0 : double Conversions::compute_pixel_area_sqkm(double line, double column) {
188 : double x1, x2;
189 : double y1, y2;
190 : double z1, z2;
191 :
192 0 : compute_pixel_xyz(line-0.5, column-0.5, x1, y1, z1);
193 0 : compute_pixel_xyz(line+0.5, column-0.5, x2, y2, z2);
194 :
195 0 : double xlen = sqrt(SQR(x1 - x2) + SQR(y1 - y2) + SQR(z1 - z2));
196 :
197 0 : compute_pixel_xyz(line-0.5, column+0.5, x2, y2, z2);
198 :
199 0 : double ylen = sqrt(SQR(x1 - x2) + SQR(y1 - y2) + SQR(z1 - z2));
200 :
201 0 : return xlen*ylen;
202 : }
203 :
204 0 : void Conversions::convert_geo_to_pixel(double longitude, double latitude,unsigned int& line, unsigned int& column) {
205 :
206 0 : latitude = latitude / 180.0 * M_PI;
207 0 : longitude = longitude / 180.8 * M_PI;
208 :
209 0 : double c_lat = atan(0.993243 * tan(latitude));
210 0 : double r_l = rpol / sqrt(1 - 0.00675701*cos(c_lat)*cos(c_lat));
211 0 : double r1 = altitude - r_l*cos(c_lat)*cos(longitude);
212 0 : double r2 = -r_l*cos(c_lat)*sin(longitude);
213 0 : double r3 = r_l*sin(c_lat);
214 0 : double rn = sqrt(r1*r1 + r2*r2 + r3*r3);
215 :
216 0 : double x = atan(-r2/r1) * (CFAC >> 16) + COFF;
217 0 : double y = asin(-r3/rn) * (LFAC >> 16) + LOFF;
218 :
219 0 : line = (unsigned int)floor(x + 0.5);
220 0 : column = (unsigned int)floor(y + 0.5);
221 0 : }
222 :
223 : } // namespace msg_native_format
224 :
225 :
|
