geo/tcLandsatData.cpp

   1 /***************************************************************************
   2  *   This file is part of Tecorrec.                                        *
   3  *   Copyright 2008 James Hogan <james@albanarts.com>                      *
   4  *                                                                         *
   5  *   Tecorrec is free software: you can redistribute it and/or modify      *
   6  *   it under the terms of the GNU General Public License as published by  *
   7  *   the Free Software Foundation, either version 2 of the License, or     *
   8  *   (at your option) any later version.                                   *
   9  *                                                                         *
  10  *   Tecorrec is distributed in the hope that it will be useful,           *
  11  *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
  12  *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
  13  *   GNU General Public License for more details.                          *
  14  *                                                                         *
  15  *   You should have received a copy of the GNU General Public License     *
  16  *   along with Tecorrec.  If not, write to the Free Software Foundation,  *
  17  *   Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.   *
  18  ***************************************************************************/
  19
  20 /**
  21  * @file tcLandsatData.cpp
  22  * @brief Landsat image data.
  23  */
  24
  25 #include "tcLandsatData.h"
  26 #include "tcLandsatMetaData.h"
  27 #include "tcSrtmModel.h"
  28 #include "tcChannelManager.h"
  29 #include "tcChannelFile.h"
  30 #include "tcChannelChromaticity.h"
  31 #include "tcShadowFreeChromaticities.h"
  32 #include "tcIlluminantDiscontinuity.h"
  33 #include "tcShadowClassification.h"
  34 #include "tcNegativeProduct.h"
  35 #include "tcLambertianShading.h"
  36 #include "tcRaytracedShadowMap.h"
  37 #include "tcShadowDepth.h"
  38 #include "tcPixelOp.h"
  39
  40 #include <QDir>
  41 #include <QFile>
  42 #include <QColor>
  43
  44 /*
  45  * Constructors + destructor
  46  */
  47
  48 /// Load landsat data given the path.
  49 tcLandsatData::tcLandsatData(const QString& path, tcSrtmModel* dem)
  50 : tcGeoImageData()
  51 , m_meta(0)
  52 , m_shadowClassification(0)
  53 , m_shading(0)
  54 {
  55   // Find name of txt file automatically
  56   QDir dir(path);
  57   QStringList textFiles = dir.entryList(QStringList() << "*_MTL.txt" << "*.met");
  58   foreach (QString textFile, textFiles)
  59   {
  60     QFile metaFile(path + textFile);
  61     if (metaFile.open(QIODevice::ReadOnly))
  62     {
  63       m_meta = new tcLandsatMetaData(metaFile);
  64       tcLandsatMetaData::Reference meta = m_meta;
  65
  66       // Landsat 7 data looks like this
  67       tcLandsatMetaData::Reference l1_metadata_file = meta("L1_METADATA_FILE");
  68       if (l1_metadata_file)
  69       {
  70         tcLandsatMetaData::Reference productMetadata = l1_metadata_file("PRODUCT_METADATA");
  71
  72         setName(QObject::tr("%1 %2 (%3)").arg(productMetadata["SPACECRAFT_ID"].toString())
  73                                          .arg(productMetadata["SENSOR_ID"].toString())
  74                                          .arg(productMetadata["ACQUISITION_DATE"].toString()));
  75
  76         // Load coordinates
  77         m_coordinates[0][0] = tcGeo(productMetadata["PRODUCT_UL_CORNER_LON"].toDouble() * M_PI/180.0,
  78                                     productMetadata["PRODUCT_UL_CORNER_LAT"].toDouble() * M_PI/180.0);
  79         m_coordinates[1][0] = tcGeo(productMetadata["PRODUCT_UR_CORNER_LON"].toDouble() * M_PI/180.0,
  80                                     productMetadata["PRODUCT_UR_CORNER_LAT"].toDouble() * M_PI/180.0);
  81         m_coordinates[0][1] = tcGeo(productMetadata["PRODUCT_LL_CORNER_LON"].toDouble() * M_PI/180.0,
  82                                     productMetadata["PRODUCT_LL_CORNER_LAT"].toDouble() * M_PI/180.0);
  83         m_coordinates[1][1] = tcGeo(productMetadata["PRODUCT_LR_CORNER_LON"].toDouble() * M_PI/180.0,
  84                                     productMetadata["PRODUCT_LR_CORNER_LAT"].toDouble() * M_PI/180.0);
  85
  86         tcLandsatMetaData::Reference productParameters = l1_metadata_file("PRODUCT_PARAMETERS");
  87         setSunDirection(tcGeo(M_PI - M_PI/180.0 * productParameters["SUN_AZIMUTH"].toDouble(),
  88                               M_PI/180.0 * productParameters["SUN_ELEVATION"].toDouble()));
  89
  90         // Names and descriptions of channels
  91         QStringList channelNames;
  92         channelNames << QObject::tr("Blue Green");
  93         channelNames << QObject::tr("Green");
  94         channelNames << QObject::tr("Red");
  95         channelNames << QObject::tr("Near IR");
  96         channelNames << QObject::tr("Mid IR 1");
  97         channelNames << QObject::tr("Thermal IR Low");
  98         channelNames << QObject::tr("Thermal IR High");
  99         channelNames << QObject::tr("Mid IR 2");
 100         channelNames << QObject::tr("Panchromatic");
 101
 102         // Get the filenames of the bands
 103         QString bandCombo = productMetadata["BAND_COMBINATION"].toString();
 104         QChar last = 0;
 105         int sameCount = 0;
 106         for (int band = 0; band < bandCombo.size(); ++band)
 107         {
 108           QChar code = bandCombo.at(band);
 109           QChar next = (band+1 < bandCombo.size() ? bandCombo.at(band+1) : 0);
 110           QString bandName;
 111           QString description;
 112           if (code == next || last == code)
 113           {
 114             ++sameCount;
 115             bandName = QString("BAND%1%2_FILE_NAME").arg(code).arg(sameCount);
 116             description = QObject::tr("Band %1.%2").arg(code).arg(sameCount);
 117           }
 118           else
 119           {
 120             sameCount = 0;
 121             bandName = QString("BAND%1_FILE_NAME").arg(code);
 122             description = QObject::tr("Band %1").arg(code);
 123           }
 124           tcChannelFile* channel = new tcChannelFile(path + productMetadata[bandName].toString(),
 125                                                          channelNames[band], description);
 126           channelManager()->addChannel(channel);
 127           // Obtain the transformation from the first of the files.
 128           if (0 == band)
 129           {
 130             setTexToGeo(channel->texToGeo());
 131           }
 132           last = code;
 133         }
 134       }
 135       // Landsat 5 data is a little different
 136       else
 137       {
 138         tcLandsatMetaData::Reference inventoryMetaData = meta("INVENTORYMETADATA");
 139         tcLandsatMetaData::Reference sensorCharacteristics = inventoryMetaData("SENSORCHARACTERISTIC")("SENSORCHARACTERISTICCONTAINER",0);
 140         tcLandsatMetaData::Reference singleDateTime = inventoryMetaData("SINGLEDATETIME");
 141
 142         setName(QObject::tr("%1 %2 (%3)").arg(sensorCharacteristics("PLATFORMSHORTNAME",0)["VALUE"].toString())
 143                                          .arg(sensorCharacteristics("SENSORSHORTNAME",0)["VALUE"].toString())
 144                                          .arg(singleDateTime("CALENDARDATE",0)["VALUE"].toString()));
 145
 146         // Load additional attributes into a nice hash
 147         tcLandsatMetaData::Reference additionalAttributes = inventoryMetaData("ADDITIONALATTRIBUTES");
 148         QHash<QString, QVariant> attributes;
 149         int numAttributes = additionalAttributes.hasObject("ADDITIONALATTRIBUTESCONTAINER");
 150         for (int i = 0; i < numAttributes; ++i)
 151         {
 152           tcLandsatMetaData::Reference container = additionalAttributes("ADDITIONALATTRIBUTESCONTAINER", i);
 153           QString name = container("ADDITIONALATTRIBUTENAME",0)["VALUE"].toString();
 154           attributes[name] = container("INFORMATIONCONTENT")("PARAMETERVALUE", 0)["VALUE"];
 155         }
 156
 157         // Sun direction
 158         setSunDirection(tcGeo(M_PI - M_PI/180.0 * attributes["SolarAzimuth"].toDouble(),
 159                                      M_PI/180.0 * attributes["SolarElevation"].toDouble()));
 160
 161         // Lon/Lat coordinates
 162         tcLandsatMetaData::Reference gRingPoint = inventoryMetaData("SPATIALDOMAINCONTAINER")
 163                                                                    ("HORIZONTALSPATIALDOMAINCONTAINER")
 164                                                                    ("GPOLYGON")
 165                                                                    ("GPOLYGONCONTAINER",0)
 166                                                                    ("GRINGPOINT");
 167         const QList<QVariant>& lons = gRingPoint("GRINGPOINTLONGITUDE",0)["VALUE"].toList();
 168         const QList<QVariant>& lats = gRingPoint("GRINGPOINTLATITUDE",0)["VALUE"].toList();
 169
 170         m_coordinates[0][0] = tcGeo(lons[0].toDouble() * M_PI/180.0,
 171                                     lats[0].toDouble() * M_PI/180.0);
 172         m_coordinates[1][0] = tcGeo(lons[1].toDouble() * M_PI/180.0,
 173                                     lats[1].toDouble() * M_PI/180.0);
 174         m_coordinates[0][1] = tcGeo(lons[3].toDouble() * M_PI/180.0,
 175                                     lats[3].toDouble() * M_PI/180.0);
 176         m_coordinates[1][1] = tcGeo(lons[2].toDouble() * M_PI/180.0,
 177                                     lats[2].toDouble() * M_PI/180.0);
 178
 179         // Names and descriptions of channels
 180         QStringList channelNames;
 181         channelNames << QObject::tr("Blue Green");
 182         channelNames << QObject::tr("Green");
 183         channelNames << QObject::tr("Red");
 184         channelNames << QObject::tr("Near IR");
 185         channelNames << QObject::tr("Mid IR 1");
 186         channelNames << QObject::tr("Thermal IR");
 187         channelNames << QObject::tr("Mid IR 2");
 188
 189         // Get input files and filter tif files
 190         int band = 0;
 191         foreach (QVariant file, inventoryMetaData("INPUTGRANULE")("INPUTPOINTER",0)["VALUE"].toList())
 192         {
 193           static QRegExp tifExt("^.*\\.tif$");
 194           if (tifExt.exactMatch(file.toString()))
 195           {
 196             tcChannelFile* channel = new tcChannelFile(path + file.toString(),
 197                                                        channelNames[band], QString("Band %1").arg(band));
 198             channelManager()->addChannel(channel);
 199             // Obtain the transformation from the first of the files.
 200             if (0 == band)
 201             {
 202               setTexToGeo(channel->texToGeo());
 203             }
 204             ++band;
 205           }
 206         }
 207       }
 208
 209
 210
 211       QList<tcChannel*> mainColourChannels;
 212       for (int i = 0; i < channelManager()->numChannels(); ++i)
 213       {
 214         mainColourChannels += channelManager()->channel(i);
 215       }
 216       m_shadowClassification = new tcNegativeProduct(mainColourChannels);
 217       channelManager()->addChannel(m_shadowClassification);
 218
 219       tcShadowDepth* shadowDepth1 = new tcShadowDepth(m_shadowClassification, dem, this, false);
 220       channelManager()->addChannel(shadowDepth1);
 221       tcShadowDepth* shadowDepth2 = new tcShadowDepth(m_shadowClassification, dem, this, true);
 222       channelManager()->addChannel(shadowDepth2);
 223
 224       tcLambertianShading* diffuse = new tcLambertianShading(channelManager()->channel(0), dem, this);
 225       channelManager()->addChannel(diffuse);
 226
 227       tcRaytracedShadowMap* raytrace = new tcRaytracedShadowMap(channelManager()->channel(0), dem, this);
 228       channelManager()->addChannel(raytrace);
 229
 230       m_shading = new tcPixelOp(tcPixelOp::Add, QList<tcChannel*>()
 231         //<< channelManager()->channel(3)
 232         //<< channelManager()->channel(4)
 233         //<< channelManager()->channel(5)
 234         << channelManager()->channel(6)
 235         //<< channelManager()->channel(7)
 236         //<< channelManager()->channel(8)
 237         //, 1.0f / 800
 238         , 1.0f / 125
 239         );
 240       channelManager()->addChannel(m_shading);
 241
 242       /*
 243       QList<tcChannel*> chromaticities;
 244       for (int i = 0; i <= 8; ++i)
 245       {
 246         chromaticities += new tcChannelChromaticity(channelManager()->channel(0), channelManager()->channel(i));
 247         chromaticities += new tcChannelChromaticity(channelManager()->channel(1), channelManager()->channel(i));
 248         chromaticities += new tcChannelChromaticity(channelManager()->channel(2), channelManager()->channel(i));
 249         chromaticities += new tcChannelChromaticity(channelManager()->channel(3), channelManager()->channel(i));
 250         chromaticities += new tcChannelChromaticity(channelManager()->channel(4), channelManager()->channel(i));
 251         chromaticities += new tcChannelChromaticity(channelManager()->channel(5), channelManager()->channel(i));
 252         chromaticities += new tcChannelChromaticity(channelManager()->channel(6), channelManager()->channel(i));
 253         chromaticities += new tcChannelChromaticity(channelManager()->channel(7), channelManager()->channel(i));
 254         chromaticities += new tcChannelChromaticity(channelManager()->channel(8), channelManager()->channel(i));
 255       }
 256       channelManager()->addChannels(chromaticities);
 257       */
 258
 259       /*
 260       tcShadowFreeChromaticities* sfc = new tcShadowFreeChromaticities(chromaticities);
 261       channelManager()->addChannels(sfc->channels());
 262       tcIlluminantDiscontinuity* idm = new tcIlluminantDiscontinuity(chromaticities);
 263       channelManager()->addChannels(idm->channels());
 264       tcShadowClassification* sc = new tcShadowClassification(idm, channelManager()->channel(4));
 265       channelManager()->addChannel(sc);
 266       */
 267
 268       // Only load one metafile
 269       break;
 270     }
 271   }
 272 }
 273
 274 /// Destructor.
 275 tcLandsatData::~tcLandsatData()
 276 {
 277 }
 278
 279 /*
 280  * Accessors
 281  */
 282
 283 /// Get the shadow classification channel.
 284 tcChannel* tcLandsatData::shadowClassification()
 285 {
 286   return m_shadowClassification;
 287 }
 288
 289 /// Get shading channel.
 290 tcChannel* tcLandsatData::shading()
 291 {
 292   return m_shading;
 293 }