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 #include <QtDebug>
  49 /// Load landsat data given the path.
  50 tcLandsatData::tcLandsatData(const QString& path, tcSrtmModel* dem)
  51 : tcShadowClassifyingData()
  52 , m_meta(0)
  53 {
  54   // Find name of txt file automatically
  55   QDir dir(path);
  56   QStringList textFiles = dir.entryList(QStringList() << "*_MTL.txt" << "*.met");
  57   foreach (QString textFile, textFiles)
  58   {
  59     QFile metaFile(path + textFile);
  60     if (metaFile.open(QIODevice::ReadOnly))
  61     {
  62       m_meta = new tcLandsatMetaData(metaFile);
  63       tcLandsatMetaData::Reference meta = m_meta;
  64
  65       // Landsat 7 data looks like this
  66       tcLandsatMetaData::Reference l1_metadata_file = meta("L1_METADATA_FILE");
  67       if (l1_metadata_file)
  68       {
  69         tcLandsatMetaData::Reference productMetadata = l1_metadata_file("PRODUCT_METADATA");
  70         tcLandsatMetaData::Reference minMaxRadiance = l1_metadata_file("MIN_MAX_RADIANCE");
  71         tcLandsatMetaData::Reference minMaxPixelValue = l1_metadata_file("MIN_MAX_PIXEL_VALUE");
  72
  73         setName(QObject::tr("%1 %2 (%3)").arg(productMetadata["SPACECRAFT_ID"].toString())
  74                                          .arg(productMetadata["SENSOR_ID"].toString())
  75                                          .arg(productMetadata["ACQUISITION_DATE"].toString()));
  76
  77         // Load coordinates
  78         m_coordinates[0][0] = tcGeo(productMetadata["PRODUCT_UL_CORNER_LON"].toDouble() * M_PI/180.0,
  79                                     productMetadata["PRODUCT_UL_CORNER_LAT"].toDouble() * M_PI/180.0);
  80         m_coordinates[1][0] = tcGeo(productMetadata["PRODUCT_UR_CORNER_LON"].toDouble() * M_PI/180.0,
  81                                     productMetadata["PRODUCT_UR_CORNER_LAT"].toDouble() * M_PI/180.0);
  82         m_coordinates[0][1] = tcGeo(productMetadata["PRODUCT_LL_CORNER_LON"].toDouble() * M_PI/180.0,
  83                                     productMetadata["PRODUCT_LL_CORNER_LAT"].toDouble() * M_PI/180.0);
  84         m_coordinates[1][1] = tcGeo(productMetadata["PRODUCT_LR_CORNER_LON"].toDouble() * M_PI/180.0,
  85                                     productMetadata["PRODUCT_LR_CORNER_LAT"].toDouble() * M_PI/180.0);
  86
  87         tcLandsatMetaData::Reference productParameters = l1_metadata_file("PRODUCT_PARAMETERS");
  88         setSunDirection(tcGeo(M_PI - M_PI/180.0 * productParameters["SUN_AZIMUTH"].toDouble(),
  89                               M_PI/180.0 * productParameters["SUN_ELEVATION"].toDouble()),
  90                         m_coordinates[1][0]);
  91
  92         // Names and descriptions of channels
  93         QStringList channelNames;
  94         channelNames << QObject::tr("Blue Green");
  95         channelNames << QObject::tr("Green");
  96         channelNames << QObject::tr("Red");
  97         channelNames << QObject::tr("Near IR");
  98         channelNames << QObject::tr("Mid IR 1");
  99         channelNames << QObject::tr("Thermal IR Low");
 100         channelNames << QObject::tr("Thermal IR High");
 101         channelNames << QObject::tr("Mid IR 2");
 102         channelNames << QObject::tr("Panchromatic");
 103
 104         // Get the filenames of the bands
 105         QString bandCombo = productMetadata["BAND_COMBINATION"].toString();
 106         QChar last = 0;
 107         int sameCount = 0;
 108         for (int band = 0; band < bandCombo.size(); ++band)
 109         {
 110           QChar code = bandCombo.at(band);
 111           QChar next = (band+1 < bandCombo.size() ? bandCombo.at(band+1) : 0);
 112           QString bandNumber;
 113           QString bandName;
 114           QString description;
 115           if (code == next || last == code)
 116           {
 117             ++sameCount;
 118             bandNumber = QString("%1%2").arg(code).arg(sameCount);
 119             description = QObject::tr("Band %1.%2").arg(code).arg(sameCount);
 120           }
 121           else
 122           {
 123             sameCount = 0;
 124             bandNumber = QString("%1").arg(code);
 125             description = QObject::tr("Band %1").arg(code);
 126           }
 127           bandName = QString("BAND%1_FILE_NAME").arg(bandNumber);
 128           float lMin = minMaxRadiance[QString("LMIN_BAND%1").arg(bandNumber)].toDouble();
 129           float lMax = minMaxRadiance[QString("LMAX_BAND%1").arg(bandNumber)].toDouble();
 130           float qCalMin = minMaxPixelValue[QString("QCALMIN_BAND%1").arg(bandNumber)].toDouble();
 131           float qCalMax = minMaxPixelValue[QString("QCALMAX_BAND%1").arg(bandNumber)].toDouble();
 132           float offset = lMin;
 133           float gain = (lMax - lMin) / (qCalMax - qCalMin);
 134           tcChannelFile* channel = new tcChannelFile(path + productMetadata[bandName].toString(),
 135                                                      channelNames[band], description,
 136                                                      offset, gain);
 137           channelManager()->addChannel(channel);
 138           // Obtain the transformation from the first of the files.
 139           if (0 == band)
 140           {
 141             setTexToGeo(channel->texToGeo());
 142           }
 143           last = code;
 144         }
 145       }
 146       // Landsat 5 data is a little different
 147       else
 148       {
 149         tcLandsatMetaData::Reference inventoryMetaData = meta("INVENTORYMETADATA");
 150         tcLandsatMetaData::Reference sensorCharacteristics = inventoryMetaData("SENSORCHARACTERISTIC")("SENSORCHARACTERISTICCONTAINER",0);
 151         tcLandsatMetaData::Reference singleDateTime = inventoryMetaData("SINGLEDATETIME");
 152
 153         setName(QObject::tr("%1 %2 (%3)").arg(sensorCharacteristics("PLATFORMSHORTNAME",0)["VALUE"].toString())
 154                                          .arg(sensorCharacteristics("SENSORSHORTNAME",0)["VALUE"].toString())
 155                                          .arg(singleDateTime("CALENDARDATE",0)["VALUE"].toString()));
 156
 157         // Load additional attributes into a nice hash
 158         tcLandsatMetaData::Reference additionalAttributes = inventoryMetaData("ADDITIONALATTRIBUTES");
 159         QHash<QString, QVariant> attributes;
 160         int numAttributes = additionalAttributes.hasObject("ADDITIONALATTRIBUTESCONTAINER");
 161         for (int i = 0; i < numAttributes; ++i)
 162         {
 163           tcLandsatMetaData::Reference container = additionalAttributes("ADDITIONALATTRIBUTESCONTAINER", i);
 164           QString name = container("ADDITIONALATTRIBUTENAME",0)["VALUE"].toString();
 165           attributes[name] = container("INFORMATIONCONTENT")("PARAMETERVALUE", 0)["VALUE"];
 166         }
 167
 168         // Lon/Lat coordinates
 169         tcLandsatMetaData::Reference gRingPoint = inventoryMetaData("SPATIALDOMAINCONTAINER")
 170                                                                    ("HORIZONTALSPATIALDOMAINCONTAINER")
 171                                                                    ("GPOLYGON")
 172                                                                    ("GPOLYGONCONTAINER",0)
 173                                                                    ("GRINGPOINT");
 174         const QList<QVariant>& lons = gRingPoint("GRINGPOINTLONGITUDE",0)["VALUE"].toList();
 175         const QList<QVariant>& lats = gRingPoint("GRINGPOINTLATITUDE",0)["VALUE"].toList();
 176
 177         m_coordinates[0][0] = tcGeo(lons[0].toDouble() * M_PI/180.0,
 178                                     lats[0].toDouble() * M_PI/180.0);
 179         m_coordinates[1][0] = tcGeo(lons[1].toDouble() * M_PI/180.0,
 180                                     lats[1].toDouble() * M_PI/180.0);
 181         m_coordinates[0][1] = tcGeo(lons[3].toDouble() * M_PI/180.0,
 182                                     lats[3].toDouble() * M_PI/180.0);
 183         m_coordinates[1][1] = tcGeo(lons[2].toDouble() * M_PI/180.0,
 184                                     lats[2].toDouble() * M_PI/180.0);
 185
 186         // Sun direction
 187         setSunDirection(tcGeo(M_PI - M_PI/180.0 * attributes["SolarAzimuth"].toDouble(),
 188                                      M_PI/180.0 * attributes["SolarElevation"].toDouble()),
 189                         m_coordinates[1][0]);
 190
 191         // Names and descriptions of channels
 192         QStringList channelNames;
 193         channelNames << QObject::tr("Blue Green");
 194         channelNames << QObject::tr("Green");
 195         channelNames << QObject::tr("Red");
 196         channelNames << QObject::tr("Near IR");
 197         channelNames << QObject::tr("Mid IR 1");
 198         channelNames << QObject::tr("Thermal IR");
 199         channelNames << QObject::tr("Mid IR 2");
 200
 201         // Get input files and filter tif files
 202         int band = 0;
 203         foreach (QVariant file, inventoryMetaData("INPUTGRANULE")("INPUTPOINTER",0)["VALUE"].toList())
 204         {
 205           static QRegExp tifExt("^.*\\.tif$");
 206           if (tifExt.exactMatch(file.toString()))
 207           {
 208             float bias = attributes[QString("Band%1GainSetting").arg(band+1)].toDouble();
 209             float gain = attributes[QString("Band%1BiasSetting").arg(band+1)].toDouble();
 210             if (gain == 0.0f)
 211             {
 212               gain = 1.0f;
 213             }
 214
 215             tcChannelFile* channel = new tcChannelFile(path + file.toString(),
 216                                                        channelNames[band], QString("Band %1").arg(band),
 217                                                        bias, gain);
 218             channelManager()->addChannel(channel);
 219             // Obtain the transformation from the first of the files.
 220             if (0 == band)
 221             {
 222               setTexToGeo(channel->texToGeo());
 223             }
 224             ++band;
 225           }
 226         }
 227       }
 228
 229
 230
 231       QList<tcChannel*> mainColourChannels;
 232       for (int i = 0; i < channelManager()->numChannels(); ++i)
 233       {
 234         mainColourChannels += channelManager()->channel(i);
 235       }
 236       tcNegativeProduct* shadowClassification = new tcNegativeProduct(mainColourChannels);
 237       channelManager()->addChannel(shadowClassification);
 238       setShadowClassification(shadowClassification);
 239
 240       tcShadowDepth* shadowDepth1 = new tcShadowDepth(shadowClassification, dem, this, false);
 241       channelManager()->addChannel(shadowDepth1);
 242       tcShadowDepth* shadowDepth2 = new tcShadowDepth(shadowClassification, dem, this, true);
 243       channelManager()->addChannel(shadowDepth2);
 244
 245       tcLambertianShading* diffuse = new tcLambertianShading(channelManager()->channel(0), dem, this);
 246       channelManager()->addChannel(diffuse);
 247
 248       tcRaytracedShadowMap* raytrace = new tcRaytracedShadowMap(channelManager()->channel(0), dem, this);
 249       channelManager()->addChannel(raytrace);
 250
 251       tcChannel* shading = new tcPixelOp(tcPixelOp::Add, QList<tcChannel*>()
 252         << channelManager()->channel(0) // b
 253         << channelManager()->channel(1) // g
 254         << channelManager()->channel(2) // r
 255         << channelManager()->channel(3) // nir
 256         //<< channelManager()->channel(4) // mir
 257         //<< channelManager()->channel(5) // tir
 258         //<< channelManager()->channel(6) // tir
 259         //<< channelManager()->channel(channelManager()->numChannels()-2) // mir
 260         //<< channelManager()->channel(channelManager()->numChannels()-1) // pan
 261         //, 1.0f / 800
 262         , false
 263         );
 264       channelManager()->addChannel(shading);
 265       setShading(shading);
 266
 267       /*
 268       QList<tcChannel*> chromaticities;
 269       for (int i = 0; i <= 8; ++i)
 270       {
 271         chromaticities += new tcChannelChromaticity(channelManager()->channel(0), channelManager()->channel(i));
 272         chromaticities += new tcChannelChromaticity(channelManager()->channel(1), channelManager()->channel(i));
 273         chromaticities += new tcChannelChromaticity(channelManager()->channel(2), channelManager()->channel(i));
 274         chromaticities += new tcChannelChromaticity(channelManager()->channel(3), channelManager()->channel(i));
 275         chromaticities += new tcChannelChromaticity(channelManager()->channel(4), channelManager()->channel(i));
 276         chromaticities += new tcChannelChromaticity(channelManager()->channel(5), channelManager()->channel(i));
 277         chromaticities += new tcChannelChromaticity(channelManager()->channel(6), channelManager()->channel(i));
 278         chromaticities += new tcChannelChromaticity(channelManager()->channel(7), channelManager()->channel(i));
 279         chromaticities += new tcChannelChromaticity(channelManager()->channel(8), channelManager()->channel(i));
 280       }
 281       channelManager()->addChannels(chromaticities);
 282       */
 283
 284       /*
 285       tcShadowFreeChromaticities* sfc = new tcShadowFreeChromaticities(chromaticities);
 286       channelManager()->addChannels(sfc->channels());
 287       tcIlluminantDiscontinuity* idm = new tcIlluminantDiscontinuity(chromaticities);
 288       channelManager()->addChannels(idm->channels());
 289       tcShadowClassification* sc = new tcShadowClassification(idm, channelManager()->channel(4));
 290       channelManager()->addChannel(sc);
 291       */
 292
 293       shadowClassification->loadParams(path + "negativeproduct.npp");
 294
 295       // Only load one metafile
 296       break;
 297     }
 298   }
 299 }
 300
 301 /// Destructor.
 302 tcLandsatData::~tcLandsatData()
 303 {
 304 }