pktools 2.6.7
Processing Kernel for geospatial data
pkdsm2shadow.cc
1/**********************************************************************
2pkdsm2shadow.cc: program to calculate sun shadow based on digital surface model and sun angles
3Copyright (C) 2008-2014 Pieter Kempeneers
4
5This file is part of pktools
6
7pktools is free software: you can redistribute it and/or modify
8it under the terms of the GNU General Public License as published by
9the Free Software Foundation, either version 3 of the License, or
10(at your option) any later version.
11
12pktools is distributed in the hope that it will be useful,
13but WITHOUT ANY WARRANTY; without even the implied warranty of
14MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15GNU General Public License for more details.
16
17You should have received a copy of the GNU General Public License
18along 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/******************************************************************************/
79using namespace std;
80
81/*------------------
82 Main procedure
83 ----------------*/
84int 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
180 filter2d::Filter2d filter2d;
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}
int nrOfRow(void) const
Get the number of rows of this dataset.
std::string getProjection(void) const
Get the projection string (deprecated, use getProjectionRef instead)
std::string getGeoTransform() const
Get the geotransform data for this dataset as a string.
GDALColorTable * getColorTable(int band=0) const
Get the GDAL color table for this dataset as an instance of the GDALColorTable class.
CPLErr setGeoTransform(double *gt)
Set the geotransform data for this dataset.
int nrOfBand(void) const
Get the number of bands of this dataset.
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
double getDeltaX(void) const
Get the pixel cell spacing in x.
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
std::string getInterleave() const
Get the band coding (interleave)
GDALDataType getDataType(int band=0) const
Get the GDAL datatype for this dataset.
CPLErr setProjection(const std::string &projection)
Set the projection for this dataset in well known text (wkt) format.
void close(void)
Set the memory (in MB) to cache a number of rows in memory.
void open(const std::string &filename, const GDALAccess &readMode=GA_ReadOnly)
Open an image.
void close(void)
Close the image.
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...
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)