pktools  2.6.7
Processing Kernel for geospatial data
pkdsm2shadow.cc
1 /**********************************************************************
2 pkdsm2shadow.cc: program to calculate sun shadow based on digital surface model and sun angles
3 Copyright (C) 2008-2014 Pieter Kempeneers
4 
5 This file is part of pktools
6 
7 pktools is free software: you can redistribute it and/or modify
8 it under the terms of the GNU General Public License as published by
9 the Free Software Foundation, either version 3 of the License, or
10 (at your option) any later version.
11 
12 pktools is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 GNU General Public License for more details.
16 
17 You should have received a copy of the GNU General Public License
18 along with pktools. If not, see <http://www.gnu.org/licenses/>.
19 ***********************************************************************/
20 #include <assert.h>
21 #include <iostream>
22 #include <string>
23 #include <fstream>
24 #include <math.h>
25 #include <sys/types.h>
26 #include <stdio.h>
27 #include "base/Optionpk.h"
28 #include "base/Vector2d.h"
29 #include "algorithms/Filter2d.h"
30 #include "imageclasses/ImgReaderGdal.h"
31 #include "imageclasses/ImgWriterGdal.h"
32 
33 /******************************************************************************/
79 using namespace std;
80 
81 /*------------------
82  Main procedure
83  ----------------*/
84 int main(int argc,char **argv) {
85  Optionpk<std::string> input_opt("i","input","input image file");
86  Optionpk<std::string> output_opt("o", "output", "Output image file");
87  Optionpk<double> sza_opt("sza", "sza", "Sun zenith angle.");
88  Optionpk<double> saa_opt("saa", "saa", "Sun azimuth angle (N=0 E=90 S=180 W=270).");
89  Optionpk<int> flag_opt("f", "flag", "Flag to put in image if pixel shadow", 0);
90  Optionpk<std::string> otype_opt("ot", "otype", "Data type for output image ({Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/CInt16/CInt32/CFloat32/CFloat64}). Empty string: inherit type from input image", "");
91  Optionpk<string> oformat_opt("of", "oformat", "Output image format (see also gdal_translate).","GTiff");
92  Optionpk<string> colorTable_opt("ct", "ct", "color table (file with 5 columns: id R G B ALFA (0: transparent, 255: solid)");
93  Optionpk<string> option_opt("co", "co", "Creation option for output file. Multiple options can be specified.");
94  Optionpk<double> scale_opt("s", "scale", "scale used for input dsm: height=scale*input+offset");
95  Optionpk<double> offset_opt("off", "offset", "offset used for input dsm: height=scale*input+offset");
96  Optionpk<short> verbose_opt("v", "verbose", "verbose mode if > 0", 0,2);
97 
98  scale_opt.setHide(1);
99  offset_opt.setHide(1);
100 
101  bool doProcess;//stop process when program was invoked with help option (-h --help)
102  try{
103  doProcess=input_opt.retrieveOption(argc,argv);
104  output_opt.retrieveOption(argc,argv);
105  sza_opt.retrieveOption(argc,argv);
106  saa_opt.retrieveOption(argc,argv);
107  flag_opt.retrieveOption(argc,argv);
108  scale_opt.retrieveOption(argc,argv);
109  offset_opt.retrieveOption(argc,argv);
110  option_opt.retrieveOption(argc,argv);
111  otype_opt.retrieveOption(argc,argv);
112  oformat_opt.retrieveOption(argc,argv);
113  colorTable_opt.retrieveOption(argc,argv);
114  verbose_opt.retrieveOption(argc,argv);
115  }
116  catch(string predefinedString){
117  std::cout << predefinedString << std::endl;
118  exit(0);
119  }
120  if(!doProcess){
121  cout << endl;
122  cout << "Usage: pkdsm2shadow -i input.txt -o output [-sza angle] [-saa angle]" << endl;
123  cout << endl;
124  std::cout << "short option -h shows basic options only, use long option --help to show all options" << std::endl;
125  exit(0);//help was invoked, stop processing
126  }
127 
128  ImgReaderGdal input;
129  ImgWriterGdal output;
130  assert(input_opt.size());
131  assert(output_opt.size());
132  input.open(input_opt[0]);
133  if(scale_opt.size())
134  input.setScale(scale_opt[0]);
135  if(offset_opt.size())
136  input.setOffset(offset_opt[0]);
137 
138  // output.open(output_opt[0],input);
139  GDALDataType theType=GDT_Unknown;
140  if(verbose_opt[0])
141  cout << "possible output data types: ";
142  for(int iType = 0; iType < GDT_TypeCount; ++iType){
143  if(verbose_opt[0])
144  cout << " " << GDALGetDataTypeName((GDALDataType)iType);
145  if( GDALGetDataTypeName((GDALDataType)iType) != NULL
146  && EQUAL(GDALGetDataTypeName((GDALDataType)iType),
147  otype_opt[0].c_str()))
148  theType=(GDALDataType) iType;
149  }
150  if(theType==GDT_Unknown)
151  theType=input.getDataType();
152 
153  if(verbose_opt[0])
154  std::cout << std::endl << "Output pixel type: " << GDALGetDataTypeName(theType) << endl;
155 
156  string imageType;//=input.getImageType();
157  if(oformat_opt.size())
158  imageType=oformat_opt[0];
159 
160  if(option_opt.findSubstring("INTERLEAVE=")==option_opt.end()){
161  string theInterleave="INTERLEAVE=";
162  theInterleave+=input.getInterleave();
163  option_opt.push_back(theInterleave);
164  }
165  try{
166  output.open(output_opt[0],input.nrOfCol(),input.nrOfRow(),input.nrOfBand(),theType,imageType,option_opt);
167  }
168  catch(string errorstring){
169  cout << errorstring << endl;
170  exit(4);
171  }
172  output.setProjection(input.getProjection());
173  double gt[6];
174  input.getGeoTransform(gt);
175  output.setGeoTransform(gt);
176 
177  if(input.getColorTable()!=NULL)
178  output.setColorTable(input.getColorTable());
179 
181  if(verbose_opt[0])
182  std::cout<< "class values: ";
183  if(colorTable_opt.size())
184  output.setColorTable(colorTable_opt[0]);
185  filter2d.shadowDsm(input,output,sza_opt[0],saa_opt[0],input.getDeltaX(),flag_opt[0]);
186  input.close();
187  output.close();
188  return 0;
189 }
void setOffset(double theOffset, int band=0)
Set offset for a specific band when writing the raster data values. The scaling and offset are applie...
Definition: ImgRasterGdal.h:84
GDALColorTable * getColorTable(int band=0) const
Get the GDAL color table for this dataset as an instance of the GDALColorTable class.
int nrOfBand(void) const
Get the number of bands of this dataset.
STL namespace.
GDALDataType getDataType(int band=0) const
Get the GDAL datatype for this dataset.
void close(void)
Set the memory (in MB) to cache a number of rows in memory.
int nrOfRow(void) const
Get the number of rows of this dataset.
void setColorTable(const std::string &filename, int band=0)
Set the color table using an (ASCII) file with 5 columns (value R G B alpha)
CPLErr setGeoTransform(double *gt)
Set the geotransform data for this dataset.
std::string getGeoTransform() const
Get the geotransform data for this dataset as a string.
void open(const std::string &filename, const ImgReaderGdal &imgSrc, const std::vector< std::string > &options=std::vector< std::string >())
Open an image for writing, copying image attributes from a source image. Image is directly written to...
double getDeltaX(void) const
Get the pixel cell spacing in x.
CPLErr setProjection(const std::string &projection)
Set the projection for this dataset in well known text (wkt) format.
std::string getInterleave() const
Get the band coding (interleave)
void close(void)
Close the image.
std::string getProjection(void) const
Get the projection string (deprecated, use getProjectionRef instead)
void open(const std::string &filename, const GDALAccess &readMode=GA_ReadOnly)
Open an image.
void setScale(double theScale, int band=0)
Set scale for a specific band when writing the raster data values. The scaling and offset are applied...
Definition: ImgRasterGdal.h:75
int nrOfCol(void) const
Get the number of columns of this dataset.
Definition: ImgRasterGdal.h:98