1 /*
2 The MIT License (MIT)
3 
4 Copyright (C) 2015-2020 Shashank Khare, skhare at hotmail dot com
5 Copyright (C) 2015-2020 Philip Nienhuis <prnienhuis@users.sf.net>
6 Large parts of the below code originate from http://www.gdal.org/gdal_tutorial.html
7 (nowadays https://gdal.org/tutorials/raster_api_tut.html)
8 
9 Acknowledgement: Snow and Avalanche Studies Establishment, DRDO, Chandigarh, India
10 
11 Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated
12 documentation files (the "Software"), to deal in the Software without restriction, including without limitation the
13 rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit
14 persons to whom the Software is furnished to do so, subject to the following conditions:
15 
16 The above copyright notice and this permission notice shall be included in all copies or substantial portions of the
17 Software.
18 
19 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE
20 WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
21 COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
22 OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
23 
24 */
25 
26 #ifdef HAVE_CONFIG_H
27 #include "./config.h"
28 #endif
29 
30 #include <octave/oct.h>
31 #include <octave/ov-struct.h>
32 #include <octave/oct-map.h>
33 #include <octave/error.h>
34 
35 #include <cstring>
36 #include "misc.h"
37 #include "stdio.h"
38 
39 #ifdef HAVE_GDAL
40 #include <gdal_priv.h>
41 #include <cpl_conv.h>
42 #endif
43 
44 
45 /*
46 INPUT : string filepath
47 RETURNS: int containing return code,
48          struct array containing band data,
49          struct array containing metadata
50 */
51 
52 /*
53 See code examples here [1]:
54 http://www.gdal.org/gdal_tutorial.html
55 */
56 
57 DEFUN_DLD (gdalread, args, ,
58            "-*- texinfo -*-\n\
59 @deftypefn  {} {[@var{rstat}, @var{rinfo}, @var{bands}] =} gdalread (@var{fname})\n\
60 @deftypefnx {} {[@var{rstat}, @var{rinfo}, @var{bands}] =} gdalread (@var{fname}, @var{info})\n\
61 Get raster data and info or (optionally) just info from a GIS raster file.\n\
62 For normal raster data reading it is better to use rasterread.m as that \n\
63 takes care of postprocessing gdalread output.\n\
64 \n\
65 @var{fname} is the name (or full path name) of a raster file to be read.\n\
66 \n\
67 @var{info} is any non-null value.  If present, only raster info is\n\
68 returned.  This option can be useful if only raster information is\n\
69 required and reading the raster data can be skipped.\n\
70 \n\
71 @var{rstat} is set to zero if reading was successful, -1 otherwise.\n\
72 \n\
73 @var{rinfo} is a scalar struct containing information about the raster\n\
74 data: datatype_index (a numerical GDAL datatype index); datatype_name\n\
75 (GDAL type name); BitDepth; geotransformation (1x6 double vector);\n\
76 projection (a string containing projection information); Width (nr. of\n\
77 raster data columns); Height (nr. of raster data rows); FileType; and\n\
78 nbands (the number of raster bands in the file).\n\
79 \n\
80 @var{bands} is a struct array containing raster data for each band in\n\
81 the file. If only raster info was requested, bands is empty.\n\
82 When data is read, each band is a struct element containing: bbox\n\
83 (bounding box: [xmin ymin; xmax ymax]), data (a matrix with the actual\n\
84 raster data); min (minimum raster data value); max (maximum raster\n\
85 data value); has_ndv (indicates whether the band as a special value\n\
86 indicating pixels that shouldn't be processed); ndv_val (the special\n\
87 value for improper pixels), and colortable.  If the band has a\n\
88 colortable, field ``colortable'' contains an Mx4 array of colortable\n\
89 entries (see GDAL reference), otherwise this field is empty.\n\
90 Note that the actual raster data array (field ``data'' is rotated\n\
91 90 degrees clockwise.\n\
92 @seealso{rasterread,rot90}\n\
93 @end deftypefn")
94 
95 {
96 #ifndef HAVE_GDAL
97   error ("gdalread: reading of raster file with GDAL was disabled during installation");
98   return octave_value_list ();
99 #else
100 
101   int nargin = args.length ();
102   octave_value_list ret_list(3);
103   int num_items = 0;
104   int XSize = 0, YSize = 0;
105 
106   if (nargin < 1)
107   {
108     get_null_values (3, &ret_list);
109     print_usage ();
110     return ret_list;
111   };
112 
113   // [1] "Opening the file"
114   std::string arg_file = args(0).string_value ();
115   GDALDataset* poDS;
116   GDALAllRegister ();
117   poDS = (GDALDataset*) GDALOpen (arg_file.c_str(), GA_ReadOnly);
118   if (poDS == NULL)
119   {
120     octave_stdout << "Error: Open failed.\n";
121     get_null_values (3, &ret_list);
122     return ret_list;
123   }
124 
125   // [1] "Getting dataset information"
126   double adfGeoTransform[6];
127   int rasterX;                        // = poDS->GetRasterXSize();
128   int rasterY;                        // = poDS->GetRasterYSize();
129   octave_idx_type band_count = poDS->GetRasterCount ();
130 
131   if (!(poDS->GetGeoTransform(adfGeoTransform) == CE_None))
132   {
133     // error could not get geotransformation data
134     // put origin as 0,0 and pixel size as 1,1
135     // should not happen but if it does proceed with 0,0 as origin
136     octave_stdout << "Warning: GetGeoTransform failed.\n";
137     adfGeoTransform[0] = 0;
138     adfGeoTransform[3] = 0;
139     adfGeoTransform[1] = 1;
140     adfGeoTransform[5] = 1;
141   }
142 
143   /************************READING RASTER DATA **********************/
144   GDALRasterBand* poBand;
145   int raster_data_type, raster_type_size;
146   const char* raster_type_name;
147 
148   // Only get raster data when no info request was given
149   if (nargin <= 1)
150   {
151     static const char *fields[] = {"bbox", "data", "min", "max",
152                                    "has_ndv", "ndv_val", "colortable", 0};
153     // Output struct
154     octave_map m_band (dim_vector (band_count, 1), string_vector (fields));
155     octave_scalar_map band_struct = (string_vector (fields));
156 
157     // temp variables in below loop
158     int nBlockXSize, nBlockYSize;
159     int bGotMin, bGotMax;
160     double adfMinMax[2];
161     Matrix raster_data, raster_data_tmp, X_tmp, Y_tmp, colortable;
162     Matrix bbox(2, 2);
163 
164     int clrs;
165     GDALColorTable* GDALclr_table;
166     const GDALColorEntry* clr_entry;
167 
168     for (octave_idx_type curr_band = 0; curr_band < band_count; curr_band++)
169     {
170       // [1] "Fetch a raster band"
171       // Get actual raster data
172       poBand = poDS->GetRasterBand (curr_band + 1);
173       poBand->GetBlockSize (&nBlockXSize, &nBlockYSize);
174 
175       raster_data_type = (int)poBand->GetRasterDataType ();
176       raster_type_name = GDALGetDataTypeName (poBand->GetRasterDataType ());
177       raster_type_size = (int)GDALGetDataTypeSize (poBand->GetRasterDataType ());
178 
179       adfMinMax[0] = poBand->GetMinimum (&bGotMin);
180       adfMinMax[1] = poBand->GetMaximum (&bGotMax);
181       if (!(bGotMin && bGotMax))
182         GDALComputeRasterMinMax ((GDALRasterBandH)poBand, TRUE, adfMinMax);
183 
184       // Get colortable, if present
185       clrs = 0;
186       GDALclr_table = poBand->GetColorTable ();
187       if (GDALclr_table != NULL)
188       {
189         clrs = GDALclr_table->GetColorEntryCount ();
190         colortable = Matrix (dim_vector(clrs, 4));
191         for (octave_idx_type c_entry = 0; c_entry < clrs; c_entry++)
192         {
193           clr_entry = GDALclr_table->GetColorEntry (c_entry);
194           colortable(c_entry, 0) = (double) clr_entry->c1;
195           colortable(c_entry, 1) = (double) clr_entry->c2;
196           colortable(c_entry, 2) = (double) clr_entry->c3;
197           colortable(c_entry, 3) = (double) clr_entry->c4;
198         }
199       }
200       else
201         int colortable = -1;
202 
203       /******************* READING RASTER DATA *******************************************/
204       rasterX = poBand->GetXSize ();
205       rasterY = poBand->GetYSize ();
206       ///////////////////////////////////////////////////
207       raster_data = Matrix(dim_vector(rasterX, rasterY));
208 
209       long int nXBlocks = (rasterX + nBlockXSize - 1) / nBlockXSize;
210       long int nYBlocks = (rasterY + nBlockYSize - 1) / nBlockYSize;
211 
212       // read the data into memory
213       GByte* data_ptr = (GByte*) raster_data.fortran_vec ();
214       long int* pafScanline = (long int*) CPLMalloc (sizeof (double) * rasterX * rasterY);
215       long int read_size_x = 0, read_size_y = 0, nXOff, nYOff;
216       CPLErr read_err =
217         poBand->RasterIO (GF_Read, 0, 0, rasterX, rasterY, pafScanline, rasterX, rasterY, GDT_Float64, 0, 0);
218       if (read_err != CE_None)
219       {
220         ret_list(0) = -1;
221         return ret_list;
222       }
223       // TODO: memcpy is not needed
224       memcpy (data_ptr, pafScanline, sizeof (double) * rasterX * rasterY);
225 
226       /************ calculate bounds for the band  *********************/
227       double minx = adfGeoTransform[0];
228       double maxy = adfGeoTransform[3];
229       double maxx = minx + adfGeoTransform[1] * rasterX;
230       double miny = maxy + adfGeoTransform[5] * rasterY;
231       bbox(0, 0) = minx;
232       bbox(0, 1) = miny;
233       bbox(1, 0) = maxx;
234       bbox(1, 1) = maxy;
235 
236       band_struct.setfield ("bbox", bbox);
237       band_struct.setfield ("data", raster_data);
238       band_struct.setfield ("min",  adfMinMax[0]);
239       band_struct.setfield ("max",  adfMinMax[1]);
240 
241       int ndv_present = 0;
242       double ndvalue = poBand->GetNoDataValue(&ndv_present);
243       if (ndv_present)
244         band_struct.setfield ("has_ndv", ndv_present);
245       else
246         band_struct.setfield ("has_ndv", -1);
247       band_struct.setfield ("ndv_val", ndvalue);
248       band_struct.setfield ("colortable", colortable);
249 
250       // Assign temp struct to output struct
251       m_band.fast_elem_insert (curr_band, band_struct);
252 
253       CPLFree (pafScanline);
254     }
255     ret_list(2) = octave_value(m_band);
256   }
257   else
258   {
259     // Get raster info
260     poBand = poDS->GetRasterBand (1);
261     raster_data_type = (int)poBand->GetRasterDataType ();
262     raster_type_name = GDALGetDataTypeName (poBand->GetRasterDataType ());
263     raster_type_size = (int)GDALGetDataTypeSize (poBand->GetRasterDataType ());
264     rasterX = poBand->GetXSize ();
265     rasterY = poBand->GetYSize ();
266     ret_list(2) = -1;
267   }
268 
269   /*********************** generate meta data **************/
270   Matrix geom(1, 6);
271   geom(0, 0) = adfGeoTransform[0];
272   geom(0, 1) = adfGeoTransform[1];
273   geom(0, 2) = adfGeoTransform[2];
274   geom(0, 3) = adfGeoTransform[3];
275   geom(0, 4) = adfGeoTransform[4];
276   geom(0, 5) = adfGeoTransform[5];
277 
278   // Output info struct
279   octave_scalar_map meta_data;
280 
281   meta_data.assign(std::string("FileType"),
282     poDS->GetDriver()->GetMetadataItem (GDAL_DMD_LONGNAME));
283   meta_data.assign(std::string("datatype_index"), raster_data_type);
284   meta_data.assign(std::string("datatype_name"), raster_type_name);
285   meta_data.assign(std::string("BitDepth"), raster_type_size);
286   meta_data.assign(std::string("GeoTransformation"), geom);
287   meta_data.assign(std::string("Projection"), poDS->GetProjectionRef());
288   meta_data.assign(std::string("Width"), rasterX);
289   meta_data.assign(std::string("Height"), rasterY);
290   meta_data.assign(std::string("nbands"), band_count);
291 
292   ret_list(0) = octave_value(0);
293   ret_list(1) = octave_value(meta_data);
294 
295   GDALClose ((GDALDatasetH) poDS);
296 
297   return ret_list;
298 #endif
299 }
300