1
2 ///////////////////////////////////////////////////////////
3 // //
4 // SAGA //
5 // //
6 // System for Automated Geoscientific Analyses //
7 // //
8 // Tool Library //
9 // //
10 // io_gdal //
11 // //
12 //-------------------------------------------------------//
13 // //
14 // gdal_import_aster.cpp //
15 // //
16 // Copyright (C) 2018 O. Conrad //
17 // //
18 //-------------------------------------------------------//
19 // //
20 // This file is part of 'SAGA - System for Automated //
21 // Geoscientific Analyses'. SAGA is free software; you //
22 // can redistribute it and/or modify it under the terms //
23 // of the GNU General Public License as published by the //
24 // Free Software Foundation, either version 2 of the //
25 // License, or (at your option) any later version. //
26 // //
27 // SAGA is distributed in the hope that it will be //
28 // useful, but WITHOUT ANY WARRANTY; without even the //
29 // implied warranty of MERCHANTABILITY or FITNESS FOR A //
30 // PARTICULAR PURPOSE. See the GNU General Public //
31 // License for more details. //
32 // //
33 // You should have received a copy of the GNU General //
34 // Public License along with this program; if not, see //
35 // <http://www.gnu.org/licenses/>. //
36 // //
37 //-------------------------------------------------------//
38 // //
39 // e-mail: oconrad@saga-gis.de //
40 // //
41 // contact: Olaf Conrad //
42 // Institute of Geography //
43 // University of Hamburg //
44 // Germany //
45 // //
46 ///////////////////////////////////////////////////////////
47
48 //---------------------------------------------------------
49 #include "gdal_import_aster.h"
50
51
52 ///////////////////////////////////////////////////////////
53 // //
54 // //
55 // //
56 ///////////////////////////////////////////////////////////
57
58 //---------------------------------------------------------
CGDAL_Import_ASTER(void)59 CGDAL_Import_ASTER::CGDAL_Import_ASTER(void)
60 {
61 //-----------------------------------------------------
62 Set_Name (_TL("Import ASTER Scene"));
63
64 Set_Author ("O.Conrad (c) 2018");
65
66 Set_Description (_TW(
67 "Import ASTER scene from Hierarchical Data Format version 4 (HDF4). "
68 ));
69
70 Add_Reference("https://asterweb.jpl.nasa.gov/", _TL("ASTER Homepage Jet Propulsion Laboratory"));
71 Add_Reference("https://lpdaac.usgs.gov/dataset_discovery/aster", _TL("ASTER Overview at NASA/USGS"));
72
73 //-----------------------------------------------------
74 Parameters.Add_FilePath("",
75 "FILE" , _TL("File"),
76 _TL(""),
77 CSG_String::Format("%s (*.hdf)|*.hdf|%s|*.*",
78 _TL("HDF4 Files"),
79 _TL("All Files")
80 ), NULL, false
81 );
82
83 Parameters.Add_Choice("",
84 "FORMAT" , _TL("Format"),
85 _TL(""),
86 CSG_String::Format("%s|%s|",
87 _TL("single grids"),
88 _TL("grid collections")
89 ), 1
90 );
91
92 Parameters.Add_Grids_Output("", "VNIR", _TL("Visible and Near Infrared"), _TL(""));
93 Parameters.Add_Grids_Output("", "SWIR", _TL("Short Wave Infrared" ), _TL(""));
94 Parameters.Add_Grids_Output("", "TIR", _TL("Thermal Infrared" ), _TL(""));
95
96 Parameters.Add_Grid_List("", "BANDS", _TL("Bands"), _TL(""), PARAMETER_OUTPUT, false);
97
98 Parameters.Add_Table("", "METADATA" , _TL("Metadata"), _TL(""), PARAMETER_OUTPUT_OPTIONAL);
99 }
100
101
102 ///////////////////////////////////////////////////////////
103 // //
104 ///////////////////////////////////////////////////////////
105
106 //---------------------------------------------------------
On_Execute(void)107 bool CGDAL_Import_ASTER::On_Execute(void)
108 {
109 //-----------------------------------------------------
110 const double Wave[14][2] =
111 {
112 { 0.520, 0.600 },
113 { 0.630, 0.690 },
114 { 0.760, 0.860 },
115 { 1.600, 1.700 },
116 { 2.145, 2.185 },
117 { 2.185, 2.225 },
118 { 2.235, 2.285 },
119 { 2.295, 2.365 },
120 { 2.360, 2.430 },
121 { 8.125, 8.475 },
122 { 8.475, 8.825 },
123 { 8.925, 9.275 },
124 { 10.250, 10.950 },
125 { 10.950, 11.650 }
126 };
127
128 //-----------------------------------------------------
129 CSG_GDAL_DataSet DataSet;
130
131 if( !DataSet.Open_Read(Parameters("FILE")->asString()) )
132 {
133 Error_Fmt("%s [%s]", _TL("could not open file"), Parameters("FILE")->asString());
134
135 return( false );
136 }
137
138 CSG_String FileName = SG_File_Get_Name(Parameters("FILE")->asString(), false);
139
140 if( DataSet.Get_DriverID().Cmp("HDF4") )
141 {
142 Message_Fmt("\n%s: %s [%s]\n", _TL("Warning"), _TL("Driver"), DataSet.Get_DriverID().c_str());
143 }
144
145 //-----------------------------------------------------
146 CSG_MetaData MetaData;
147
148 if( !DataSet.Get_MetaData(MetaData) )
149 {
150 Error_Fmt("%s [%s]", _TL("failed to read metadata"), FileName.c_str());
151
152 return( false );
153 }
154 else if( Parameters("METADATA")->asTable() )
155 {
156 CSG_Table *pTable = Parameters("METADATA")->asTable();
157
158 pTable->Destroy();
159 pTable->Fmt_Name("%s [%s]", FileName.c_str(), _TL("Metadata"));
160 pTable->Add_Field("KEY" , SG_DATATYPE_String);
161 pTable->Add_Field("VALUE", SG_DATATYPE_String);
162
163 for(int i=0; i<MetaData.Get_Children_Count(); i++)
164 {
165 CSG_Table_Record *pRecord = pTable->Add_Record();
166
167 pRecord->Set_Value(0, MetaData[i].Get_Name ());
168 pRecord->Set_Value(1, MetaData[i].Get_Content());
169 }
170 }
171
172 //-----------------------------------------------------
173 TSG_Rect Extent; CSG_Projection Projection;
174
175 if( !Get_System(MetaData, Extent, Projection) )
176 {
177 Error_Fmt("%s [%s]", _TL("failed to read spatial reference"), FileName.c_str());
178
179 return( false );
180 }
181
182 //-----------------------------------------------------
183 CSG_Strings SubDataSets = DataSet.Get_SubDataSets();
184
185 if( SubDataSets.Get_Count() < 1 )
186 {
187 Error_Fmt("%s [%s]", _TL("failed to read subdatasets"), FileName.c_str());
188
189 return( false );
190 }
191
192 //-----------------------------------------------------
193 CSG_Parameter_Grid_List *pBands = NULL; CSG_Grids *pVNIR, *pSWIR, *pTIR; CSG_Table Attributes;
194
195 if( Parameters("FORMAT")->asInt() == 0 )
196 {
197 pBands = Parameters("BANDS")->asGridList();
198
199 pBands->Del_Items();
200 }
201 else
202 {
203 Attributes.Add_Field("ID" , SG_DATATYPE_Char );
204 Attributes.Add_Field("NAME" , SG_DATATYPE_String);
205 Attributes.Add_Field("WAVMIN", SG_DATATYPE_Double);
206 Attributes.Add_Field("WAVMID", SG_DATATYPE_Double);
207 Attributes.Add_Field("WAVMAX", SG_DATATYPE_Double);
208 Attributes.Add_Record();
209
210 Parameters("VNIR")->Set_Value(pVNIR = SG_Create_Grids());
211 pVNIR->Get_Attributes_Ptr()->Create(&Attributes);
212 pVNIR->Fmt_Name("%s %s", FileName.c_str(), SG_T("VNIR"));
213 pVNIR->Get_MetaData().Add_Child(MetaData)->Set_Name("ASTER");
214 pVNIR->Set_Z_Attribute (3);
215 pVNIR->Set_Z_Name_Field(1);
216
217 Parameters("SWIR")->Set_Value(pSWIR = SG_Create_Grids());
218 pSWIR->Get_Attributes_Ptr()->Create(&Attributes);
219 pSWIR->Fmt_Name("%s %s", FileName.c_str(), SG_T("SWIR"));
220 pSWIR->Get_MetaData().Add_Child(MetaData)->Set_Name("ASTER");
221 pSWIR->Set_Z_Attribute (3);
222 pSWIR->Set_Z_Name_Field(1);
223
224 Parameters( "TIR")->Set_Value(pTIR = SG_Create_Grids());
225 pTIR ->Get_Attributes_Ptr()->Create(&Attributes);
226 pTIR ->Fmt_Name("%s %s", FileName.c_str(), SG_T( "TIR"));
227 pTIR ->Get_MetaData().Add_Child(MetaData)->Set_Name("ASTER");
228 pTIR ->Set_Z_Attribute (3);
229 pTIR ->Set_Z_Name_Field(1);
230 }
231
232 //-----------------------------------------------------
233 for(int i=0; i<SubDataSets.Get_Count() && Process_Get_Okay(); i++)
234 {
235 CSG_Grid *pDataSet = SubDataSets[i].Find("EOS_SWATH") > 0 && DataSet.Open_Read(SubDataSets[i]) ? DataSet.Read(0) : NULL;
236
237 if( pDataSet )
238 {
239 CSG_Grid *pBand = SG_Create_Grid(
240 pDataSet->Get_Type(), pDataSet->Get_NX(), pDataSet->Get_NY(),
241 (Extent.xMax - Extent.xMin) / (pDataSet->Get_NX() - 1),
242 Extent.xMin, Extent.yMin
243 );
244
245 if( pBand )
246 {
247 CSG_String Band = SubDataSets[i].AfterLast(':'); Band.Replace("ImageData", "");
248 int iBand = Band.asInt();
249
250 Band.Printf("Band %02d %s", iBand, iBand <= 3 ? SG_T("VNIR") : iBand <= 9 ? SG_T("SWIR") : SG_T("TIR"));
251
252 pBand->Set_Name(Band);
253 pBand->Get_Projection().Create(Projection);
254 pBand->Set_NoData_Value(0.0);
255 pBand->Set_Description(pDataSet->Get_MetaData().asText());
256
257 #pragma omp parallel for
258 for(int y=0; y<pDataSet->Get_NY(); y++)
259 {
260 for(int x=0; x<pDataSet->Get_NX(); x++)
261 {
262 pBand->Set_Value(x, y, pDataSet->asDouble(x, y));
263 }
264 }
265
266 if( pBands )
267 {
268 Parameters("BANDS")->asGridList()->Add_Item(pBand);
269 }
270 else
271 {
272 Attributes[0].Set_Value(0, iBand);
273 Attributes[0].Set_Value(1, Band);
274 Attributes[0].Set_Value(2, Wave[iBand - 1][0]);
275 Attributes[0].Set_Value(3, (Wave[iBand - 1][0] + Wave[iBand - 1][1]) / 2.0);
276 Attributes[0].Set_Value(4, Wave[iBand - 1][1]);
277
278 if( iBand <= 3 )
279 {
280 pVNIR->Add_Grid(Attributes[0], pBand, true);
281 }
282 else if( iBand <= 9 )
283 {
284 pSWIR->Add_Grid(Attributes[0], pBand, true);
285 }
286 else // if( iBand <= 14 )
287 {
288 pTIR ->Add_Grid(Attributes[0], pBand, true);
289 }
290 }
291 }
292
293 delete(pDataSet);
294 }
295 }
296
297 //-----------------------------------------------------
298 return( true );
299 }
300
301
302 ///////////////////////////////////////////////////////////
303 // //
304 ///////////////////////////////////////////////////////////
305
306 //---------------------------------------------------------
Get_System(const CSG_MetaData & MetaData,TSG_Rect & Extent,CSG_Projection & Projection)307 bool CGDAL_Import_ASTER::Get_System(const CSG_MetaData &MetaData, TSG_Rect &Extent, CSG_Projection &Projection)
308 {
309 if( !MetaData("UPPERLEFTM")
310 || !MetaData("LOWERRIGHTM")
311 || !MetaData("UTMZONENUMBER")
312 || !MetaData("NORTHBOUNDINGCOORDINATE") )
313 {
314 return( false );
315 }
316
317 Extent.xMin = MetaData["UPPERLEFTM" ].Get_Content().AfterFirst (',').asDouble();
318 Extent.yMax = MetaData["UPPERLEFTM" ].Get_Content().BeforeFirst(',').asDouble();
319 Extent.xMax = MetaData["LOWERRIGHTM" ].Get_Content().AfterFirst (',').asDouble();
320 Extent.yMin = MetaData["LOWERRIGHTM" ].Get_Content().BeforeFirst(',').asDouble();
321
322 int Zone = MetaData["UTMZONENUMBER"].Get_Content().asInt(), EPSG_ID;
323
324 bool bSouth = MetaData["NORTHBOUNDINGCOORDINATE"].Get_Content().asDouble() < 0.0;
325
326 if( bSouth )
327 {
328 Extent.yMin += 10000000;
329 Extent.yMax += 10000000;
330
331 EPSG_ID = 32700 + Zone;
332 }
333 else
334 {
335 EPSG_ID = 32600 + Zone;
336 }
337
338 if( !Projection.Create(EPSG_ID) )
339 {
340 CSG_String Proj4;
341
342 Proj4.Printf("+proj=utm +zone=%d +datum=WGS84 +units=m +no_defs ");
343
344 if( bSouth )
345 {
346 Proj4 += "+south ";
347 }
348
349 Projection.Create(Proj4, SG_PROJ_FMT_Proj4);
350 }
351
352 return( true );
353 }
354
355
356 ///////////////////////////////////////////////////////////
357 // //
358 // //
359 // //
360 ///////////////////////////////////////////////////////////
361
362 //---------------------------------------------------------
363