core/vil/file_formats/vil_geotiff_header.cxx
Go to the documentation of this file.
00001 #include "vil_geotiff_header.h"
00002 //:
00003 // \file
00004 #include <vcl_cassert.h>
00005 #include <vcl_cstdlib.h>
00006 #include <vcl_iostream.h>
00007 #include <geo_tiffp.h>
00008 #include <geotiffio.h>
00009 #include <geovalues.h>
00010 
00011 vil_geotiff_header::vil_geotiff_header(TIFF* tif) : tif_(tif)
00012 {
00013   if (tif) {
00014     gtif_ = GTIFNew(tif);
00015     if (gtif_) {
00016       GTIFPrint(gtif_, 0, 0);
00017 
00018       // read the header of the GeoDirectoryKey Tag
00019       int version[3];
00020       GTIFDirectoryInfo(gtif_, version, &number_of_geokeys_);
00021       key_directory_version_ = (unsigned short)version[0];
00022       key_revision_ =          (unsigned short)version[1];
00023       minor_revision_ =        (unsigned short)version[2];
00024     }
00025   }
00026 }
00027 
00028 bool vil_geotiff_header::gtif_tiepoints(vcl_vector<vcl_vector<double> > &tiepoints)
00029 {
00030   double* points=0;
00031   short count;
00032   if (TIFFGetField(tif_, GTIFF_TIEPOINTS, &count, &points) < 0)
00033     return false;
00034 
00035   // tiepoints are stored as 3d points (I,J,K)->(X,Y,Z)
00036   // where the point at location (I,J) at raster space with pixel value K
00037   // and (X,Y,Z) is a vector in model space
00038 
00039   // the number of values should be K*6
00040   assert((count % 6) == 0);
00041   for (unsigned short i=0; i<count; ) {
00042     vcl_vector<double> tiepoint(6);
00043     tiepoint[0] = points[i++];
00044     tiepoint[1] = points[i++];
00045     tiepoint[2] = points[i++];
00046     tiepoint[3] = points[i++];
00047     tiepoint[4] = points[i++];
00048     tiepoint[5] = points[i++];
00049     tiepoints.push_back(tiepoint);
00050   }
00051   return true;
00052 }
00053 
00054 bool vil_geotiff_header::gtif_pixelscale(double &scale_x, double &scale_y, double &scale_z)
00055 {
00056   double *data;
00057   short count;
00058   if (TIFFGetField(tif_, GTIFF_PIXELSCALE, &count, &data )) {
00059     assert (count == 3);
00060     scale_x = data[0];
00061     scale_y = data[1];
00062     scale_z = data[2];
00063     return true;
00064   }
00065   else
00066     return false;
00067 }
00068 
00069 bool vil_geotiff_header::gtif_trans_matrix (double* &trans_matrix)
00070 {
00071   short count;
00072   if (TIFFGetField(tif_, GTIFF_TRANSMATRIX, &count, &trans_matrix )) {
00073     assert (count == 16);
00074     return true;
00075   }
00076   else
00077     return false;
00078 }
00079 
00080 bool vil_geotiff_header::gtif_modeltype (modeltype_t& type)
00081 {
00082   geocode_t model;
00083   if (!GTIFKeyGet(gtif_, GTModelTypeGeoKey, &model, 0, 1))  {
00084     vcl_cerr << "NO Model Type defined!!!!\n";
00085     return false;
00086   }
00087   else {
00088     type = static_cast<modeltype_t> (model);
00089     return true;
00090   }
00091 }
00092 
00093 bool vil_geotiff_header::gtif_rastertype (rastertype_t& type)
00094 {
00095   geocode_t raster;
00096   if (!GTIFKeyGet(gtif_, GTRasterTypeGeoKey, &raster, 0, 1)) {
00097     vcl_cerr << "NO Raster Type, failure!!!!\n";
00098     return false;
00099   }
00100   else {
00101     type = static_cast<rastertype_t> (raster);
00102     return true;
00103   }
00104 }
00105 
00106 bool vil_geotiff_header::geounits (geounits_t& units)
00107 {
00108   short nGeogUOMLinear;
00109   if (!GTIFKeyGet(gtif_, GeogLinearUnitsGeoKey, &nGeogUOMLinear, 0, 1 )) {
00110     vcl_cerr << "NO GEOUNITS, failure!!!!\n";
00111     return false;
00112   }
00113   else {
00114     units = static_cast<geounits_t> (nGeogUOMLinear);
00115     return true;
00116   }
00117 }
00118 
00119 bool vil_geotiff_header::PCS_WGS84_UTM_zone(int &zone, GTIF_HEMISPH &hemisph) // hemisph is O for N, 1 for S
00120 {
00121   modeltype_t type;
00122   if (gtif_modeltype(type) && type == ModelTypeProjected)
00123   {
00124     void *value;
00125     int size;
00126     int length;
00127     tagtype_t ttype;
00128     get_key_value(ProjectedCSTypeGeoKey, &value, size, length, ttype);
00129 
00130     // WGS84 / UTM northern hemisphere:  326zz where zz is UTM zone number
00131     // WGS84 / UTM southern hemisphere:  327zz where zz is UTM zone number
00132     // we are looking for a short value, only one
00133     assert ((length == 1) && (ttype == TYPE_SHORT));
00134 
00135     short *val = static_cast<short*> (value);
00136     if ((*val < PCS_WGS84_UTM_zone_1N ) || ((*val > PCS_WGS84_UTM_zone_60S ))) {
00137       vcl_cerr << "NOT A WGS84_UTM!!!!\n";
00138       return false;
00139     }
00140 
00141 #if 0
00142      int zone;
00143      int hemisph; // O for N, 1 for S
00144 #endif // 0
00145     if ((*val >= PCS_WGS84_UTM_zone_1N) && (*val <= PCS_WGS84_UTM_zone_60N)) {
00146       zone = *val - 32600;
00147       hemisph = NORTH;
00148     }
00149     else if ((*val >= PCS_WGS84_UTM_zone_1S) && (*val <= PCS_WGS84_UTM_zone_60S)) {
00150       zone = *val - 32700;
00151       hemisph = SOUTH;
00152     }
00153     return true;
00154   }
00155   else {
00156     hemisph = UNDEF;
00157     return false;
00158   }
00159 }
00160 
00161 bool vil_geotiff_header::get_key_value(geokey_t key, void** value,
00162                                        int& size, int& length, tagtype_t& type)
00163 {
00164   length = GTIFKeyInfo(gtif_, key, &size, &type);
00165 
00166   // check if it is a valid key id
00167   if (length == 0) {
00168     // key is not defined
00169     vcl_cerr << "KeyID=" << (short int)key << " is not valid\n";
00170     return false;
00171   }
00172   else {
00173     *value = vcl_malloc(size*length);
00174     GTIFKeyGet(gtif_, key, *value, 0, length);
00175     return true;
00176   }
00177 }