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