updated GDAl writing

impls/square-cell/vn-neighbours/1-levels/simple/ESRI_ASCIIGrid.hpp

@@ -61,9 +61,12 @@ namespace CA {
 
 		//};
 
-
-
-		
+#ifdef CAAPI_GDAL_OPTION
+		//static std::string driver = "AAIGrid";
+		static std::string driver;
+		//static int epsgCode = 27700;
+		static int epsgCode;
+#endif	
 
 
 
@@ -267,7 +270,7 @@ namespace CA {
 			AsciiGridGeneral<T>::nrows = dataSet->GetRasterYSize();
 			AsciiGridGeneral<T>::xllcorner = adfGeoTransform[0];
 			// check for that negative?
-			AsciiGridGeneral<T>::yllcorner = adfGeoTransform[3] + (dataSet->GetRasterYSize() * adfGeoTransform[5]);
+			AsciiGridGeneral<T>::yllcorner = adfGeoTransform[3] - (dataSet->GetRasterYSize() * abs(adfGeoTransform[5]));
 
 			// check for negative and that they are both the same value?
 			AsciiGridGeneral<T>::cellsize = adfGeoTransform[1];
@@ -304,6 +307,8 @@ namespace CA {
 			// Get starting time.
 			//Clock total_timer;
 
+#ifndef CAAPI_GDAL_OPTION
+
 #ifdef CA2D_MPI
 
 			//MPI_Comm_size(MPI_COMM_WORLD, &_commSize);
@@ -324,6 +329,7 @@ namespace CA {
 			String _filename = filename + ".asc";
 #endif // CA2D_MPI
 
+
 			// Open file
 			std::ofstream onfile(_filename.c_str());
 
@@ -365,9 +371,106 @@ namespace CA {
 
 			// Close the file.
 			onfile.close();
+			
 
-			//std::cout<<"Write run time taken (s) = " << total_timer.millisecond()/1000.0 << std::endl;    
-		}
+#else
+
+			GDALAllRegister();
+
+			//GDALDataset* dataSet;
+			//dataSet = (GDALDataset*)GDALOpen(filename.c_str(), GDALAccess::GA_Update);
+
+			GDALDriver* poDriver;
+			poDriver = GetGDALDriverManager()->GetDriverByName(ESRI_ASCIIGrid<T>::driver.c_str());
+
+		
+			std::string extension = ".";
+
+			extension += poDriver->GetMetadataItem(GDAL_DMD_EXTENSION);
+
+
+
+			// Add file extension
+			String _filename = filename + extension;
+
+			GDALDataset* dataSet;
+			char** papszOptions = NULL;
+			dataSet = poDriver->Create(_filename.c_str(), AsciiGridGeneral<T>::ncols, AsciiGridGeneral<T>::nrows, 1,
+				GDT_Float64,
+				papszOptions);
+
+
+			// ulx, xres, xskew, uly, yskew, yres 
+			double adfGeoTransform[6];
+			//dataSet->GetGeoTransform(adfGeoTransform);
+
+
+			
+			adfGeoTransform[1] = AsciiGridGeneral<T>::cellsize;
+			adfGeoTransform[5] = -AsciiGridGeneral<T>::cellsize;
+			
+			adfGeoTransform[0] = AsciiGridGeneral<T>::xllcorner;
+
+			adfGeoTransform[3] = AsciiGridGeneral<T>::yllcorner + (AsciiGridGeneral<T>::nrows* AsciiGridGeneral<T>::cellsize);
+
+
+			// check for negative and that they are both the same value?
+			//adfGeoTransform[1] = AsciiGridGeneral<T>::cellsize;
+			//adfGeoTransform[4] = AsciiGridGeneral<T>::cellsize;
+			adfGeoTransform[2] = 0;
+			adfGeoTransform[4] = 0;
+
+			dataSet->SetGeoTransform(adfGeoTransform);
+
+
+			GDALRasterBand* poBand;
+
+			// might be zero in c++?
+			poBand = dataSet->GetRasterBand(1);
+
+			poBand->SetNoDataValue(AsciiGridGeneral<T>::nodata);
+
+			// Allocate data and set the default value to nodata.
+			//AsciiGridGeneral<T>::data.clear();
+			//AsciiGridGeneral<T>::data.resize(AsciiGridGeneral<T>::ncols * AsciiGridGeneral<T>::nrows, AsciiGridGeneral<T>::nodata);
+
+
+			
+			int   nXSize = poBand->GetXSize();
+			int   nYSize = poBand->GetYSize();
+			//pafScanline = (CA::Real*)CPLMalloc(sizeof(CA::Real) * nXSize);
+
+			OGRSpatialReference oSRS;
+			//OSRImportFromEPSG(osr, AsciiGridGeneral<T>::epsgCode);
+			char* pszSRS_WKT = NULL;
+
+			//oSRS.SetWellKnownGeogCS("NAD27");
+			oSRS.importFromEPSG(ESRI_ASCIIGrid<T>::epsgCode);
+			oSRS.exportToWkt(&pszSRS_WKT);
+			dataSet->SetProjection(pszSRS_WKT);
+
+			free(pszSRS_WKT);
+
+
+			poBand->RasterIO(GDALRWFlag::GF_Write, 0, 0, nXSize, nYSize,
+				static_cast<void*>(AsciiGridGeneral<T>::data.data()), nXSize, nYSize, GDT_Float64,
+				0, 0);
+
+
+			//GDAlClose((GDALDataSet)dataSet);
+
+
+			poBand->FlushCache();
+
+			dataSet->FlushCache();
+
+			//delete(poBand);
+			delete(dataSet);
+			//delete(poDriver);
+
+#endif
+						//std::cout<<"Write run time taken (s) = " << total_timer.millisecond()/1000.0 << std::endl;    
+		}// end of write function
 
 
 
@@ -1094,10 +1197,17 @@ namespace CA {
 
 	};// end of class
 
-	
-
+#ifdef CAAPI_GDAL_OPTION
+	template<typename T>
+	std::string CA::ESRI_ASCIIGrid<T>::driver = "AAIGrid";
+	template<typename T>
+	int CA::ESRI_ASCIIGrid<T>::epsgCode = 27700;
+#endif
 
 }// end of namespace
 
 
+
+
+
 #endif	// _CA_ESRI_ASCIIGRID_HPP_