change hgt to ConvertLidarTile

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 {
@@ -26,7 +30,7 @@ public class ConvertLidarTile {
       for (; ; ) {
         ZipEntry ze = zis.getNextEntry();
         if (ze.getName().endsWith(".hgt")) {
-          readHgtFromStream(zis, rowOffset, colOffset);
+          readHgtFromStream(zis, rowOffset, colOffset, 1200);
           return;
         }
       }
@@ -35,13 +39,40 @@ 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!
+        int i0 = dis.read();
+
+        if (i0 == -1 || i1 == -1)
+          throw new RuntimeException("unexcepted end of file reading hgt entry!");
+
+        short val = (short) ((i1 << 8) | i0);
+
+        if (val == NODATA2) {
+          val = NODATA;
+        }
+
+        setPixel(row, col, val);
+      }
+    }
+  }
+
+  private static void readHgtFromFile(File f, int rowOffset, int colOffset)
+    throws Exception {
+    DataInputStream dis = new DataInputStream(new BufferedInputStream(new FileInputStream(f)));
+
+    for (int ir = 0; ir < ROW_LENGTH + 1; ir++) {
+      int row = rowOffset + ir;
+
+      for (int ic = 0; ic < ROW_LENGTH + 1; ic++) {
         int col = colOffset + ic;
 
         int i1 = dis.read(); // msb first!
@@ -186,4 +217,91 @@ public class ConvertLidarTile {
     doConvert(args[1], ilon_base, ilat_base, filename30);
   }
 
+  public SrtmRaster getRasterDirect(File f, double lon, double lat) throws Exception {
+    if (f.length() > ((SRTM3_ROW_LENGTH + 1) * (SRTM3_ROW_LENGTH + 1) * 2)) {
+      ROW_LENGTH = SRTM1_ROW_LENGTH;
+    } else {
+      ROW_LENGTH = SRTM3_ROW_LENGTH;
+    }
+    System.out.println("read file " + f + " rl " + 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;
+      }
+    }
+
+    readHgtFromFile(f, 0, 0);
+
+    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;
+  }
+
+  public SrtmRaster getRasterZip(File f, double lon, double lat) throws Exception {
+
+    ZipInputStream zis = new ZipInputStream(new BufferedInputStream(new FileInputStream(f)));
+    try {
+      for (; ; ) {
+        ZipEntry ze = zis.getNextEntry();
+        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;
+          }
+          System.out.println("read file " + f + " rl " + 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();
+    }
+
+
+    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/PosUnifier.java

@@ -32,7 +32,6 @@ public class PosUnifier extends MapCreatorBase {
   private int lastSrtmLatIdx;
   private SrtmRaster lastSrtmRaster;
   private String srtmdir;
-  private HgtReader hgtReader;
 
   private CompactLongSet borderNids;
 
@@ -46,13 +45,12 @@ public class PosUnifier extends MapCreatorBase {
       double lat = Double.parseDouble(args[2]);
 
       NodeData n = new NodeData(1, lon, lat);
-      SrtmRaster srtm = posu.srtmForNode(n.ilon, n.ilat);
+      SrtmRaster srtm = posu.hgtForNode(n.ilon, n.ilat);
       short selev = Short.MIN_VALUE;
       if (srtm == null) {
-        selev = posu.hgtForNode(n.ilon, n.ilat);
+        srtm = posu.srtmForNode(n.ilon, n.ilat);
-      } else {
-        selev = srtm.getElevation(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;
@@ -98,14 +96,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;
+    SrtmRaster srtm = hgtForNode(n.ilon, n.ilat);
     if (srtm == null) {
-      n.selev = hgtForNode(n.ilon, n.ilat);
+      srtm = srtmForNode(n.ilon, n.ilat);
-    } else {
-      n.selev = srtm.getElevation(n.ilon, n.ilat);
     }
 
+    if (srtm != null) n.selev = srtm.getElevation(n.ilon, n.ilat);
     findUniquePos(n);
 
     n.writeTo(nodesOutStream);
@@ -182,7 +179,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 {
@@ -201,6 +197,8 @@ public class PosUnifier extends MapCreatorBase {
       if (f.exists()) {
         try {
           lastSrtmRaster = new SrtmData(f).getRaster();
+          srtmmap.put(filename, lastSrtmRaster);
+          return lastSrtmRaster;
         } catch (Exception e) {
           System.out.println("**** ERROR reading " + f + " ****");
         }
@@ -210,21 +208,45 @@ public class PosUnifier extends MapCreatorBase {
     return lastSrtmRaster;
   }
 
-  private short hgtForNode(int ilon, int ilat) throws Exception {
+  private SrtmRaster hgtForNode(int ilon, int ilat) throws Exception {
     double lon = (ilon - 180000000) / 1000000.;
     double lat = (ilat - 90000000) / 1000000.;
 
-    try {
+    String filename = buildHgtFilename(lat, lon);
-      if (hgtReader == null) hgtReader = new HgtReader(srtmdir);
+    lastSrtmRaster = srtmmap.get(filename);
-      double res = hgtReader.getElevationFromHgt(lat, lon);
+    if (lastSrtmRaster == null) {
-      // System.out.println("**** reading " + res + " ****");
+      File f = new File(new File(srtmdir), filename + ".hgt");
-      if (Double.isNaN(res)) return Short.MIN_VALUE;
+      if (f.exists()) {
-      return (short) (res*4);
+        lastSrtmRaster = new ConvertLidarTile().getRasterDirect(f, lon, lat);
-    } catch (Exception e) {
+        srtmmap.put(filename, lastSrtmRaster);
-      System.out.println("**** ERROR reading hgt " + lon + " " + lat + " ****");
+        return lastSrtmRaster;
+      }
+      f = new File(new File(srtmdir), filename + ".zip");
+      if (f.exists()) {
+        lastSrtmRaster = new ConvertLidarTile().getRasterZip(f, lon, lat);
+        srtmmap.put(filename, lastSrtmRaster);
+        return lastSrtmRaster;
+      }
+    }
+    return lastSrtmRaster;
   }
 
-    return Short.MIN_VALUE;
+  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("%s%02d%s%03d", latPref, lat, lonPref, lon);
   }
 
   private void resetSrtm() {
@@ -232,7 +254,6 @@ public class PosUnifier extends MapCreatorBase {
     lastSrtmLonIdx = -1;
     lastSrtmLatIdx = -1;
     lastSrtmRaster = null;
-    if (hgtReader != null) hgtReader.clear();
   }
 
 }