First Mosaic Dem implementation

Dockerfile

@@ -1,6 +1,6 @@
 FROM ghcr.io/sen2vm/sen2vm-build-env:latest AS launcher
 
-ENV SEN2VM_VERSION=1.2.0
+ENV SEN2VM_VERSION=1.3.0.rc1
 
 WORKDIR /Sen2vm
 

documentation/Input/DEM_CDSE_Download.md

@@ -23,7 +23,7 @@
 
 Digital Elevation Models (DEM) are essential for geolocation and orthorectification in Sen2VM. Users can use different types of DEM (cf [Inputs Description](../Input/input_description.md)).
 
-Sen2VM requires DEM data **organized per square degree** (see §[DEM format requirements](../Input/input_description.md#131-dem)). The recommended way to obtain Copernicus DEM in the correct format is to use the **[CDSE-Copernicus-DEM-downloader](https://github.com/senbox-org/CDSE-Copernicus-DEM-downloader)** tool, which downloads individual 1°×1° geocells from the Copernicus Data Space Ecosystem.
+Sen2VM requires DEM data **organized per square degree** (see §[DEM format requirements](../Input/input_description.md#131-dem)) but can now also handle mosaic of square degrees. The recommended way to obtain Copernicus DEM in the correct format is to use the **[CDSE-Copernicus-DEM-downloader](https://github.com/senbox-org/CDSE-Copernicus-DEM-downloader)** tool, which downloads individual 1°×1° geocells from the Copernicus Data Space Ecosystem.
 
 ---
 

documentation/Input/input_description.md

@@ -174,8 +174,11 @@ For this, as Sen2VM uses SXGEO (OREKIT/RUGGED), a GEOID and a DEM shall be used.
 
 Access to the DEM is provided via a path to a folder containing the dataset. 
 The DEM must meet the following requirements:
- * it should be split into files or folders (dynamically read) per square degrees,
- * each DEM file (per square degree) shall be readable by gdal.  
+
+ * it can be split into files or folders (dynamically read),
+ * each DEM file shall be readable by gdal.  
+
+Sen2VM also handles mosaic of square degree DEM.
 
  Examples of DEM structures can be found in [/src/test/resources/DEM/](/src/test/resources/DEM)
 

pom.xml

@@ -7,7 +7,7 @@
     <groupId>esa.sen2vm</groupId>
     <artifactId>sen2vm-core</artifactId>
     <!-- NOTE: remind to change also the SEN2VM_VERSION in the Dokerfile-->
-    <version>1.2.0</version>
+    <version>1.3.0.rc1</version>
 
     <name>sen2vm-core</name>
     <url>https://github.com/sen2vm/sen2vm-core</url>

src/main/java/esa/sen2vm/Sen2VM.java

@@ -22,6 +22,9 @@
 import java.io.IOException;
 import java.io.InputStream;
 import java.util.ArrayList;
+
+import java.util.Arrays;
+
 import java.util.Vector;
 import java.util.HashMap;
 import java.util.List;
@@ -44,7 +47,7 @@
 import esa.sen2vm.enums.DetectorInfo;
 import esa.sen2vm.exception.Sen2VMException;
 import esa.sen2vm.input.Configuration;
-import esa.sen2vm.input.GenericDemFileManager;
+import esa.sen2vm.input.DEM.GenericDemFileManager;
 import esa.sen2vm.input.OptionManager;
 import esa.sen2vm.input.Params;
 import esa.sen2vm.input.datastrip.DataStripManager;
@@ -181,7 +184,6 @@ public static void main( String[] args ) throws Sen2VMException, Exception
 
             //Using Sen2VM FileManager
             GenericDemFileManager demFileManager = new GenericDemFileManager(config.getDem());
-            demFileManager.buildMap(config.getDem());
 
             GeoidManager geoidManager = new GeoidManager(config.getGeoid(), isOverlappingTiles);
             DemManager demManager = new DemManager(
@@ -260,11 +262,11 @@ public static void main( String[] args ) throws Sen2VMException, Exception
             // Test if no grids exists already
             if (config.getOperation().equals(Sen2VMConstants.DIRECT))
             {
-                safeManager.testifDirectGridsToComputeAlreadyExist(detectors, bands) ;
+                safeManager.testifDirectGridsToComputeAlreadyExist(detectors, bands);
             }
             else
             {
-                safeManager.testifInverseGridsToComputeAlreadyExist(detectors, bands, config.getInverseLocOutputFolder()) ;
+                safeManager.testifInverseGridsToComputeAlreadyExist(detectors, bands, config.getInverseLocOutputFolder());
             }
 
             for (BandInfo bandInfo: bands)
@@ -339,7 +341,6 @@ public static void main( String[] args ) throws Sen2VMException, Exception
                         // Correction post build VRT
                         outputFileManager.correctGeoGrid(inputTIFs);
                         outputFileManager.correctVRT(vrtFileName);
-
                     }
 
                     // Inverse Loc case
@@ -352,6 +353,7 @@ public static void main( String[] args ) throws Sen2VMException, Exception
                         double[][] groundGrid = invGrid.get2DgridLatLon();
 
                         double[][] inverseLocGrid = simpleLocEngine.computeInverseLoc(sensorList.get(bandInfo.getNameWithB() + "/" + detectorInfo.getNameWithD()),  groundGrid, "EPSG:4326");
+
                         double[][][] grid3D = invGrid.get3Dgrid(inverseLocGrid, georefConventionOffsetPixel, -georefConventionOffsetLine);
 
                         String invFileName = datastrip.getCorrespondingInverseLocGrid(detectorInfo, bandInfo, config.getInverseLocOutputFolder());
@@ -372,11 +374,11 @@ public static void main( String[] args ) throws Sen2VMException, Exception
             }
             outputFileManager.writeInfoJson(config, bands, detectors, outputConfigPath);
         }
-        catch ( IOException exception )
+        catch (IOException exception)
         {
             throw new Sen2VMException(exception);
         }
-        catch ( SXGeoException exception )
+        catch (SXGeoException exception)
         {
             String newMessage = "";
             if(exception.toString().contains("Cant find bundle for base name S2GeoMessages"))

src/main/java/esa/sen2vm/input/DEM/DemTile.java

@@ -0,0 +1,31 @@
+package esa.sen2vm.input.DEM;
+
+public class DemTile
+{
+    public double minX;
+    public double maxX;
+    public double minY;
+    public double maxY;
+    public String filePath;
+
+    public DemTile(double a_minX, double a_maxX, double a_minY, double a_maxY, String filePathString)
+    {
+        minX = a_minX;
+        maxX = a_maxX;
+        minY = a_minY;
+        maxY = a_maxY;
+        filePath = filePathString;
+    }
+
+    public boolean containPoint(double x, double y)
+    {
+        return x >= minX && x <= maxX &&
+            y >= minY && y <= maxY;
+    }
+
+    @Override
+    public String toString()
+    {
+        return "DemTile (" + minX + ";" + maxX + "),(" + minY + ";" + maxY + ")" + filePath;
+    }
+}

src/main/java/esa/sen2vm/input/GenericDemFileManager.java → src/main/java/esa/sen2vm/input/DEM/GenericDemFileManager.java

@@ -14,7 +14,7 @@
 * See the License for the specific language governing permissions and
 * limitations under the License.*/
 
-package esa.sen2vm.input;
+package esa.sen2vm.input.DEM;
 
 import org.hipparchus.util.FastMath;
 import java.io.File;
@@ -24,6 +24,8 @@
 import java.nio.file.Files;
 import java.nio.file.Path;
 import java.util.Map;
+import java.util.ArrayList;
+import java.util.List;
 import java.util.HashMap;
 import java.util.logging.Logger;
 
@@ -32,9 +34,11 @@
 import org.gdal.gdalconst.gdalconstConstants;
 
 import org.sxgeo.input.dem.SrtmFileManager;
+
 import org.sxgeo.exception.SXGeoException;
 
 import esa.sen2vm.exception.Sen2VMException;
+import esa.sen2vm.input.DEM.DemTile;
 
 //Extend SrtmFileManager and not DemManager, as an check on the type of instanciation is made if isInstance of SrtmFileManager
 
@@ -50,22 +54,48 @@ public class GenericDemFileManager extends SrtmFileManager
 
     // Map dem filepath with a string that represents longitude/latitude
     // Example: with a SRTM tile on Madeira island located at longitude -16 and latitude 30
-    // the correponding map entry will be ("-16/30"="/DEMDIR/DEM_SRTM/w016/n30.dt1")
+    // the corresponding map entry will be ("-16/30"="/DEMDIR/DEM_SRTM/w016/n30.dt1")
     // Key corresponding to "longitude/latitude" and value corresponding to the dem filepath.
-    Map<String, String> demFilePathMap = new HashMap<>();
+    Map<Long, List<DemTile>> demGridMap = new HashMap<>();
 
     /**
      * {@inheritDoc}
      */
-    public GenericDemFileManager(String demRootDir)
+    public GenericDemFileManager(String demRootDir) throws Sen2VMException
     {
         super(demRootDir);
+        LOGGER.info("Loading DEM from  : " + demRootDir);
+        buildMap(demRootDir);
+    }
+
+
+    private long mapKey(int x, int y)
+    {
+        return (((long) x) << 32) | (y & 0xffffffffL);
+    }
+
+    private void addDemTile(DemTile d)
+    {
+        int xMin = (int)FastMath.floor(d.minX);
+        int xMax = (int)FastMath.floor(d.maxX);
+        int yMin = (int)FastMath.floor(d.minY);
+        int yMax = (int)FastMath.floor(d.maxY);
+
+        for (int x = xMin; x < xMax; x++)
+        {
+            for (int y = yMin; y < yMax; y++)
+            {
+                LOGGER.finer("Loading DEM from  : " + x + " " + y);
+                long key = mapKey(x, y);
+                demGridMap.computeIfAbsent(key, k -> new ArrayList<>()).add(d);
+            }
+        }
     }
 
     /**
      * Build a map that contains dem files
      */
-    public void buildMap(String directory) throws Sen2VMException
+    private void buildMap(String directory) throws Sen2VMException
     {
         try
         {
@@ -81,10 +111,14 @@ public void buildMap(String directory) throws Sen2VMException
                 else
                 {
                     String filePath = currentFile.getAbsolutePath();
-                    String lonlat = getLonLatFromFile(filePath);
-                    if (lonlat != null)
+                    LOGGER.finer("Loading DEM Tile  : " + filePath);
+                    LOGGER.info("Loading DEM Tile  : " + filePath);
+                    DemTile newDemTile = getDemTileFromFile(filePath);
+                    if (newDemTile != null)
                     {
-                        demFilePathMap.put(lonlat, filePath);
+                        LOGGER.finer("DEM Tile loaded : " + newDemTile.toString());
+                        LOGGER.info("DEM Tile loaded : " + newDemTile.toString());
+                        addDemTile(newDemTile);
                     }
                 }
             }
@@ -104,34 +138,7 @@ public boolean findRasterFile(String directory) throws SXGeoException
     {
         try
         {
-            Path dir = FileSystems.getDefault().getPath(directory);
-            DirectoryStream<Path> stream = Files.newDirectoryStream(dir);
-            boolean found = false;
-            for (Path path : stream)
-            {
-                if (!found)
-                {
-                    File currentFile = path.toFile();
-                    if (currentFile.isDirectory())
-                    {
-                        found = findRasterFile(currentFile.getAbsolutePath());
-                    }
-                    else
-                    {
-                        String filePath = currentFile.getAbsolutePath();
-                        if ( ( filePath.matches(".*.dt1") ) || ( filePath.matches(".*.dt2") )) {
-                            found = true;
-                        }
-                    }
-                    if (found)
-                    {
-                        stream.close();
-                        return true;
-                    }
-                }
-            }
-            stream.close();
-            throw new SXGeoException("NO_RASTER_FILE_FOUND_IN_DEM");
+            return (demGridMap.size() != 0);
         }
         catch (Exception e)
         {
@@ -145,8 +152,10 @@ public boolean findRasterFile(String directory) throws SXGeoException
     @Override
     protected String getRasterFilePath(double latitude, double longitude)
     {
-        double latFloor = FastMath.floor(FastMath.toDegrees(latitude));
-        double lonFloor = FastMath.floor(FastMath.toDegrees(longitude));
+        int latFloor = (int)FastMath.floor(FastMath.toDegrees(latitude));
+        int lonFloor = (int)FastMath.floor(FastMath.toDegrees(longitude));
+
+        LOGGER.finer("Searching DEM Tile for lat " + FastMath.toDegrees(latitude) + " lon " + FastMath.toDegrees(longitude));
 
         // when close to the anti-meridian
         if (lonFloor >= 180)
@@ -158,26 +167,38 @@ else if (lonFloor < -180)
             lonFloor += 360;
         }
 
-        String lonlat = (int) lonFloor + "/" + (int) latFloor;
-        String filePath = this.demFilePathMap.get(lonlat);
-        if (filePath == null)
+        long key = mapKey(lonFloor, latFloor);
+        List<DemTile> candidates = demGridMap.get(key);
+
+        if (candidates != null)
         {
-            filePath = "";
+            for (DemTile d: candidates)
+            {
+                if(d.containPoint(FastMath.toDegrees(longitude), FastMath.toDegrees(latitude)))
+                {
+                    // LOGGER.finer("DEM TILE FOUND FOR THIS LAT/LON");
+                    // LOGGER.finer(d.toString());
+                    LOGGER.info("DEM TILE FOUND FOR THIS LAT/LON");
+                    LOGGER.info(d.toString());
+                    return d.filePath;
+                }
+            }
         }
-        return filePath;
+
+        return "";
     }
 
     /**
-     * Get footprint information from file
+     * Get footprint information and create a DemTile object from file
      */
-    public String getLonLatFromFile(String filePath)
+    public DemTile getDemTileFromFile(String filePath)
     {
         gdal.AllRegister();
 
         Dataset dataset = gdal.Open(filePath, gdalconstConstants.GA_ReadOnly);
         if (dataset == null)
         {
-            LOGGER.severe("Error when reading  : " + gdal.GetLastErrorMsg());
+            LOGGER.severe("Error when reading : " + gdal.GetLastErrorMsg());
             //System.err.println("Error when reading  : " + gdal.GetLastErrorMsg());
             return null;
         }
@@ -186,17 +207,23 @@ public String getLonLatFromFile(String filePath)
 
         double minX = geoTransform[0];
         double maxY = geoTransform[3];
-        double pixelWidth = geoTransform[1];
-        double pixelHeight = geoTransform[5];
-
+        final double pixelWidth = geoTransform[1];
+        final double pixelHeight = geoTransform[5];
+        // DEM convention int Lat/Lon at center
+        minX += pixelWidth/2;
+        maxY += pixelHeight/2;
         int imageWidth = dataset.getRasterXSize();
         int imageHeight = dataset.getRasterYSize();
-
-        double maxX = minX + imageWidth * pixelWidth;
-        double minY = maxY + imageHeight * pixelHeight;
+        double maxX = minX + (imageWidth - 1) * pixelWidth;
+        double minY = maxY + (imageHeight - 1) * pixelHeight;
         dataset.delete();
 
-        String lonlat = Math.round(minX) + "/" + Math.round(minY);
-        return lonlat;
+        minX = Math.round(minX);
+        maxX = Math.round(maxX);
+        minY = Math.round(minY);
+        maxY = Math.round(maxY);
+
+        return new DemTile(minX, maxX, minY, maxY, filePath);
     }
+
 }

src/test/java/esa/sen2vm/Sen2VMDirectTest.java

@@ -349,9 +349,9 @@ public void testDirectParallelisation()
     @Test
     public void testDirectDem()
     {
-        String[] detectors = new String[]{"01"};
+        String[] detectors = new String[]{"08"};
         String[] bands = new String[]{"B01"};
-        String[] testsDem = new String[]{"dem_1", "dem_2", "dem_3", "dem_4"};
+        String[] testsDem = new String[]{"dem_1", "dem_2", "dem_3", "dem_4", "dem_5", "dem_6"};
         int stepBand10m = 600; // corresponding to 6 kms
 
         try