diff --git a/brouter-map-creator/src/main/java/btools/mapcreator/ConvertLidarTile.java b/brouter-map-creator/src/main/java/btools/mapcreator/ConvertLidarTile.java index 65f18aa..65307ce 100644 --- a/brouter-map-creator/src/main/java/btools/mapcreator/ConvertLidarTile.java +++ b/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; + + } + } diff --git a/brouter-map-creator/src/main/java/btools/mapcreator/HgtReader.java b/brouter-map-creator/src/main/java/btools/mapcreator/HgtReader.java new file mode 100644 index 0000000..b232bf3 --- /dev/null +++ b/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 + *
+ * Class HgtReader reads data from SRTM HGT files. Currently this class is restricted to a resolution of 3 arc seconds. + *
+ * SRTM data files are available at the NASA SRTM site
+ *
+ * @author Oliver Wieland <oliver.wieland@online.de>
+ */
+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
+ * See also stackexchange.com
+ *
+ * @param lat, lon the coordinate to get the elevation data for
+ * @return the elevation value or Double.NaN
, 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 [N|S]nn[W|E]mmm.hgt where nn is the integral latitude
+ * without decimals and mmm 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();
+ }
+ }
+}
diff --git a/brouter-map-creator/src/main/java/btools/mapcreator/PosUnifier.java b/brouter-map-creator/src/main/java/btools/mapcreator/PosUnifier.java
index 2943052..6fade10 100644
--- a/brouter-map-creator/src/main/java/btools/mapcreator/PosUnifier.java
+++ b/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