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 65307ce..46f9c57 100644 --- a/brouter-map-creator/src/main/java/btools/mapcreator/ConvertLidarTile.java +++ b/brouter-map-creator/src/main/java/btools/mapcreator/ConvertLidarTile.java @@ -5,22 +5,26 @@ import java.io.BufferedOutputStream; import java.io.DataInputStream; import java.io.File; import java.io.FileInputStream; +import java.io.FileNotFoundException; import java.io.FileOutputStream; import java.io.InputStream; import java.io.OutputStream; +import java.util.Arrays; import java.util.zip.ZipEntry; import java.util.zip.ZipInputStream; public class ConvertLidarTile { - public static int NROWS; - public static int NCOLS; - public static int ROW_LENGTH; + private static int NROWS; + private static int NCOLS; 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 + private static final String HGT_FILE_EXT = ".hgt"; + private static final int HGT_BORDER_OVERLAP = 1; + private static final int HGT_3ASEC_ROWS = 1201; // 3 arc second resolution (90m) + private static final int HGT_3ASEC_FILE_SIZE = HGT_3ASEC_ROWS * HGT_3ASEC_ROWS * Short.BYTES; + private static final int HGT_1ASEC_ROWS = 3601; // 1 arc second resolution (30m) static short[] imagePixels; @@ -30,8 +34,8 @@ public class ConvertLidarTile { for (; ; ) { ZipEntry ze = zis.getNextEntry(); if (ze == null) break; - if (ze.getName().endsWith(".hgt")) { - readHgtFromStream(zis, rowOffset, colOffset, 1200); + if (ze.getName().toLowerCase().endsWith(HGT_FILE_EXT)) { + readHgtFromStream(zis, rowOffset, colOffset, HGT_3ASEC_ROWS); return; } } @@ -40,20 +44,20 @@ public class ConvertLidarTile { } } - private static void readHgtFromStream(InputStream is, int rowOffset, int colOffset, int row_length) + private static void readHgtFromStream(InputStream is, int rowOffset, int colOffset, int rowLength) throws Exception { DataInputStream dis = new DataInputStream(new BufferedInputStream(is)); - for (int ir = 0; ir < row_length + 1; ir++) { + for (int ir = 0; ir < rowLength; ir++) { int row = rowOffset + ir; - for (int ic = 0; ic < row_length + 1; ic++) { + for (int ic = 0; ic < rowLength; 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!"); + throw new RuntimeException("unexpected end of file reading hgt entry!"); short val = (short) ((i1 << 8) | i0); @@ -191,80 +195,56 @@ public class ConvertLidarTile { } public SrtmRaster getRaster(File f, double lon, double lat) throws Exception { + long fileSize; + InputStream inputStream; 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; - } + for (; ; ) { + ZipEntry ze = zis.getNextEntry(); + if (ze == null) { + throw new FileNotFoundException(f.getName() + " doesn't contain a " + HGT_FILE_EXT + " file."); + } + if (ze.getName().toLowerCase().endsWith(HGT_FILE_EXT)) { + fileSize = ze.getSize(); + inputStream = zis; + 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(); - } + fileSize = f.length(); + inputStream = new FileInputStream(f); } - System.out.println("read file " + f + " rl " + ROW_LENGTH); - - boolean halfCol5 = false; // no halfcol tiles in lidar data (?) + int rowLength; + if (fileSize > HGT_3ASEC_FILE_SIZE) { + rowLength = HGT_1ASEC_ROWS; + } else { + rowLength = HGT_3ASEC_ROWS; + } + + // stay at 1 deg * 1 deg raster + NROWS = rowLength; + NCOLS = rowLength; + + imagePixels = new short[NROWS * NCOLS]; + + // prefill as NODATA + Arrays.fill(imagePixels, NODATA); + readHgtFromStream(inputStream, 0, 0, rowLength); + inputStream.close(); SrtmRaster raster = new SrtmRaster(); raster.nrows = NROWS; raster.ncols = NCOLS; - raster.halfcol = halfCol5; + raster.halfcol = false; // assume full resolution raster.noDataValue = NODATA; - raster.cellsize = 1. / (double) ROW_LENGTH; + raster.cellsize = 1. / (double) (rowLength - HGT_BORDER_OVERLAP); 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 deleted file mode 100644 index b232bf3..0000000 --- a/brouter-map-creator/src/main/java/btools/mapcreator/HgtReader.java +++ /dev/null @@ -1,269 +0,0 @@ -// 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 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. - *

- * 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(); - } - } -}