Opened 19 years ago
Last modified 18 years ago
#879 closed enhancement (fixed)
Meteosat Native Archive format (.nat) patch
| Reported by: | Owned by: | warmerdam | |
|---|---|---|---|
| Priority: | high | Milestone: | |
| Component: | default | Version: | unspecified |
| Severity: | minor | Keywords: | |
| Cc: | retsios@… |
Description
diff -Nur gdal/GDALmake.opt.in gdal_patched/GDALmake.opt.in
--- gdal/GDALmake.opt.in 2005-05-19 22:44:02.000000000 +0200
+++ gdal_patched/GDALmake.opt.in 2005-07-01 16:14:22.680937048 +0200
@@ -236,7 +236,7 @@
#
GDAL_FORMATS = gxf gtiff hfa aigrid aaigrid ceos ceos2 iso8211 xpm \
sdts raw dted mem jdem envisat elas fit vrt usgsdem l1b \
- nitf bmp pcidsk airsar rs2 ilwis rmf \
+ nitf bmp pcidsk airsar rs2 ilwis rmf msgn \
@OPT_GDAL_FORMATS@
#
diff -Nur gdal/frmts/gdalallregister.cpp gdal_patched/frmts/gdalallregister.cpp
--- gdal/frmts/gdalallregister.cpp 2005-05-19 22:41:14.000000000 +0200
+++ gdal_patched/frmts/gdalallregister.cpp 2005-07-01 16:07:23.548654840 +0200
@@ -292,6 +292,10 @@
GDALRegister_RMF();
#endif
+#ifdef FRMT_msgn
+ GDALRegister_MSGN();
+#endif
+
/* -------------------------------------------------------------------- */
/* Our test for the following is weak or expensive so we try */
/* them last. */
diff -Nur gdal/frmts/makefile.vc gdal_patched/frmts/makefile.vc
--- gdal/frmts/makefile.vc 2005-05-19 22:41:14.000000000 +0200
+++ gdal_patched/frmts/makefile.vc 2005-07-01 16:12:32.914624064 +0200
@@ -6,7 +6,8 @@
-DFRMT_dted -DFRMT_mem -DFRMT_jdem -DFRMT_gif \
-DFRMT_envisat -DFRMT_aaigrid -DFRMT_usgsdem -DFRMT_l1b \
-DFRMT_fit -DFRMT_vrt -DFRMT_xpm -DFRMT_bmp -DFRMT_rmf \
- -DFRMT_nitf -DFRMT_pcidsk -DFRMT_airsar -DFRMT_rs2 -DFRMT_ilwis
+ -DFRMT_nitf -DFRMT_pcidsk -DFRMT_airsar -DFRMT_rs2 -DFRMT_ilwis \
+ -DFRMT_msgn
MOREEXTRA =
DIRLIST = $(EXTRAFLAGS:-DFRMT_=)
diff -Nur gdal/frmts/msgn/GNUmakefile gdal_patched/frmts/msgn/GNUmakefile
--- gdal/frmts/msgn/GNUmakefile 1970-01-01 02:00:00.000000000 +0200
+++ gdal_patched/frmts/msgn/GNUmakefile 2005-07-01 16:07:23.550654536 +0200
@@ -0,0 +1,17 @@
+GDAL_ROOT = ../..
+
+include $(GDAL_ROOT)/GDALmake.opt
+
+OBJ = msgndataset.o msg_basic_types.o msg_reader_core.o
+
+CPPFLAGS := $(GDAL_INCLUDE) $(CPPFLAGS) -I. -Imodule-msg_reader
-DGDAL_SUPPORT
+
+default: $(OBJ)
+
+clean:
+ rm -f *.o
+
+../o/%.o: module-msg_reader/src/%.cc
+ $(CC) -c $(CPPFLAGS) $(CFLAGS) $< -o $@
+
+install-obj: $(O_OBJ)
diff -Nur gdal/frmts/msgn/frmt_msgn.html gdal_patched/frmts/msgn/frmt_msgn.html
--- gdal/frmts/msgn/frmt_msgn.html 1970-01-01 02:00:00.000000000 +0200
+++ gdal_patched/frmts/msgn/frmt_msgn.html 2005-07-01 16:26:06.216983312 +0200
@@ -0,0 +1,22 @@
+<html>
+<head>
+<title>MSGN -- Meteosat Second Generation (MSG) Native Archive Format
(.nat)</title>
+</head>
+
+<body bgcolor="#ffffff">
+
+<h1>MSGN -- Meteosat Second Generation (MSG) Native Archive Format (.nat)</h1>
+
+GDAL supports reading only of MSG native files. These files may have
+anything from 1 to 12 bands, all at 10-bit resolution.
+
+Currently the driver ignores the 12th band (HRV) because it does not have
+the same resolution as the others, so it can not be accommodated within the
+same dataset. Adding support for this band would be straightforward, though.
+
+Georeferencing is currently supported, but the results may not be
+acceptable, depending on your requirements. The current workaround is to
+implement the CGMS Geostationary projection directly.
+
+</body>
+</html>
diff -Nur gdal/frmts/msgn/makefile.vc gdal_patched/frmts/msgn/makefile.vc
--- gdal/frmts/msgn/makefile.vc 1970-01-01 02:00:00.000000000 +0200
+++ gdal_patched/frmts/msgn/makefile.vc 2005-07-01 16:22:29.534923992 +0200
@@ -0,0 +1,15 @@
+
+OBJ = msgndataset.o msg_basic_types.o msg_reader_core.o
+
+EXTRAFLAGS = -I..\iso8211 -I. -Imodule-msg_reader -DGDAL_SUPPORT
+
+GDAL_ROOT = ..\..
+
+!INCLUDE $(GDAL_ROOT)\nmake.opt
+
+default: $(OBJ)
+ copy *.obj ..\o
+
+clean:
+ -del *.obj
+
diff -Nur gdal/frmts/msgn/module-msg_reader/include/msg_basic_types.h
gdal_patched/frmts/msgn/module-msg_reader/include/msg_basic_types.h
--- gdal/frmts/msgn/module-msg_reader/include/msg_basic_types.h 1970-01-01
02:00:00.000000000 +0200
+++ gdal_patched/frmts/msgn/module-msg_reader/include/msg_basic_types.h
2005-07-01 16:07:23.564652408 +0200
@@ -0,0 +1,218 @@
+/* $Id: msg_basic_types.h 10 2005-03-11 05:41:22Z fvdbergh $
+*/
+
+#ifndef MSG_BASIC_TYPES
+#define MSG_BASIC_TYPES
+
+#include <math.h>
+
+namespace msg_native_format {
+
+const unsigned int SATELLITESTATUS_RECORD_LENGTH = 60134;
+const unsigned int IMAGEACQUISITION_RECORD_LENGTH = 700;
+const unsigned int CELESTIALEVENTS_RECORD_LENGTH = 326058; // should be
56258 according to ICD105 ??
+const unsigned int IMAGEDESCRIPTION_RECORD_LENGTH = 101;
+
+const unsigned int RADIOMETRICPROCESSING_RECORD_OFFSET =
+ SATELLITESTATUS_RECORD_LENGTH +
+ IMAGEACQUISITION_RECORD_LENGTH +
+ CELESTIALEVENTS_RECORD_LENGTH +
+ IMAGEDESCRIPTION_RECORD_LENGTH;
+
+typedef int INTEGER; // 32 bits
+typedef unsigned int UNSIGNED; // 32 bits
+typedef unsigned short USHORT; // 16 bits
+typedef unsigned char TIME_CDS_SHORT[6];
+typedef unsigned char TIME_CDS_EXPANDED[10];
+typedef unsigned char EBYTE; // enumerated byte
+typedef unsigned char UBYTE; // enumerated byte
+typedef float REAL;
+
+typedef unsigned short int GP_SC_ID; // 16 bits, enumerated
+typedef unsigned char GP_SC_CHAN_ID; // 8 bits, enumerated
+typedef unsigned char GP_FAC_ID; // 8 bits, enumerated
+typedef unsigned char GP_FAC_ENV; // 8 bits, enumerated
+typedef unsigned int GP_SU_ID; // 32 bits, interval partition
+typedef unsigned char GP_SVCE_TYPE; // 8 bits, enumerated
+
+// all structures must be packed on byte boundaries
+#pragma pack(1)
+
+typedef struct {
+ unsigned char qualifier1;
+ unsigned char qualifier2;
+ unsigned char qualifier3;
+ unsigned char qualifier4;
+} GP_CPU_ID;
+
+typedef struct {
+ char name[30];
+ char value[50];
+} PH_DATA;
+
+typedef struct {
+ char name[30];
+ char size[16];
+ char address[16];
+} PH_DATA_ID;
+
+typedef struct {
+ PH_DATA formatName;
+ PH_DATA formatDocumentName;
+ PH_DATA formatDocumentMajorVersion;
+ PH_DATA formatDocumentMinorVersion;
+ PH_DATA creationDateTime;
+ PH_DATA creatingCentre;
+ PH_DATA_ID dataSetIdentification[5];
+ UBYTE slack[1364]; // what is this? This is not in the
documentation?
+ PH_DATA totalFileSize;
+ PH_DATA gort;
+ PH_DATA asti;
+ PH_DATA llos;
+ PH_DATA snit;
+ PH_DATA aiid;
+ PH_DATA ssbt;
+ PH_DATA ssst;
+ PH_DATA rrcc;
+ PH_DATA rrbt;
+ PH_DATA rrst;
+ PH_DATA pprc;
+ PH_DATA ppdt;
+ PH_DATA gplv;
+ PH_DATA apnm;
+ PH_DATA aarf;
+ PH_DATA uudt;
+ PH_DATA qqov;
+ PH_DATA udsp;
+} MAIN_PROD_HEADER;
+
+typedef struct {
+ PH_DATA abid;
+ PH_DATA smod;
+ PH_DATA apxs;
+ PH_DATA avpa;
+ PH_DATA lscd;
+ PH_DATA lmap;
+ PH_DATA qdlc;
+ PH_DATA qdlp;
+ PH_DATA qqai;
+ PH_DATA selectedBandIds;
+ PH_DATA southLineSelectedRectangle;
+ PH_DATA northLineSelectedRectangle;
+ PH_DATA eastColumnSelectedRectangle;
+ PH_DATA westColumnSelectedRectangle;
+} SECONDARY_PROD_HEADER;
+
+typedef struct {
+ UBYTE visirlineVersion;
+ GP_SC_ID satelliteId;
+ TIME_CDS_EXPANDED trueRepeatCycleStart;
+ INTEGER lineNumberInVisirGrid;
+ GP_SC_CHAN_ID channelId;
+ TIME_CDS_SHORT l10LineMeanAcquisitionTime;
+ EBYTE lineValidity;
+ EBYTE lineRadiometricQuality;
+ EBYTE lineGeometricQuality;
+ // actual line data not represented here
+} SUB_VISIRLINE;
+
+typedef struct {
+ UBYTE headerVersionNo;
+ EBYTE packetType; // 2 = mission data
+ EBYTE subHeaderType; // 0 = no subheader, 1 = GP_PK_SH1, 2 =
GP_PK_SH2
+ GP_FAC_ID sourceFacilityId;
+ GP_FAC_ENV sourceEnvId;
+ UBYTE sourceInstanceId;
+ GP_SU_ID sourceSUId;
+ GP_CPU_ID sourceCPUId;
+ GP_FAC_ID destFacilityId;
+ GP_FAC_ENV destEnvId;
+ USHORT sequenceCount;
+ UNSIGNED packetLength;
+} GP_PK_HEADER;
+
+typedef struct {
+ UBYTE subHeaderVersionNo;
+ EBYTE checksumFlag;
+ UBYTE acknowledgement[4];
+ GP_SVCE_TYPE serviceType;
+ UBYTE serviceSubType;
+ TIME_CDS_SHORT packetTime;
+ GP_SC_ID spacecraftId;
+} GP_PK_SH1;
+
+typedef struct {
+ double cal_slope;
+ double cal_offset;
+} CALIBRATION;
+
+typedef struct {
+ EBYTE radianceLinearisation[12];
+ EBYTE detectorEqualisation[12];
+ EBYTE onboardCalibrationResult[12];
+ EBYTE MPEFCalFeedback[12];
+ EBYTE MTFAdaption[12];
+ EBYTE straylightCorrectionFlag[12];
+ CALIBRATION level1_5ImageCalibration[12];
+ // rest of structure omitted for now
+} RADIOMETRIC_PROCCESSING_RECORD;
+
+typedef struct {
+ INTEGER numberOfLines;
+ INTEGER numberOfColumns;
+ REAL lineDirGridStep;
+ REAL columnDirGridStep;
+ EBYTE gridOrigin;
+} REFERENCEGRID_VISIR;
+
+typedef struct {
+ EBYTE typeOfProjection;
+ REAL longitudeOfSSP;
+ REFERENCEGRID_VISIR referencegrid_visir;
+ // rest of record omitted, for now
+} IMAGE_DESCRIPTION_RECORD;
+
+// disable byte-packing
+#pragma pack()
+
+// endian conversion routines
+void to_native(GP_PK_HEADER& h);
+void to_native(GP_PK_SH1& h);
+void to_native(SUB_VISIRLINE& v);
+void to_native(RADIOMETRIC_PROCCESSING_RECORD& r);
+void to_native(IMAGE_DESCRIPTION_RECORD& r);
+
+// utility function, alters string fields permanently
+void to_string(PH_DATA& d);
+
+// unit tests on structures, returns true on success
+bool perform_type_size_check(void);
+
+class Conversions {
+public:
+ static void convert_pixel_to_geo(double line, double column, double&
longitude, double& latitude);
+ static void convert_geo_to_pixel(double longitude, double latitude,
unsigned int& line, unsigned int& column);
+
+ static void compute_pixel_xyz(double line, double column, double& x,
double& y, double& z);
+ static double compute_pixel_area_sqkm(double line, double column);
+
+ static const double altitude; // from origin
+ static const double req; // earth equatorial radius
+ static const double rpol; // earth polar radius
+ static const double oblate; // oblateness of earth
+ static const double deg_to_rad;
+ static const double rad_to_deg;
+ static const double step; // pixel / line step in degrees
+ static const double nlines; // number of lines in an image
+
+ static const int CFAC; // Column scale factor
+ static const int LFAC; // Line scale factor
+ static const int COFF; // Column offset
+ static const int LOFF; // Line offset
+
+};
+
+} // msg_native_format
+
+#endif
+
diff -Nur gdal/frmts/msgn/module-msg_reader/include/msg_reader_core.h
gdal_patched/frmts/msgn/module-msg_reader/include/msg_reader_core.h
--- gdal/frmts/msgn/module-msg_reader/include/msg_reader_core.h 1970-01-01
02:00:00.000000000 +0200
+++ gdal_patched/frmts/msgn/module-msg_reader/include/msg_reader_core.h
2005-07-01 16:07:23.564652408 +0200
@@ -0,0 +1,118 @@
+/* $Id: msg_reader_core.h 13 2005-03-14 06:44:22Z fvdbergh $
+ Purpose: Base class for reading in the headers of MSG native images
+*/
+
+#ifndef MSG_READER_CORE_H
+#define MSG_READER_CORE_H
+
+#include "module-msg_reader/include/msg_basic_types.h"
+#include <stdio.h>
+
+namespace msg_native_format {
+
+const unsigned int MSG_NUM_CHANNELS = 12;
+
+typedef struct {
+ double vc;
+ double A;
+ double B;
+} Blackbody_lut_type;
+
+typedef enum {
+ VIS0_6 = 2,
+ VIS0_8 = 4,
+ NIR1_6 = 8,
+ IR3_9 = 16,
+ IR6_2 = 32,
+ IR7_3 = 64,
+ IR8_7 = 128,
+ IR9_7 = 256,
+ IR10_8 = 512,
+ IR12_0 = 1024,
+ IR13_4 = 2048,
+ HRV = 4096
+} Msg_channel_names;
+
+class Msg_reader_core {
+public:
+ Msg_reader_core(const char* fname);
+ Msg_reader_core(FILE* fp);
+ virtual ~Msg_reader_core(void) {};
+
+ #ifndef GDAL_SUPPORT
+ virtual void radiance_to_blackbody(int using_chan_no = 0) = 0; // can
override which channel's parameters to use
+ virtual double* get_data(int chan_no=0) = 0;
+ #endif
+
+ unsigned int get_lines(void) { return _lines; }
+ unsigned int get_columns(void) { return _columns; }
+
+ void get_pixel_geo_coordinates(unsigned int line, unsigned int column,
double& longitude, double& latitude); // x and y relative to this image, not
full disc image
+ void get_pixel_geo_coordinates(double line, double column, double&
longitude, double& latitude); // x and y relative to this image, not full disc image
+ double compute_pixel_area_sqkm(double line, double column);
+
+ static const Blackbody_lut_type Blackbody_LUT[MSG_NUM_CHANNELS+1];
+
+ unsigned int get_year(void) { return _year; }
+ unsigned int get_month(void) { return _month; }
+ unsigned int get_day(void) { return _day; }
+ unsigned int get_hour(void) { return _hour; }
+ unsigned int get_minute(void) { return _minute; }
+
+ unsigned int get_line_start(void) { return _line_start; }
+ unsigned int get_col_start(void) { return _col_start; }
+
+ float get_col_dir_step(void) { return _col_dir_step; }
+ float get_line_dir_step(void) { return _line_dir_step; }
+
+ unsigned int get_f_data_offset(void) { return _f_data_offset; }
+ unsigned int get_visir_bytes_per_line(void) { return _visir_bytes_per_line; }
+ unsigned int get_visir_packet_size(void) { return _visir_packet_size; }
+ unsigned int get_interline_spacing(void) { return _interline_spacing; }
+
+ unsigned char* get_band_map(void) { return _bands; }
+
+ CALIBRATION* get_calibration_parameters(void) { return _calibration; }
+
+private:
+ void read_metadata_block(FILE* fp);
+
+protected:
+
+ int _chan_to_idx(Msg_channel_names channel);
+
+ unsigned int _lines;
+ unsigned int _columns;
+
+ unsigned int _line_start;
+ unsigned int _col_start;
+
+ float _col_dir_step;
+ float _line_dir_step;
+
+ MAIN_PROD_HEADER _main_header;
+ SECONDARY_PROD_HEADER _sec_header;
+ CALIBRATION _calibration[MSG_NUM_CHANNELS];
+
+ unsigned int _f_data_offset;
+ unsigned int _f_data_size;
+ unsigned int _f_header_offset;
+ unsigned int _f_header_size;
+
+ unsigned int _visir_bytes_per_line; // packed length of a VISIR line,
wihtout headers
+ unsigned int _visir_packet_size; // effectively, the spacing between
lines of consecutive bands in bytes
+ unsigned int _hrv_packet_size;
+ unsigned int _interline_spacing;
+
+ unsigned char _bands[MSG_NUM_CHANNELS];
+
+ unsigned int _year;
+ unsigned int _month;
+ unsigned int _day;
+ unsigned int _hour;
+ unsigned int _minute;
+};
+
+}// namespace msg_native_format
+
+#endif
diff -Nur gdal/frmts/msgn/module-msg_reader/src/msg_basic_types.cc
gdal_patched/frmts/msgn/module-msg_reader/src/msg_basic_types.cc
--- gdal/frmts/msgn/module-msg_reader/src/msg_basic_types.cc 1970-01-01
02:00:00.000000000 +0200
+++ gdal_patched/frmts/msgn/module-msg_reader/src/msg_basic_types.cc
2005-07-01 16:28:12.592771272 +0200
@@ -0,0 +1,197 @@
+/* $Id: msg_basic_types.cc 10 2005-03-11 05:41:22Z fvdbergh $
+*/
+
+#include "include/msg_basic_types.h"
+
+#include <netinet/in.h>
+#include <stdio.h>
+
+namespace msg_native_format {
+
+#ifndef SQR
+#define SQR(x) ((x)*(x))
+#endif
+
+// endian conversion routines
+void to_native(GP_PK_HEADER& h) {
+ h.sourceSUId = ntohl(h.sourceSUId);
+ h.sequenceCount = ntohs(h.sequenceCount);
+ h.packetLength = ntohl(h.packetLength);
+}
+
+void to_native(GP_PK_SH1& h) {
+ h.spacecraftId = ntohs(h.spacecraftId);
+}
+
+void to_native(SUB_VISIRLINE& v) {
+ v.satelliteId = ntohs(v.satelliteId);
+ v.lineNumberInVisirGrid = ntohl(v.lineNumberInVisirGrid);
+}
+
+static void swap_64_bits(unsigned char* b) {
+ for (int i=0; i < 4; i++) {
+ unsigned char t = b[i];
+ b[i] = b[7-i];
+ b[7-i] = t;
+ }
+}
+
+void to_native(RADIOMETRIC_PROCCESSING_RECORD& r) {
+ for (int i=0; i < 12; i++) {
+ swap_64_bits((unsigned char*)&r.level1_5ImageCalibration[i].cal_slope);
+ swap_64_bits((unsigned char*)&r.level1_5ImageCalibration[i].cal_offset);
+ }
+}
+
+void to_native(IMAGE_DESCRIPTION_RECORD& r) {
+ r.referencegrid_visir.numberOfLines =
ntohl(r.referencegrid_visir.numberOfLines);
+ r.referencegrid_visir.numberOfColumns =
ntohl(r.referencegrid_visir.numberOfColumns);
+ // should floats be swapped too?
+ unsigned int t;
+
+ // convert float using ntohl
+ t = *(unsigned int *)&r.referencegrid_visir.lineDirGridStep;
+ t = ntohl(t);
+ r.referencegrid_visir.lineDirGridStep = *(float *)&t;
+
+ // convert float using ntohl
+ t = *(unsigned int *)&r.referencegrid_visir.columnDirGridStep;
+ t = ntohl(t);
+ r.referencegrid_visir.columnDirGridStep = *(float *)&t;
+}
+
+void to_string(PH_DATA& d) {
+ d.name[29] = 0;
+ d.value[49] = 0;
+}
+
+// unit tests on structures
+bool perform_type_size_check(void) {
+ bool success = true;
+ if (sizeof(MAIN_PROD_HEADER) != 3674) {
+ fprintf(stderr, "MAIN_PROD_HEADER size not 3674 (%d)\n",
sizeof(MAIN_PROD_HEADER));
+ success = false;
+ }
+ if (sizeof(SECONDARY_PROD_HEADER) != 1120) {
+ fprintf(stderr, "SECONDARY_PROD_HEADER size not 1120 (%d)\n",
sizeof(SECONDARY_PROD_HEADER));
+ success = false;
+ }
+ if (sizeof(SUB_VISIRLINE) != 27) {
+ fprintf(stderr, "SUB_VISIRLINE size not 17 (%d)\n", sizeof(SUB_VISIRLINE));
+ success = false;
+ }
+ if (sizeof(GP_PK_HEADER) != 22) {
+ fprintf(stderr, "GP_PK_HEADER size not 22 (%d)\n", sizeof(GP_PK_HEADER));
+ success = false;
+ }
+ if (sizeof(GP_PK_SH1) != 16) {
+ fprintf(stderr, "GP_PK_SH1 size not 16 (%d)\n", sizeof(GP_PK_SH1));
+ success = false;
+ }
+ return success;
+}
+
+const double Conversions::altitude = 42164; // from origin
+const double Conversions::req = 6378.1690; // earth
equatorial radius
+const double Conversions::rpol = 6356.5838; // earth polar radius
+const double Conversions::oblate = 1.0/298.257; // oblateness of earth
+const double Conversions::deg_to_rad = M_PI/180.0;
+const double Conversions::rad_to_deg = 180.0/M_PI;
+const double Conversions::nlines = 3712; // number of lines
in an image
+const double Conversions::step = 17.83/nlines; // pixel / line
step in degrees
+
+const int Conversions::CFAC = -781648343;
+const int Conversions::LFAC = -781648343;
+const int Conversions::COFF = 1856;
+const int Conversions::LOFF = 1856;
+
+#define SQR(x) ((x)*(x))
+
+void Conversions::convert_pixel_to_geo(double line, double column, double&
longitude, double& latitude) {
+ double x = (column - COFF - 0.0) / double(CFAC >> 16);
+ double y = (line - LOFF - 0.0) / double(LFAC >> 16);
+
+ double sd = sqrt(SQR(altitude*cos(x)*cos(y)) - (SQR(cos(y)) +
1.006803*SQR(sin(y)))*1737121856);
+ double sn = (altitude*cos(x)*cos(y) - sd)/(SQR(cos(y)) + 1.006803*SQR(sin(y)));
+ double s1 = altitude - sn*cos(x)*cos(y);
+ double s2 = sn*sin(x)*cos(y);
+ double s3 = -sn*sin(y);
+ double sxy = sqrt(s1*s1 + s2*s2);
+
+ longitude = atan(s2/s1);
+ latitude = atan(1.006803*s3/sxy);
+
+ longitude = longitude / M_PI * 180.0;
+ latitude = latitude / M_PI * 180.0;
+}
+
+void Conversions::compute_pixel_xyz(double line, double column, double& x,
double& y, double& z) {
+ double asamp = -(column - (nlines/2.0 + 0.5)) * step;
+ double aline = (line - (nlines/2.0 + 0.5)) * step;
+
+ asamp *= deg_to_rad;
+ aline *= deg_to_rad;
+
+ double tanal = tan(aline);
+ double tanas = tan(asamp);
+
+ double p = -1;
+ double q = tanas;
+ double r = tanal * sqrt(1 + q*q);
+
+ double a = q*q + (r*req/rpol)*(r*req/rpol) + p*p;
+ double b = 2 * altitude * p;
+ double c = altitude * altitude - req*req;
+
+ double det = b*b - 4*a*c;
+
+ if (det > 0) {
+ double k = (-b - sqrt(det))/(2*a);
+ x = altitude + k*p;
+ y = k * q;
+ z = k * r;
+
+ } else {
+ fprintf(stderr, "Warning: pixel not visible\n");
+ }
+}
+
+double Conversions::compute_pixel_area_sqkm(double line, double column) {
+ double x1, x2;
+ double y1, y2;
+ double z1, z2;
+
+ compute_pixel_xyz(line-0.5, column-0.5, x1, y1, z1);
+ compute_pixel_xyz(line+0.5, column-0.5, x2, y2, z2);
+
+ double xlen = sqrt(SQR(x1 - x2) + SQR(y1 - y2) + SQR(z1 - z2));
+
+ compute_pixel_xyz(line-0.5, column+0.5, x2, y2, z2);
+
+ double ylen = sqrt(SQR(x1 - x2) + SQR(y1 - y2) + SQR(z1 - z2));
+
+ return xlen*ylen;
+}
+
+void Conversions::convert_geo_to_pixel(double longitude, double latitude,
unsigned int& line, unsigned int& column) {
+
+ latitude = latitude / 180.0 * M_PI;
+ longitude = longitude / 180.8 * M_PI;
+
+ double c_lat = atan(0.993243 * tan(latitude));
+ double r_l = rpol / sqrt(1 - 0.00675701*cos(c_lat)*cos(c_lat));
+ double r1 = altitude - r_l*cos(c_lat)*cos(longitude);
+ double r2 = -r_l*cos(c_lat)*sin(longitude);
+ double r3 = r_l*sin(c_lat);
+ double rn = sqrt(r1*r1 + r2*r2 + r3*r3);
+
+ double x = atan(-r2/r1) * (CFAC >> 16) + COFF;
+ double y = asin(-r3/rn) * (LFAC >> 16) + LOFF;
+
+ line = (unsigned int)floor(x + 0.5);
+ column = (unsigned int)floor(y + 0.5);
+}
+
+} // namespace msg_native_format
+
+
diff -Nur gdal/frmts/msgn/module-msg_reader/src/msg_reader_core.cc
gdal_patched/frmts/msgn/module-msg_reader/src/msg_reader_core.cc
--- gdal/frmts/msgn/module-msg_reader/src/msg_reader_core.cc 1970-01-01
02:00:00.000000000 +0200
+++ gdal_patched/frmts/msgn/module-msg_reader/src/msg_reader_core.cc
2005-07-01 16:07:23.567651952 +0200
@@ -0,0 +1,246 @@
+/* $Id: msg_reader_core.cc 10 2005-03-11 05:41:22Z fvdbergh $
+ Purpose: Base class for reading in the headers of MSG native images
+*/
+
+#include "include/msg_reader_core.h"
+#include "include/msg_basic_types.h"
+#include <stdio.h>
+#include <string.h>
+#include <math.h>
+
+#include <sys/time.h>
+
+#ifdef DEBUG
+#ifdef GDAL_SUPPORT
+#undef DEBUG
+#endif
+#endif
+
+#ifdef GDAL_SUPPORT
+#include "cpl_vsi.h"
+#else
+#define VSIFSeek(fp, pos, ref) fseek(fp, pos, ref)
+#define VSIFRead(p, bs, nb, fp) fread(p, bs, nb, ref)
+#endif
+
+namespace msg_native_format {
+
+const Blackbody_lut_type Msg_reader_core::Blackbody_LUT[MSG_NUM_CHANNELS+1] = {
+ {0,0,0}, // dummy channel
+ {0,0,0}, // N/A
+ {0,0,0}, // N/A
+ {0,0,0}, // N/A
+ {2569.094, 0.9959, 3.471},
+ {1598.566, 0.9963, 2.219},
+ {1362.142, 0.9991, 0.485},
+ {1149.083, 0.9996, 0.181},
+ {1034.345, 0.9999, 0.060},
+ { 930.659, 0.9983, 0.627},
+ { 839.661, 0.9988, 0.397},
+ { 752.381, 0.9981, 0.576},
+ {0,0,0} // N/A
+};
+
+
+Msg_reader_core::Msg_reader_core(const char* fname) {
+
+ FILE* fin = fopen(fname, "rb");
+ if (!fin) {
+ fprintf(stderr, "Could not open file %s\n", fname);
+ return;
+ }
+ read_metadata_block(fin);
+}
+
+Msg_reader_core::Msg_reader_core(FILE* fp) {
+ read_metadata_block(fp);
+}
+
+
+void Msg_reader_core::read_metadata_block(FILE* fin) {
+
+ VSIFRead(&_main_header, sizeof(_main_header), 1, fin);
+ VSIFRead(&_sec_header, sizeof(_sec_header), 1, fin);
+
+#ifdef DEBUG
+ // print out all the fields in the header
+ PH_DATA* hd = (PH_DATA*)&_main_header;
+ for (int i=0; i < 6; i++) {
+ to_string(*hd);
+ printf("[%02d] %s %s", i, hd->name, hd->value);
+ hd++;
+ }
+ PH_DATA_ID* hdi = (PH_DATA_ID*)&_main_header.dataSetIdentification;
+ for (int i=0; i < 5; i++) {
+ printf("%s %s %s", hdi->name, hdi->size, hdi->address);
+ hdi++;
+ }
+ hd = (PH_DATA*)(&_main_header.totalFileSize);
+ for (int i=0; i < 19; i++) {
+ to_string(*hd);
+ printf("[%02d] %s %s", i, hd->name, hd->value);
+ hd++;
+ }
+#endif // DEBUG
+
+ // extract data & header positions
+ for (int i=0; i < 5; i++) {
+ PH_DATA_ID* hdi = (PH_DATA_ID*)&_main_header.dataSetIdentification[i];
+ if (strncmp(hdi->name, "15Header", strlen("15Header")) == 0) {
+ sscanf(hdi->size, "%d", &_f_header_size);
+ sscanf(hdi->address, "%d", &_f_header_offset);
+ } else
+ if (strncmp(hdi->name, "15Data", strlen("15Data")) == 0) {
+ sscanf(hdi->size, "%d", &_f_data_size);
+ sscanf(hdi->address, "%d", &_f_data_offset);
+ }
+ }
+#ifdef DEBUG
+ printf("Data: %d %d\n", _f_data_offset, _f_data_size);
+ printf("Header: %d %d\n", _f_header_offset, _f_header_size);
+#endif // DEBUG
+
+ unsigned int lines;
+ sscanf(_sec_header.northLineSelectedRectangle.value, "%d", &_lines);
+ sscanf(_sec_header.southLineSelectedRectangle.value, "%d", &lines);
+ _line_start = lines;
+ _lines -= lines - 1;
+
+ unsigned int cols;
+ sscanf(_sec_header.westColumnSelectedRectangle.value, "%d", &_columns);
+ sscanf(_sec_header.eastColumnSelectedRectangle.value, "%d", &cols);
+ _col_start = cols;
+ _columns -= cols - 1;
+
+#ifdef DEBUG
+ printf("lines = %d, cols = %d\n", _lines, _columns);
+#endif // DEBUG
+
+ int records_per_line = 0;
+ for (unsigned int i=0; i < MSG_NUM_CHANNELS; i++) {
+ if (_sec_header.selectedBandIds.value[i] == 'X') {
+ _bands[i] = 1;
+ records_per_line += (i == (MSG_NUM_CHANNELS-1)) ? 3 : 1;
+ } else {
+ _bands[i] = 0;
+ }
+ }
+
+#ifdef DEBUG
+ printf("reading a total of %d records per line\n", records_per_line);
+#endif // DEBUG
+
+ // extract time fields, assume that SNIT is the correct field:
+ sscanf(_main_header.snit.value + 0, "%04d", &_year);
+ sscanf(_main_header.snit.value + 4, "%02d", &_month);
+ sscanf(_main_header.snit.value + 6, "%02d", &_day);
+ sscanf(_main_header.snit.value + 8, "%02d", &_hour);
+ sscanf(_main_header.snit.value + 10, "%02d", &_minute);
+
+ // read radiometric block
+ RADIOMETRIC_PROCCESSING_RECORD rad;
+ off_t offset = RADIOMETRICPROCESSING_RECORD_OFFSET + _f_header_offset +
sizeof(GP_PK_HEADER) + sizeof(GP_PK_SH1) + 1;
+ VSIFSeek(fin, offset, SEEK_SET);
+ VSIFRead(&rad, sizeof(RADIOMETRIC_PROCCESSING_RECORD), 1, fin);
+ to_native(rad);
+ memcpy((void*)_calibration, (void*)&rad.level1_5ImageCalibration,
sizeof(_calibration));
+
+#ifdef DEBUG
+ for (unsigned int i=0; i < MSG_NUM_CHANNELS; i++) {
+ if (_calibration[i].cal_slope < 0 || _calibration[i].cal_slope > 0.4) {
+ printf("Warning: calibration slope (%lf) out of nominal range. MSG
reader probably broken\n", _calibration[i].cal_slope);
+
+ }
+ if (_calibration[i].cal_offset > 0 || _calibration[i].cal_offset < -20) {
+ printf("Warning: calibration offset (%lf) out of nominal range. MSG
reader probably broken\n", _calibration[i].cal_offset);
+ }
+ }
+#endif
+
+ // read image description block
+ IMAGE_DESCRIPTION_RECORD idr;
+ offset = RADIOMETRICPROCESSING_RECORD_OFFSET -
IMAGEDESCRIPTION_RECORD_LENGTH + _f_header_offset + sizeof(GP_PK_HEADER) +
sizeof(GP_PK_SH1) + 1;
+ VSIFSeek(fin, offset, SEEK_SET);
+ VSIFRead(&idr, sizeof(IMAGE_DESCRIPTION_RECORD), 1, fin);
+ to_native(idr);
+ _line_dir_step = idr.referencegrid_visir.lineDirGridStep;
+ _col_dir_step = idr.referencegrid_visir.columnDirGridStep;
+
+
+ // Rather convoluted, but this code is required to compute the real data
block sizes
+ // It does this by reading in the first line of every band, to get to the
packet size field
+ GP_PK_HEADER gp_header;
+ GP_PK_SH1 sub_header;
+ SUB_VISIRLINE visir_line;
+
+ VSIFSeek(fin, _f_data_offset, SEEK_SET);
+
+ _hrv_packet_size = 0;
+ _interline_spacing = 0;
+ visir_line.channelId = 0;
+
+ int scanned_bands[MSG_NUM_CHANNELS];
+ int band_count = 0;
+ for (unsigned int i=0; i < MSG_NUM_CHANNELS; i++) {
+ scanned_bands[i] = _bands[i];
+ band_count += _bands[i];
+ }
+
+ unsigned int inter_line_spacing = 0;
+ do {
+ VSIFRead(&gp_header, sizeof(GP_PK_HEADER), 1, fin);
+ VSIFRead(&sub_header, sizeof(GP_PK_SH1), 1, fin);
+ VSIFRead(&visir_line, sizeof(SUB_VISIRLINE), 1, fin);
+ to_native(visir_line);
+ to_native(gp_header);
+ // skip over the actual line data
+ VSIFSeek(fin,
+ gp_header.packetLength - (sizeof(GP_PK_SH1) + sizeof(SUB_VISIRLINE)
- 1),
+ SEEK_CUR
+ );
+
+ if (scanned_bands[visir_line.channelId - 1]) {
+ scanned_bands[visir_line.channelId - 1] = 0;
+ band_count--;
+
+ if (visir_line.channelId != 12) { // not the HRV channel
+ _visir_bytes_per_line = gp_header.packetLength -
(sizeof(GP_PK_SH1) + sizeof(SUB_VISIRLINE) - 1);
+ _visir_packet_size = gp_header.packetLength +
sizeof(GP_PK_HEADER) + 1;
+ _interline_spacing += _visir_packet_size;
+ } else {
+ _hrv_packet_size = gp_header.packetLength +
sizeof(GP_PK_HEADER) + 1;
+ _interline_spacing += 3*_hrv_packet_size;
+ }
+ }
+ } while (band_count > 0);
+}
+
+#ifndef GDAL_SUPPORT
+
+int Msg_reader_core::_chan_to_idx(Msg_channel_names channel) {
+ unsigned int idx = 0;
+ while (idx < MSG_NUM_CHANNELS) {
+ if ( (1 << (idx+1)) == (int)channel ) {
+ return idx;
+ }
+ idx++;
+ }
+ return 0;
+}
+
+void Msg_reader_core::get_pixel_geo_coordinates(unsigned int line, unsigned int
column, double& longitude, double& latitude) {
+ Conversions::convert_pixel_to_geo((unsigned int)(line + _line_start),
(unsigned int)(column + _col_start), longitude, latitude);
+}
+
+void Msg_reader_core::get_pixel_geo_coordinates(double line, double column,
double& longitude, double& latitude) {
+ Conversions::convert_pixel_to_geo(line + _line_start, column + _col_start,
longitude, latitude);
+}
+
+double Msg_reader_core::compute_pixel_area_sqkm(double line, double column) {
+ return Conversions::compute_pixel_area_sqkm(line + _line_start, column +
_col_start);
+}
+
+#endif // GDAL_SUPPORT
+
+} // namespace msg_native_format
+
diff -Nur gdal/frmts/msgn/msgndataset.cpp gdal_patched/frmts/msgn/msgndataset.cpp
--- gdal/frmts/msgn/msgndataset.cpp 1970-01-01 02:00:00.000000000 +0200
+++ gdal_patched/frmts/msgn/msgndataset.cpp 2005-07-01 16:37:32.216695504 +0200
@@ -0,0 +1,401 @@
+/******************************************************************************
+ * $Id: $
+ *
+ * Project: MSG Native Reader
+ * Purpose: All code for EUMETSAT Archive format reader
+ * Author: Frans van den Bergh, fvdbergh@csir.co.za
+ *
+ ******************************************************************************
+ * Copyright (c) 2005, Frans van den Bergh <fvdbergh@csir.co.za>
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ ******************************************************************************
+ */
+
+#include "gdal_priv.h"
+#include "ogr_spatialref.h"
+
+#include "module-msg_reader/include/msg_reader_core.h"
+using namespace msg_native_format;
+
+CPL_CVSID("$Id: msgndataset.cpp $");
+
+CPL_C_START
+void GDALRegister_MSGN(void);
+CPL_C_END
+
+class MSGNRasterBand;
+
+/************************************************************************/
+/* ==================================================================== */
+/* MSGNDataset */
+/* ==================================================================== */
+/************************************************************************/
+
+class MSGNDataset : public GDALDataset
+{
+ friend class MSGNRasterBand;
+
+ FILE *fp;
+ GByte abyHeader[1012];
+
+ Msg_reader_core* msg_reader_core;
+ double adfGeoTransform[6];
+ char *pszProjection;
+
+ public:
+ MSGNDataset();
+ ~MSGNDataset();
+
+ static GDALDataset *Open( GDALOpenInfo * );
+
+ CPLErr GetGeoTransform( double * padfTransform );
+ const char *GetProjectionRef();
+
+ static void double2hex(double a, char* s);
+};
+
+/************************************************************************/
+/* ==================================================================== */
+/* MSGNRasterBand */
+/* ==================================================================== */
+/************************************************************************/
+
+class MSGNRasterBand : public GDALRasterBand
+{
+ friend class MSGNDataset;
+
+ unsigned int visir_packet_size;
+ unsigned int visir_bytes_per_line;
+ unsigned int interline_spacing;
+ double GetNoDataValue (int *pbSuccess=NULL) {
+ if (pbSuccess) {
+ *pbSuccess = 1;
+ }
+ return MSGN_NODATA_VALUE;
+ }
+
+ static const unsigned short MSGN_NODATA_VALUE;
+
+ public:
+
+ MSGNRasterBand( MSGNDataset *, int );
+
+ virtual CPLErr IReadBlock( int, int, void * );
+};
+
+const unsigned short MSGNRasterBand::MSGN_NODATA_VALUE = 0;
+
+/************************************************************************/
+/* MSGNRasterBand() */
+/************************************************************************/
+
+MSGNRasterBand::MSGNRasterBand( MSGNDataset *poDS, int nBand )
+
+{
+ this->poDS = poDS;
+ this->nBand = nBand;
+
+ eDataType = GDT_UInt16;
+
+ nBlockXSize = poDS->GetRasterXSize();
+ nBlockYSize = 1;
+
+ visir_packet_size = poDS->msg_reader_core->get_visir_packet_size();
+ visir_bytes_per_line = poDS->msg_reader_core->get_visir_bytes_per_line();
+ interline_spacing = poDS->msg_reader_core->get_interline_spacing();
+}
+
+/************************************************************************/
+/* IReadBlock() */
+/************************************************************************/
+
+CPLErr MSGNRasterBand::IReadBlock( int nBlockXOff, int nBlockYOff,
+ void * pImage )
+
+{
+ MSGNDataset *poGDS = (MSGNDataset *) poDS;
+ char *pszRecord;
+
+ unsigned int packet_length = visir_bytes_per_line + sizeof(SUB_VISIRLINE);
+ unsigned int packet_offset = poGDS->msg_reader_core->get_f_data_offset() +
+ interline_spacing*nBlockYOff + (nBand-1)*visir_packet_size +
+ (visir_packet_size - packet_length);
+
+ VSIFSeek( poGDS->fp, packet_offset, SEEK_SET );
+
+ pszRecord = (char *) CPLMalloc(packet_length);
+ size_t nread = VSIFRead( pszRecord, 1, packet_length, poGDS->fp );
+
+ SUB_VISIRLINE* p = (SUB_VISIRLINE*) pszRecord;
+ to_native(*p);
+
+ if (p->lineValidity != 1) {
+ for (int c=0; c < nBlockXSize; c++) {
+ ((short int *)pImage)[c] = MSGN_NODATA_VALUE;
+ }
+ }
+
+ if ( nread != packet_length ||
+ (p->lineNumberInVisirGrid - poGDS->msg_reader_core->get_line_start())
!= (unsigned int)nBlockYOff ) {
+ CPLFree( pszRecord );
+
+ CPLError( CE_Failure, CPLE_AppDefined, "MSGN Scanline corrupt." );
+
+ return CE_Failure;
+ }
+
+ // unpack the 10-bit values into 16-bit unsigned short ints
+ unsigned char *cptr = (unsigned char*)pszRecord +
+ (packet_length - visir_bytes_per_line);
+ int bitsLeft = 8;
+ unsigned short value = 0;
+ for (int c=0; c < nBlockXSize; c++) {
+ value = 0;
+ for (int bit=0; bit < 10; bit++) {
+ value <<= 1;
+ if (*cptr & 128) {
+ value |= 1;
+ }
+ *cptr <<= 1;
+ bitsLeft--;
+ if (bitsLeft == 0) {
+ cptr++;
+ bitsLeft = 8;
+ }
+ }
+ ((short int *)pImage)[c] = value;
+ }
+
+ return CE_None;
+}
+
+/************************************************************************/
+/* ==================================================================== */
+/* MSGNDataset */
+/* ==================================================================== */
+/************************************************************************/
+
+MSGNDataset::MSGNDataset() {
+ pszProjection = CPLStrdup("");
+}
+
+/************************************************************************/
+/* ~MSGNDataset() */
+/************************************************************************/
+
+MSGNDataset::~MSGNDataset()
+
+{
+ if( fp != NULL )
+ VSIFClose( fp );
+
+ CPLFree(pszProjection);
+}
+
+/************************************************************************/
+/* GetGeoTransform() */
+/************************************************************************/
+
+void MSGNDataset::double2hex(double a, char* out_str) {
+ const char hexlut[] = "0123456789abcdef";
+ char* p = (char*)&a;
+ char* s = out_str;
+ for(unsigned int i=0; i < sizeof(double); i++) {
+ *s = hexlut[(*p >> 4) & 0x0f];
+ s++;
+ *s = hexlut[*p & 0x0f];
+ s++;
+ p++;
+ }
+ *s = 0;
+}
+
+/************************************************************************/
+/* GetGeoTransform() */
+/************************************************************************/
+
+CPLErr MSGNDataset::GetGeoTransform( double * padfTransform )
+
+{
+
+ for (int i=0; i < 6; i++) {
+ padfTransform[i] = adfGeoTransform[i];
+ }
+
+ return CE_None;
+}
+
+/************************************************************************/
+/* GetProjectionRef() */
+/************************************************************************/
+
+const char *MSGNDataset::GetProjectionRef()
+
+{
+ return ( pszProjection );
+}
+
+/************************************************************************/
+/* Open() */
+/************************************************************************/
+
+GDALDataset *MSGNDataset::Open( GDALOpenInfo * poOpenInfo )
+
+{
+/* -------------------------------------------------------------------- */
+/* Before trying MSGNOpen() we first verify that there is at */
+/* least one "\n#keyword" type signature in the first chunk of */
+/* the file. */
+/* -------------------------------------------------------------------- */
+ if( poOpenInfo->fp == NULL || poOpenInfo->nHeaderBytes < 50 )
+ return NULL;
+
+ /* check if this is a "NATIVE" MSG format image */
+ if( !EQUALN((char *)poOpenInfo->pabyHeader,"FormatName :
NATIVE",36) )
+ {
+ return NULL;
+ }
+
+/* -------------------------------------------------------------------- */
+/* Create a corresponding GDALDataset. */
+/* -------------------------------------------------------------------- */
+ MSGNDataset *poDS;
+
+ poDS = new MSGNDataset();
+
+ poDS->fp = poOpenInfo->fp;
+ poOpenInfo->fp = NULL;
+
+/* -------------------------------------------------------------------- */
+/* Read the header. */
+/* -------------------------------------------------------------------- */
+ // first reset the file pointer, then hand over to the msg_reader_core
+ VSIFSeek( poDS->fp, 0, SEEK_SET );
+
+ poDS->msg_reader_core = new Msg_reader_core(poDS->fp);
+
+ poDS->nRasterXSize = poDS->msg_reader_core->get_columns();
+ poDS->nRasterYSize = poDS->msg_reader_core->get_lines();
+
+/* -------------------------------------------------------------------- */
+/* Create band information objects. */
+/* -------------------------------------------------------------------- */
+ unsigned int band_count = 1;
+ unsigned char* bands = poDS->msg_reader_core->get_band_map();
+ for (unsigned int i=0; i < MSG_NUM_CHANNELS-1; i++) {
+ if (bands[i]) {
+ poDS->SetBand( band_count, new MSGNRasterBand( poDS, band_count ));
+ band_count++;
+ }
+ }
+
+ double pixel_gsd_x = 1000 * poDS->msg_reader_core->get_col_dir_step(); //
convert from km to m
+ double pixel_gsd_y = 1000 * poDS->msg_reader_core->get_line_dir_step(); //
convert from km to m
+ double origin_x = -pixel_gsd_x * (-(Conversions::nlines / 2.0) +
poDS->msg_reader_core->get_col_start());
+ double origin_y = -pixel_gsd_y * ((Conversions::nlines / 2.0) -
poDS->msg_reader_core->get_line_start());
+
+ poDS->adfGeoTransform[0] = origin_x;
+ poDS->adfGeoTransform[1] = -pixel_gsd_x;
+ poDS->adfGeoTransform[2] = 0.0;
+
+ poDS->adfGeoTransform[3] = origin_y;
+ poDS->adfGeoTransform[4] = 0.0;
+ poDS->adfGeoTransform[5] = pixel_gsd_y;
+
+ OGRSpatialReference oSRS;
+
+ oSRS.SetProjCS("Geostationary projection (MSG)");
+ oSRS.SetGEOS( 0, 35785831, 0, 0 );
+ oSRS.SetGeogCS(
+ "MSG Ellipsoid",
+ "MSG_DATUM",
+ "MSG_SPHEROID",
+ Conversions::rpol * 1000.0,
+ 1 / ( 1 - Conversions::rpol/Conversions::req)
+ );
+
+ oSRS.exportToWkt( &(poDS->pszProjection) );
+
+ CALIBRATION* cal = poDS->msg_reader_core->get_calibration_parameters();
+ char tagname[30];
+ char field[300];
+ char hexvalue1[30];
+ char hexvalue2[30];
+
+ poDS->SetMetadataItem("Radiometric parameters format", "offset slope
hex_offset hex_slope");
+ for (unsigned int i=0; i < MSG_NUM_CHANNELS; i++) {
+ sprintf(tagname, "ch%02d_cal", i+1);
+
+ MSGNDataset::double2hex(cal[i].cal_offset, hexvalue1);
+ MSGNDataset::double2hex(cal[i].cal_slope, hexvalue2);
+
+ sprintf(field, "%.12e %.12e %s %s", cal[i].cal_offset,
cal[i].cal_slope, hexvalue1, hexvalue2);
+ poDS->SetMetadataItem(tagname, field);
+ }
+
+ poDS->SetMetadataItem("Blackbody parameters format", "vc A B");
+ for (unsigned int i=4; i < MSG_NUM_CHANNELS-1; i++) {
+ sprintf(tagname, "ch%02d_blackbody", i);
+ sprintf(field, "%.4f %.4f %.4f",
+ Msg_reader_core::Blackbody_LUT[i].vc,
+ Msg_reader_core::Blackbody_LUT[i].A,
+ Msg_reader_core::Blackbody_LUT[i].B
+ );
+ poDS->SetMetadataItem(tagname, field);
+ }
+
+ sprintf(field, "%04d%02d%02d/%02d:%02d",
+ poDS->msg_reader_core->get_year(),
+ poDS->msg_reader_core->get_month(),
+ poDS->msg_reader_core->get_day(),
+ poDS->msg_reader_core->get_hour(),
+ poDS->msg_reader_core->get_minute()
+ );
+ poDS->SetMetadataItem("Date/Time", field);
+
+
+ return( poDS );
+}
+
+/************************************************************************/
+/* GDALRegister_MSGN() */
+/************************************************************************/
+
+void GDALRegister_MSGN()
+
+{
+ GDALDriver *poDriver;
+
+ if( GDALGetDriverByName( "MSGN" ) == NULL )
+ {
+ poDriver = new GDALDriver();
+
+ poDriver->SetDescription( "MSGN" );
+ poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
+ "EUMETSAT Archive native (.nat)" );
+ poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
+ "frmt_various.html#MSGN" );
+ poDriver->SetMetadataItem( GDAL_DMD_EXTENSION, "nat" );
+
+ poDriver->pfnOpen = MSGNDataset::Open;
+
+ GetGDALDriverManager()->RegisterDriver( poDriver );
+ }
+}
diff -Nur gdal/gcore/gdal_frmts.h gdal_patched/gcore/gdal_frmts.h
--- gdal/gcore/gdal_frmts.h 2005-05-19 22:43:08.000000000 +0200
+++ gdal_patched/gcore/gdal_frmts.h 2005-07-01 16:13:31.659693456 +0200
@@ -206,6 +206,7 @@
void CPL_DLL GDALRegister_IDA(void);
void CPL_DLL GDALRegister_NDF(void);
void CPL_DLL GDALRegister_RMF(void);
+void CPL_DLL GDALRegister_MSGN(void);
CPL_C_END
#endif /* ndef GDAL_FRMTS_H_INCLUDED */
Attachments (1)
Change History (4)
comment:2 by , 18 years ago
Frank, the docs of this driver are damaged. Find attached a patch resolving the issue. Best regards Markus
Note:
See TracTickets
for help on using tickets.
