Ticket #879 (closed enhancement: fixed)

Opened 9 years ago

Last modified 9 years ago

Meteosat Native Archive format (.nat) patch

Reported by: fvdbergh@… 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

frmt_msgn.patch Download (1.2 KB) - added by neteler@… 9 years ago.
patch to fix broken documentation

Change History

Changed 9 years ago by warmerdam

Frans,

I have committed most of this after rearranging things a bit.  Everything
is now in the gdal/frmts/msgn directory instead of in subdirectories.  I also
fixed some portability issues (use of htonl and htons) and other odds and ends.

However, now I need some test data!  Can you point me to some. 

PS. patches should be handled as attachments.  garbling of the patch due
to line wrapping made it a big hassle to apply.  I'm still not certain I haven't
botched something. 

Changed 9 years ago by neteler@…

Frank,

the docs of this driver are damaged.
Find attached a patch resolving the issue.

Best regards

Markus

Changed 9 years ago by neteler@…

patch to fix broken documentation

Changed 9 years ago by warmerdam

Markus, 

Docs patch applied manually. 

Thanks!
Note: See TracTickets for help on using tickets.