Merge pull request #548 from afischerdev/elev-source

Add a hgt Reader

brouter-map-creator/src/main/java/btools/mapcreator/ConvertLidarTile.java

@@ -14,10 +14,14 @@ import java.util.zip.ZipInputStream;
 public class ConvertLidarTile {
   public static int NROWS;
   public static int NCOLS;
+  public static int ROW_LENGTH;
 
   public static final short NODATA2 = -32767; // hgt-formats nodata
   public static final short NODATA = Short.MIN_VALUE;
 
+  public static final int SRTM3_ROW_LENGTH = 1200; // number of elevation values per line
+  public static final int SRTM1_ROW_LENGTH = 3600; //-- New file resolution is 3601x3601
+
   static short[] imagePixels;
 
   private static void readHgtZip(String filename, int rowOffset, int colOffset) throws Exception {
@@ -25,8 +29,9 @@ public class ConvertLidarTile {
     try {
       for (; ; ) {
         ZipEntry ze = zis.getNextEntry();
+        if (ze == null) break;
         if (ze.getName().endsWith(".hgt")) {
-          readHgtFromStream(zis, rowOffset, colOffset);
+          readHgtFromStream(zis, rowOffset, colOffset, 1200);
           return;
         }
       }
@@ -35,13 +40,13 @@ public class ConvertLidarTile {
     }
   }
 
-  private static void readHgtFromStream(InputStream is, int rowOffset, int colOffset)
+  private static void readHgtFromStream(InputStream is, int rowOffset, int colOffset, int row_length)
     throws Exception {
     DataInputStream dis = new DataInputStream(new BufferedInputStream(is));
-    for (int ir = 0; ir < 1201; ir++) {
+    for (int ir = 0; ir < row_length + 1; ir++) {
       int row = rowOffset + ir;
 
-      for (int ic = 0; ic < 1201; ic++) {
+      for (int ic = 0; ic < row_length + 1; ic++) {
         int col = colOffset + ic;
 
         int i1 = dis.read(); // msb first!
@@ -61,7 +66,6 @@ public class ConvertLidarTile {
     }
   }
 
-
   private static void setPixel(int row, int col, short val) {
     if (row >= 0 && row < NROWS && col >= 0 && col < NCOLS) {
       imagePixels[row * NCOLS + col] = val;
@@ -186,4 +190,81 @@ public class ConvertLidarTile {
     doConvert(args[1], ilon_base, ilat_base, filename30);
   }
 
+  public SrtmRaster getRaster(File f, double lon, double lat) throws Exception {
+
+    if (f.getName().toLowerCase().endsWith(".zip")) {
+      ZipInputStream zis = new ZipInputStream(new BufferedInputStream(new FileInputStream(f)));
+      try {
+        for (; ; ) {
+          ZipEntry ze = zis.getNextEntry();
+          if (ze == null) break;
+          if (ze.getName().toLowerCase().endsWith(".hgt")) {
+            if (ze.getSize() > ((SRTM3_ROW_LENGTH + 1) * (SRTM3_ROW_LENGTH + 1) * 2)) {
+              ROW_LENGTH = SRTM1_ROW_LENGTH;
+            } else {
+              ROW_LENGTH = SRTM3_ROW_LENGTH;
+            }
+            // stay a 1x1 raster
+            NROWS = ROW_LENGTH + 1;
+            NCOLS = ROW_LENGTH + 1;
+
+            imagePixels = new short[NROWS * NCOLS]; // 650 MB !
+
+            // prefill as NODATA
+            for (int row = 0; row < NROWS; row++) {
+              for (int col = 0; col < NCOLS; col++) {
+                imagePixels[row * NCOLS + col] = NODATA;
+              }
+            }
+            readHgtFromStream(zis, 0, 0, ROW_LENGTH);
+            break;
+          }
+        }
+      } finally {
+        zis.close();
+      }
+    } else {
+      if (f.length() > ((SRTM3_ROW_LENGTH + 1) * (SRTM3_ROW_LENGTH + 1) * 2)) {
+        ROW_LENGTH = SRTM1_ROW_LENGTH;
+      } else {
+        ROW_LENGTH = SRTM3_ROW_LENGTH;
+      }
+      // stay a 1x1 raster
+      NROWS = ROW_LENGTH + 1;
+      NCOLS = ROW_LENGTH + 1;
+
+      imagePixels = new short[NROWS * NCOLS]; // 650 MB !
+
+      // prefill as NODATA
+      for (int row = 0; row < NROWS; row++) {
+        for (int col = 0; col < NCOLS; col++) {
+          imagePixels[row * NCOLS + col] = NODATA;
+        }
+      }
+      FileInputStream dis = new FileInputStream(f);
+      try {
+        readHgtFromStream(dis, 0, 0, ROW_LENGTH);
+      } finally {
+        dis.close();
+      }
+    }
+
+    System.out.println("read file " + f + " rl " + ROW_LENGTH);
+    
+    boolean halfCol5 = false; // no halfcol tiles in lidar data (?)
+
+    SrtmRaster raster = new SrtmRaster();
+    raster.nrows = NROWS;
+    raster.ncols = NCOLS;
+    raster.halfcol = halfCol5;
+    raster.noDataValue = NODATA;
+    raster.cellsize = 1. / (double) ROW_LENGTH;
+    raster.xllcorner = (int) (lon < 0 ? lon - 1 : lon); //onDegreeStart - raster.cellsize;
+    raster.yllcorner = (int) (lat < 0 ? lat - 1 : lat); //latDegreeStart - raster.cellsize;
+    raster.eval_array = imagePixels;
+
+    return raster;
+
+  }
+
 }

brouter-map-creator/src/main/java/btools/mapcreator/HgtReader.java

@@ -0,0 +1,269 @@
+// License: GPL. For details, see LICENSE file.
+package btools.mapcreator;
+
+import java.io.BufferedInputStream;
+import java.io.File;
+import java.io.FileInputStream;
+import java.io.FileNotFoundException;
+import java.io.InputStream;
+import java.nio.ByteBuffer;
+import java.nio.ByteOrder;
+import java.nio.ShortBuffer;
+import java.nio.channels.FileChannel;
+import java.util.HashMap;
+import java.util.Map;
+import java.util.zip.ZipEntry;
+import java.util.zip.ZipInputStream;
+
+/**
+ * adapted from https://github.com/JOSM/josm-plugins/blob/master/ElevationProfile/src/org/openstreetmap/josm/plugins/elevation/HgtReader.java
+ * <p>
+ * Class HgtReader reads data from SRTM HGT files. Currently this class is restricted to a resolution of 3 arc seconds.
+ * <p>
+ * SRTM data files are available at the <a href="http://dds.cr.usgs.gov/srtm/version2_1/SRTM3">NASA SRTM site</a>
+ *
+ * @author Oliver Wieland &lt;oliver.wieland@online.de&gt;
+ */
+public class HgtReader {
+  final static boolean DEBUG = false;
+
+  private static final int SECONDS_PER_MINUTE = 60;
+
+  public static final String HGT_EXT = ".hgt";
+  public static final String ZIP_EXT = ".zip";
+
+  // alter these values for different SRTM resolutions
+  public static final int HGT3_RES = 3; // resolution in arc seconds
+  public static final int HGT3_ROW_LENGTH = 1201; // number of elevation values per line
+  public static final int HGT_VOID = -32768; // magic number which indicates 'void data' in HGT file
+  public static final int HGT1_RES = 1;  // <<- The new SRTM is 1-ARCSEC
+  public static final int HGT1_ROW_LENGTH = 3601; //-- New file resolution is 3601x3601
+  /**
+   * The 'no elevation' data magic.
+   */
+  public static double NO_ELEVATION = Double.NaN;
+
+  private static String srtmFolder = "";
+
+  private static final Map<String, ShortBuffer> cache = new HashMap<>();
+
+  public HgtReader(String folder) {
+    srtmFolder = folder;
+  }
+
+  public static double getElevationFromHgt(double lat, double lon) {
+    try {
+      String file = getHgtFileName(lat, lon);
+      if (DEBUG) System.out.println("SRTM buffer " + file + " for " + lat + " " + lon);
+
+      // given area in cache?
+      if (!cache.containsKey(file)) {
+
+        // fill initial cache value. If no file is found, then
+        // we use it as a marker to indicate 'file has been searched
+        // but is not there'
+        cache.put(file, null);
+        // Try all resource directories
+        //for (String location : Main.pref.getAllPossiblePreferenceDirs())
+        {
+          String fullPath = new File(srtmFolder, file + HGT_EXT).getPath();
+          File f = new File(fullPath);
+          if (f.exists()) {
+            // found something: read HGT file...
+            ShortBuffer data = readHgtFile(fullPath);
+            // ... and store result in cache
+            cache.put(file, data);
+            //break;
+          } else {
+            fullPath = new File(srtmFolder, file + ZIP_EXT).getPath();
+            f = new File(fullPath);
+            if (f.exists()) {
+              ZipInputStream zis = new ZipInputStream(new BufferedInputStream(new FileInputStream(f)));
+              try {
+                for (; ; ) {
+                  ZipEntry ze = zis.getNextEntry();
+                  if (ze == null) break;
+                  if (ze.getName().toLowerCase().endsWith(HGT_EXT)) {
+                    // System.out.println("read zip " + ze.getName());
+                    ShortBuffer data = readHgtStream(zis);
+                    // ... and store result in cache
+                    cache.put(file, data);
+                    break;
+                  }
+                  zis.closeEntry();
+                }
+              } finally {
+                zis.close();
+              }
+            }
+          }
+        }
+      }
+
+      // read elevation value
+      return readElevation(lat, lon);
+    } catch (FileNotFoundException e) {
+      System.err.println("SRTM Get elevation from HGT " + lat + ", " + lon + " failed: => " + e.getMessage());
+      // no problem... file not there
+      return NO_ELEVATION;
+    } catch (Exception ioe) {
+      // oops...
+      ioe.printStackTrace(System.err);
+      // fallback
+      return NO_ELEVATION;
+    }
+  }
+
+  @SuppressWarnings("resource")
+  private static ShortBuffer readHgtFile(String file) throws Exception {
+    if (file == null) throw new Exception("no srmt file " + file);
+
+    FileChannel fc = null;
+    ShortBuffer sb = null;
+    try {
+      // Eclipse complains here about resource leak on 'fc' - even with 'finally' clause???
+      fc = new FileInputStream(file).getChannel();
+      // choose the right endianness
+
+      ByteBuffer bb = ByteBuffer.allocateDirect((int) fc.size());
+      while (bb.remaining() > 0) fc.read(bb);
+
+      bb.flip();
+      //sb = bb.order(ByteOrder.LITTLE_ENDIAN).asShortBuffer();
+      sb = bb.order(ByteOrder.BIG_ENDIAN).asShortBuffer();
+    } finally {
+      if (fc != null) fc.close();
+    }
+
+    return sb;
+  }
+
+  // @SuppressWarnings("resource")
+  private static ShortBuffer readHgtStream(InputStream zis) throws Exception {
+    if (zis == null) throw new Exception("no srmt stream ");
+
+    ShortBuffer sb = null;
+    try {
+      // choose the right endianness
+
+      byte[] bytes = zis.readAllBytes();
+      ByteBuffer bb = ByteBuffer.allocate(bytes.length);
+      bb.put(bytes, 0, bytes.length);
+      //while (bb.remaining() > 0) zis.read(bb, 0, size);
+
+      //ByteBuffer bb = ByteBuffer.allocate(zis.available());
+      //Channels.newChannel(zis).read(bb);
+      bb.flip();
+      //sb = bb.order(ByteOrder.LITTLE_ENDIAN).asShortBuffer();
+      sb = bb.order(ByteOrder.BIG_ENDIAN).asShortBuffer();
+    } finally {
+
+    }
+
+    return sb;
+  }
+
+  /**
+   * Reads the elevation value for the given coordinate.
+   * <p>
+   * See also <a href="http://gis.stackexchange.com/questions/43743/how-to-extract-elevation-from-hgt-file">stackexchange.com</a>
+   *
+   * @param lat, lon the coordinate to get the elevation data for
+   * @return the elevation value or <code>Double.NaN</code>, if no value is present
+   */
+  public static double readElevation(double lat, double lon) {
+    String tag = getHgtFileName(lat, lon);
+
+    ShortBuffer sb = cache.get(tag);
+
+    if (sb == null) {
+      return NO_ELEVATION;
+    }
+
+    if (DEBUG) System.out.println("SRTM buffer size " + sb.capacity() + " limit " + sb.limit());
+    try {
+      int rowLength = HGT3_ROW_LENGTH;
+      int resolution = HGT3_RES;
+      if (sb.capacity() > (HGT3_ROW_LENGTH * HGT3_ROW_LENGTH)) {
+        rowLength = HGT1_ROW_LENGTH;
+        resolution = HGT1_RES;
+      }
+      // see http://gis.stackexchange.com/questions/43743/how-to-extract-elevation-from-hgt-file
+      double fLat = frac(lat) * SECONDS_PER_MINUTE;
+      double fLon = frac(lon) * SECONDS_PER_MINUTE;
+
+      // compute offset within HGT file
+      int row = (int) Math.round((fLat) * SECONDS_PER_MINUTE / resolution);
+      int col = (int) Math.round((fLon) * SECONDS_PER_MINUTE / resolution);
+      if (lon < 0) col = rowLength - col - 1;
+      if (lat > 0) row = rowLength - row - 1;
+
+
+      //row = rowLength - row;
+      int cell = (rowLength * (row)) + col;
+      //int cell = ((rowLength * (latitude)) + longitude);
+
+      if (DEBUG)
+        System.out.println("Read SRTM elevation data from row/col/cell " + row + "," + col + ", " + cell + ", " + sb.limit());
+
+      // valid position in buffer?
+      if (cell < sb.limit()) {
+        short ele = sb.get(cell);
+        // check for data voids
+        if (ele == HGT_VOID) {
+          return NO_ELEVATION;
+        } else {
+          return ele;
+        }
+      } else {
+        return NO_ELEVATION;
+      }
+    } catch (Exception e) {
+      System.err.println("error at " + lon + " " + lat + " ");
+      e.printStackTrace();
+    }
+    return NO_ELEVATION;
+  }
+
+  /**
+   * Gets the associated HGT file name for the given way point. Usually the
+   * format is <tt>[N|S]nn[W|E]mmm.hgt</tt> where <i>nn</i> is the integral latitude
+   * without decimals and <i>mmm</i> is the longitude.
+   *
+   * @param llat,llon the coordinate to get the filename for
+   * @return the file name of the HGT file
+   */
+  public static String getHgtFileName(double llat, double llon) {
+    int lat = (int) llat;
+    int lon = (int) llon;
+
+    String latPref = "N";
+    if (lat < 0) {
+      latPref = "S";
+      lat = -lat + 1;
+    }
+    String lonPref = "E";
+    if (lon < 0) {
+      lonPref = "W";
+      lon = -lon + 1;
+    }
+
+    return String.format("%s%02d%s%03d", latPref, lat, lonPref, lon);
+  }
+
+  public static double frac(double d) {
+    long iPart;
+    double fPart;
+
+    // Get user input
+    iPart = (long) d;
+    fPart = d - iPart;
+    return Math.abs(fPart);
+  }
+
+  public static void clear() {
+    if (cache != null) {
+      cache.clear();
+    }
+  }
+}

brouter-map-creator/src/main/java/btools/mapcreator/PosUnifier.java

@@ -7,6 +7,7 @@ import java.io.File;
 import java.io.FileInputStream;
 import java.io.InputStream;
 import java.util.HashMap;
+import java.util.Locale;
 import java.util.Map;
 
 import btools.util.CompactLongSet;
@@ -37,7 +38,24 @@ public class PosUnifier extends MapCreatorBase {
 
   public static void main(String[] args) throws Exception {
     System.out.println("*** PosUnifier: Unify position values and enhance elevation");
-    if (args.length != 5) {
+    if (args.length == 3) {
+      PosUnifier posu = new PosUnifier();
+      posu.srtmdir = (args[0]);
+      posu.srtmmap = new HashMap<>();
+      double lon = Double.parseDouble(args[1]);
+      double lat = Double.parseDouble(args[2]);
+
+      NodeData n = new NodeData(1, lon, lat);
+      SrtmRaster srtm = posu.hgtForNode(n.ilon, n.ilat);
+      short selev = Short.MIN_VALUE;
+      if (srtm == null) {
+        srtm = posu.srtmForNode(n.ilon, n.ilat);
+      }
+      if (srtm != null) selev = srtm.getElevation(n.ilon, n.ilat);
+      posu.resetSrtm();
+      System.out.println("-----> selv for " + lat + ", " + lon + " = " + selev + " = " + (selev / 4.));
+      return;
+    } else if (args.length != 5) {
       System.out.println("usage: java PosUnifier <node-tiles-in> <node-tiles-out> <bordernids-in> <bordernodes-out> <srtm-data-dir>");
       return;
     }
@@ -79,9 +97,13 @@ public class PosUnifier extends MapCreatorBase {
 
   @Override
   public void nextNode(NodeData n) throws Exception {
-    SrtmRaster srtm = srtmForNode(n.ilon, n.ilat);
+    n.selev = Short.MIN_VALUE;
-    n.selev = srtm == null ? Short.MIN_VALUE : srtm.getElevation(n.ilon, n.ilat);
+    SrtmRaster srtm = hgtForNode(n.ilon, n.ilat);
+    if (srtm == null) {
+      srtm = srtmForNode(n.ilon, n.ilat);
+    }
 
+    if (srtm != null) n.selev = srtm.getElevation(n.ilon, n.ilat);
     findUniquePos(n);
 
     n.writeTo(nodesOutStream);
@@ -93,6 +115,7 @@ public class PosUnifier extends MapCreatorBase {
   @Override
   public void nodeFileEnd(File nodeFile) throws Exception {
     nodesOutStream.close();
+    resetSrtm();
   }
 
   private boolean checkAdd(int lon, int lat) {
@@ -157,7 +180,6 @@ public class PosUnifier extends MapCreatorBase {
     lastSrtmRaster = srtmmap.get(filename);
     if (lastSrtmRaster == null && !srtmmap.containsKey(filename)) {
       File f = new File(new File(srtmdir), filename + ".bef");
-      System.out.println("checking: " + f + " ilon=" + ilon + " ilat=" + ilat);
       if (f.exists()) {
         System.out.println("*** reading: " + f);
         try {
@@ -172,10 +194,12 @@ public class PosUnifier extends MapCreatorBase {
       }
 
       f = new File(new File(srtmdir), filename + ".zip");
-      System.out.println("reading: " + f + " ilon=" + ilon + " ilat=" + ilat);
+      // System.out.println("reading: " + f + " ilon=" + ilon + " ilat=" + ilat);
       if (f.exists()) {
         try {
           lastSrtmRaster = new SrtmData(f).getRaster();
+          srtmmap.put(filename, lastSrtmRaster);
+          return lastSrtmRaster;
         } catch (Exception e) {
           System.out.println("**** ERROR reading " + f + " ****");
         }
@@ -185,8 +209,49 @@ public class PosUnifier extends MapCreatorBase {
     return lastSrtmRaster;
   }
 
+  private SrtmRaster hgtForNode(int ilon, int ilat) throws Exception {
+    double lon = (ilon - 180000000) / 1000000.;
+    double lat = (ilat - 90000000) / 1000000.;
+
+    String filename = buildHgtFilename(lat, lon);
+    lastSrtmRaster = srtmmap.get(filename);
+    if (lastSrtmRaster == null) {
+      File f = new File(new File(srtmdir), filename + ".hgt");
+      if (f.exists()) {
+        lastSrtmRaster = new ConvertLidarTile().getRaster(f, lon, lat);
+        srtmmap.put(filename, lastSrtmRaster);
+        return lastSrtmRaster;
+      }
+      f = new File(new File(srtmdir), filename + ".zip");
+      if (f.exists()) {
+        lastSrtmRaster = new ConvertLidarTile().getRaster(f, lon, lat);
+        srtmmap.put(filename, lastSrtmRaster);
+        return lastSrtmRaster;
+      }
+    }
+    return lastSrtmRaster;
+  }
+
+  private String buildHgtFilename(double llat, double llon) {
+    int lat = (int) llat;
+    int lon = (int) llon;
+
+    String latPref = "N";
+    if (lat < 0) {
+      latPref = "S";
+      lat = -lat + 1;
+    }
+    String lonPref = "E";
+    if (lon < 0) {
+      lonPref = "W";
+      lon = -lon + 1;
+    }
+
+    return String.format(Locale.US, "%s%02d%s%03d", latPref, lat, lonPref, lon);
+  }
+
   private void resetSrtm() {
-    srtmmap = new HashMap<String, SrtmRaster>();
+    srtmmap = new HashMap<>();
     lastSrtmLonIdx = -1;
     lastSrtmLatIdx = -1;
     lastSrtmRaster = null;