Rebuilt dem (#15)

* RebuildDEM option now working. Additional fix of geotif handling for some coordinate systems. --------- Co-authored-by: Jake Langham <[email protected]> Co-authored-by: Jake Langham <[email protected]>
jakelangham · Dec 9, 2023 · e8569a8 · e8569a8
1 parent 2bde668
commit e8569a8
Show file tree

Hide file tree

Showing 3 changed files with 89 additions and 12 deletions.
diff --git a/src/Output.f90 b/src/Output.f90
@@ -38,7 +38,7 @@
 module output_module
 
    use set_precision_module, only: wp
-   use utilities_module, only: CheckFileExists, PathTrail, Int2String, KahanAdd, pair
+   use utilities_module, only: AddToOrderedVector, CheckFileExists, PathTrail, Int2String, KahanAdd, pair
    use grid_module, only: GridCoords, GridType, TileList, TileType
    use closures_module, only : FlowSquaredSpeedSlopeAligned, GeometricCorrectionFactor
    use runsettings_module, only : RunSet

diff --git a/src/RasterData.cpp b/src/RasterData.cpp
@@ -55,17 +55,22 @@ RasterData::RasterData(const char *filename) {
         // Get EPSG projection reference
 	    OGRSpatialReference sref = GetSpatialReference(poDataset);
         latlon = sref.EPSGTreatsAsLatLong(); // Check whether CRS uses lat-long coordinates
-
+
+        // Try to find the CRS Authority Code
+        // this approach uses the GDAL FindMatches function which could return several possible
+        // codes -- the first one is the 'best' match
+        int matches;
+        int* conf;
+        OGRSpatialReferenceH *srefs = sref.FindMatches(NULL, &matches, &conf);
+        const char* a_code = OSRGetAuthorityCode(srefs[0],NULL);
+        std::string auth_code = a_code;
+
 	    // Retrieve EPSG code and store
         try {
-	        if (strcmp(sref.GetAuthorityName(NULL), "EPSG")==0) {
-	    	    EPSG_code = atoi(sref.GetAuthorityCode(NULL));
-	        }
-            else {
-                throw (poDataset->GetFileList());
-            }
+	    	    EPSG_code = stoi(auth_code);
         }
-        catch (char **fileList) {
+        catch (std::exception const & e) {
+            char **fileList = poDataset->GetFileList();
             std::cerr << "Error: no EPSG authority in " << fileList << std::endl;
             std::cerr << "Check files and ensure they are properly georeferenced" << std::endl;
             exit(EXIT_FAILURE);

diff --git a/src/dem.f90 b/src/dem.f90
@@ -30,7 +30,7 @@ module dem_module
    use set_precision_module, only: wp, pi
    use runsettings_module, only: RunSet, FortranRasterData
    use grid_module, only: GridType
-   use utilities_module, only: pair
+   use utilities_module, only: pair, CheckFileExists
    use messages_module, only: InfoMessage, WarningMessage, FatalErrorMessage
    use varstring_module, only: varString
    use interpolate_2d_module, only: Bicubic
@@ -45,6 +45,56 @@ module dem_module
 
 contains
 
+   function CheckDEMvrtExists(RunParams) result(exists)
+      implicit none
+
+      type(RunSet), intent(in) :: RunParams
+      logical :: exists
+
+      type(varString) :: dem_file
+
+      dem_file = RunParams%out_path + 'DEM.vrt'
+      exists = CheckFileExists(dem_file)
+
+      return
+
+   end function CheckDEMvrtExists
+
+   function CheckDEMProperties(RunParams, deltaX, deltaY, dem_path, dem_file) result(good_dem)
+      implicit none
+
+      type(RunSet), intent(in) :: RunParams
+      real(kind=wp), intent(in) :: deltaX
+      real(kind=wp), intent(in) :: deltaY
+      type(varString), intent(in) :: dem_path
+      type(varString), intent(in) :: dem_file
+      logical :: good_dem
+
+      type(pair) :: res
+      real(kind=c_double) :: minE, maxE, minN, maxN
+
+      type(FortranRasterData) :: demInfo
+
+      good_dem = .TRUE.
+
+      call DomainExtent(RunParams, deltaX, deltaY, minE, maxE, minN, maxN, res, pad=.True.)
+
+      call RasterInfo(dem_path, dem_file, demInfo)
+
+      if (demInfo%ReadSuccess .ne. 1) call FatalErrorMessage("Cannot read DEM file "//dem_path%s//dem_file%s)
+
+      if (demInfo%pixelWidth .ne. deltaX) good_dem = .FALSE.
+      if (-demInfo%pixelHeight .ne. deltaY) good_dem = .FALSE.
+
+      if (demInfo%NE%first .ne. maxN) good_dem = .FALSE.
+      if (demInfo%NE%second .ne. maxE) good_dem = .FALSE.
+      if (demInfo%SW%first .ne. minN) good_dem = .FALSE.
+      if (demInfo%SW%second .ne. minE) good_dem = .FALSE.
+
+      return
+
+   end function CheckDEMProperties
+
    ! If the user has specified a topography that requires an external datafile,
    ! such as a DEM, we try to load it in preparation for running a simulation.
    ! Note that we need to have read the input file (to populate RunParams) and
@@ -57,7 +107,7 @@ subroutine LoadDEM(RunParams, grid)
 
       select case (RunParams%Topog%s)
          case ('srtm', 'SRTM', 'raster', 'Raster')
-            call InfoMessage("Preparing topographic data")
+            call InfoMessage("Preparing topographic data.")
             call PrepareDEM(RunParams, grid%deltaX, grid%deltaY)
       end select
 
@@ -85,7 +135,8 @@ subroutine PrepareDEM(RunParams, deltaX, deltaY, interp_flag_in)
 
       type(varString) :: geotif
       character(len=1,kind=c_char) :: cgeotif(len_trim(trim(RunParams%demPath%s))+1+len_trim(RunParams%RasterFile%s) + 1)
-      logical FileExists
+      logical :: FileExists
+      logical :: dem_exists, good_dem
 
       if (present(interp_flag_in)) interp_flag = interp_flag_in
 
@@ -94,6 +145,27 @@ subroutine PrepareDEM(RunParams, deltaX, deltaY, interp_flag_in)
 
       call DomainExtent(RunParams, deltaX, deltaY, minE, maxE, minN, maxN, res, pad=.True.)
 
+      dem_exists = CheckDEMvrtExists(RunParams)
+
+      if (.not. RunParams%RebuildDEM) then
+         ! Requested to use existing DEM.vrt file.  Check this exists
+         if (dem_exists) then
+            ! DEM.vrt exists, so check it covers the domain and has the correct resolution
+            good_dem = CheckDEMProperties(RunParams, deltaX, deltaY, RunParams%out_path, varString("DEM.vrt"))
+            if (good_dem) then
+               ! DEM.vrt is good -- report and return
+               call InfoMessage("Using existing DEM.vrt file.")
+               return
+            else
+               ! DEM.vrt is not good -- warn and move on to rebuild
+               call WarningMessage("DEM.vrt is not compatible with domain extent and/or resolution. Rebuilding DEM.")
+            end if
+         else
+            ! DEM.vrt does not exist, so warn and move on to rebuild.
+            call WarningMessage("DEM.vrt not found in "//RunParams%out_path%s// ". Building DEM.")
+         end if
+      end if
+
       select case (RunParams%Topog%s)
        case ('srtm')
          raster = .FALSE.