From 526bb53b703a3a187825dbcc2994e6fca74f2a4d Mon Sep 17 00:00:00 2001 From: afischerdev Date: Sat, 2 Mar 2024 10:44:38 +0100 Subject: [PATCH] added create elev image --- .../CreateElevationRasterImage.java | 266 ++++++++++++++ .../java/btools/mapcreator/HgtReader.java | 342 ++++++++++++++++++ 2 files changed, 608 insertions(+) create mode 100644 brouter-map-creator/src/main/java/btools/mapcreator/CreateElevationRasterImage.java create mode 100644 brouter-map-creator/src/main/java/btools/mapcreator/HgtReader.java diff --git a/brouter-map-creator/src/main/java/btools/mapcreator/CreateElevationRasterImage.java b/brouter-map-creator/src/main/java/btools/mapcreator/CreateElevationRasterImage.java new file mode 100644 index 0000000..50b8d81 --- /dev/null +++ b/brouter-map-creator/src/main/java/btools/mapcreator/CreateElevationRasterImage.java @@ -0,0 +1,266 @@ +package btools.mapcreator; + +import java.awt.Color; +import java.awt.image.BufferedImage; +import java.awt.image.DataBufferInt; +import java.io.BufferedInputStream; +import java.io.File; +import java.io.FileInputStream; +import java.io.FileOutputStream; +import java.io.InputStream; +import java.util.HashMap; +import java.util.Map; +import java.util.Map.Entry; +import java.util.TreeMap; + +import javax.imageio.ImageIO; + + +public class CreateElevationRasterImage { + + final static boolean DEBUG = false; + + int[] data; + ElevationRaster lastSrtmRaster; + Map srtmmap; + int lastSrtmLonIdx; + int lastSrtmLatIdx; + String srtmdir; + boolean missingData; + Map colorMap; + + private void createImage(double lon, double lat, String dir, String imageName, int maxX, int maxY, int downscale, String format) throws Exception { + srtmdir = dir; + if (format.equals("hgt")) { + createImageFromHgt(lon, lat, dir, imageName, maxX, maxY); + return; + } + srtmmap = new HashMap<>(); + lastSrtmLonIdx = -1; + lastSrtmLatIdx = -1; + lastSrtmRaster = null; + NodeData n = new NodeData(1, lon, lat); + ElevationRaster srtm = srtmForNode(n.ilon, n.ilat); + System.out.println("srtm " + srtm.toString()); + //System.out.println("srtm elev " + srtm.getElevation(n.ilon, n.ilat)); + double[] pos = getElevationPos(srtm, n.ilon, n.ilat); + System.out.println("srtm pos " + pos[0] + " " + pos[1]); // Arrays.toString(pos)); + short[] raster = srtm.eval_array; + int rasterX = srtm.ncols; + int rasterY = srtm.nrows; + + int tileSize = 1000 / downscale; + int sizeX = (maxX); + int sizeY = (maxY); + int[] imgraster = new int[sizeX * sizeY]; + System.out.println("srtm target " + sizeX + " " + sizeY + " (" + rasterX + " " + rasterY + ")"); + for (int y = 0; y < sizeY; y++) { + for (int x = 0; x < sizeX; x++) { + //short e = getElevationXY(srtm, /*pos[0] + */(sizeY - y) * downscale, /*pos[1] +*/ (x * downscale)); + short e = get(srtm, sizeY - y, x); + + if (e == Short.MIN_VALUE) { + imgraster[sizeY * y + x] = 0xffff; + } else { + //imgraster[sizeY * y + x] = getColorForHeight((short)(e/4)); //(int)(e/4.); + imgraster[sizeY * y + x] = getColorForHeight(e); + } + } + } + if (DEBUG) { + String out = "short "; + for (int i = 0; i < 100; i++) { + out += " " + get(srtm, sizeY - 0, i); + } + System.out.println(out); + } + + BufferedImage argbImage = new BufferedImage(sizeX, sizeY, BufferedImage.TYPE_INT_ARGB); + data = ((DataBufferInt) argbImage.getRaster().getDataBuffer()).getData(); + + for (int y = 0; y < sizeY; y++) { + for (int x = 0; x < sizeX; x++) { + int v0 = imgraster[sizeX * y + x]; + + int rgb; + if (v0 != 0xffff) + rgb = 0xff000000 | v0; //(v0 << 8); + else + rgb = 0xff000000; + data[y * sizeX + x] = rgb; + } + } + + ImageIO.write(argbImage, "png", new FileOutputStream(imageName)); + } + + private void createImageFromHgt(double lon, double lat, String dir, String imageName, int maxX, int maxY) throws Exception { + HgtReader rdr = new HgtReader(dir); + short[] data = rdr.getElevationDataFromHgt(lat, lon); + int size = (data != null ? data.length : 0); + int rowlen = (int) Math.sqrt(size); + System.out.println("hgt size " + rowlen); + int sizeX = (maxX); + int sizeY = (maxY); + int[] imgraster = new int[sizeX * sizeY]; + + for (int y = 0; y < sizeY; y++) { + for (int x = 0; x < sizeX; x++) { + short e = data[(rowlen * y) + x]; + if (e == HgtReader.HGT_VOID) { + imgraster[sizeY * y + x] = 0xffff; + } else if (e == 0) { + imgraster[sizeY * y + x] = 0xffff; + } else { + imgraster[sizeY * y + x] = getColorForHeight((short) (e)); //(int)(e/4.); + } + } + } + if (DEBUG) { + String out = "short "; + for (int i = 0; i < 100; i++) { + out += " " + data[i]; + } + System.out.println(out); + } + BufferedImage argbImage = new BufferedImage(sizeX, sizeY, BufferedImage.TYPE_INT_ARGB); + int[] idata = ((DataBufferInt) argbImage.getRaster().getDataBuffer()).getData(); + + for (int y = 0; y < sizeY; y++) { + for (int x = 0; x < sizeX; x++) { + int v0 = imgraster[sizeX * y + x]; + + int rgb; + if (v0 != 0xffff) + rgb = 0xff000000 | v0; //(v0 << 8); + else + rgb = 0xff000000; + idata[y * sizeX + x] = rgb; + } + } + + ImageIO.write(argbImage, "png", new FileOutputStream(imageName)); + + } + + public double[] getElevationPos(ElevationRaster srtm, int ilon, int ilat) { + double lon = ilon / 1000000. - 180.; + double lat = ilat / 1000000. - 90.; + + double dcol = (lon - srtm.xllcorner) / srtm.cellsize - 0.5; + double drow = (lat - srtm.yllcorner) / srtm.cellsize - 0.5; + int row = (int) drow; + int col = (int) dcol; + if (col < 0) col = 0; + if (row < 0) row = 0; + + return new double[]{drow, dcol}; + } + + private short get(ElevationRaster srtm, int r, int c) { + short e = srtm.eval_array[(srtm.nrows - 1 - r) * srtm.ncols + c]; + if (e == Short.MIN_VALUE) missingData = true; + return e; + } + + public short getElevationXY(ElevationRaster srtm, double drow, double dcol) { + int row = (int) drow; + int col = (int) dcol; + if (col < 0) col = 0; + if (col >= srtm.ncols - 1) col = srtm.ncols - 2; + if (row < 0) row = 0; + if (row >= srtm.nrows - 1) row = srtm.nrows - 2; + double wrow = drow - row; + double wcol = dcol - col; + missingData = false; + + double eval = (1. - wrow) * (1. - wcol) * get(srtm, row, col) + + (wrow) * (1. - wcol) * get(srtm, row + 1, col) + + (1. - wrow) * (wcol) * get(srtm, row, col + 1) + + (wrow) * (wcol) * get(srtm, row + 1, col + 1); + + return missingData ? Short.MIN_VALUE : (short) (eval * 4); + } + + + int getColorForHeight(short h) { + if (colorMap == null) { + colorMap = new TreeMap<>(); + colorMap.put((short) 0, new Color(102, 153, 153)); + colorMap.put((short) 1, new Color(0, 102, 0)); + colorMap.put((short) 500, new Color(251, 255, 128)); + colorMap.put((short) 1200, new Color(224, 108, 31)); + colorMap.put((short) 2500, new Color(200, 55, 55)); + colorMap.put((short) 4000, new Color(215, 244, 244)); + colorMap.put((short) 8000, new Color(255, 244, 244)); + } + Color lastColor = null; + short lastKey = 0; + for (Entry entry : colorMap.entrySet()) { + short key = entry.getKey(); + Color value = entry.getValue(); + if (key == h) return value.getRGB(); + if (lastColor != null && lastKey < h && key > h) { + double between = (double) (h - lastKey) / (key - lastKey); + return mixColors(value, lastColor, between); + } + lastColor = value; + lastKey = key; + } + return 0; + } + + public int mixColors(Color color1, Color color2, double percent) { + double inverse_percent = 1.0 - percent; + int redPart = (int) (color1.getRed() * percent + color2.getRed() * inverse_percent); + int greenPart = (int) (color1.getGreen() * percent + color2.getGreen() * inverse_percent); + int bluePart = (int) (color1.getBlue() * percent + color2.getBlue() * inverse_percent); + return new Color(redPart, greenPart, bluePart).getRGB(); + } + + private ElevationRaster srtmForNode(int ilon, int ilat) throws Exception { + int srtmLonIdx = (ilon + 5000000) / 5000000; + int srtmLatIdx = (654999999 - ilat) / 5000000 - 100; // ugly negative rounding... + + if (srtmLonIdx == lastSrtmLonIdx && srtmLatIdx == lastSrtmLatIdx) { + return lastSrtmRaster; + } + lastSrtmLonIdx = srtmLonIdx; + lastSrtmLatIdx = srtmLatIdx; + + String slonidx = "0" + srtmLonIdx; + String slatidx = "0" + srtmLatIdx; + String filename = "srtm_" + slonidx.substring(slonidx.length() - 2) + "_" + slatidx.substring(slatidx.length() - 2); + + lastSrtmRaster = srtmmap.get(filename); + if (lastSrtmRaster == null && !srtmmap.containsKey(filename)) { + File f = new File(new File(srtmdir), filename + ".bef"); + if (f.exists()) { + System.out.println("*** reading: " + f); + try { + InputStream isc = new BufferedInputStream(new FileInputStream(f)); + lastSrtmRaster = new ElevationRasterCoder().decodeRaster(isc); + isc.close(); + } catch (Exception e) { + System.out.println("**** ERROR reading " + f + " ****"); + } + srtmmap.put(filename, lastSrtmRaster); + return lastSrtmRaster; + } + + srtmmap.put(filename, lastSrtmRaster); + } + return lastSrtmRaster; + } + + public static void main(String[] args) throws Exception { + if (args.length < 6) { + System.out.println("usage: java CreateLidarImage [type]"); + System.out.println("\nwhere: type = [bef|hgt] downscale = [1|2|4|..]"); + return; + } + String format = args.length == 8 ? args[7] : "bef"; + new CreateElevationRasterImage().createImage(Double.parseDouble(args[0]), Double.parseDouble(args[1]), args[2], args[3], + Integer.parseInt(args[4]), Integer.parseInt(args[5]), Integer.parseInt(args[6]), format); + } +} 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..f161d7d --- /dev/null +++ b/brouter-map-creator/src/main/java/btools/mapcreator/HgtReader.java @@ -0,0 +1,342 @@ +// 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("HGT 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(); + } + } + } + System.out.println("*** reading: " + f.getName() + " " + cache.get(file)); + } + } + + // read elevation value + return readElevation(lat, lon); + } catch (FileNotFoundException e) { + System.err.println("HGT Get elevation " + 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; + } + } + + public static short[] getElevationDataFromHgt(double lat, double lon) { + try { + if (lon < 0) lon += 1; + if (lat < 0) lat += 1; + String file = getHgtFileName(lat, lon); + if (DEBUG) System.out.println("HGT buffer " + file + " for " + lat + " " + lon); + + ShortBuffer data = 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... + data = readHgtFile(fullPath); + } 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()); + data = readHgtStream(zis); + break; + } + zis.closeEntry(); + } + } finally { + zis.close(); + } + } + } + + + System.out.println("*** reading: " + f.getName() + " " + (data != null ? data.limit() : -1)); + if (data != null) { + short[] array = new short[data.limit()]; + data.get(array); + return array; + } + return null; + } catch (FileNotFoundException e) { + System.err.println("HGT Get elevation " + lat + ", " + lon + " failed: => " + e.getMessage()); + // no problem... file not there + return null; + } catch (Exception ioe) { + // oops... + ioe.printStackTrace(System.err); + // fallback + return null; + } + } + + @SuppressWarnings("resource") + private static ShortBuffer readHgtFile(String file) throws Exception { + if (file == null) throw new Exception("no hgt 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 hgt 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("HGT 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 HGT 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(); + } + } + + public static void main(String[] args) throws Exception { + System.out.println("*** HGT position values and enhance elevation"); + if (args.length == 3) { + HgtReader elevReader = new HgtReader(args[0]); + double lon = Double.parseDouble(args[1]); + double lat = Double.parseDouble(args[2]); + // check hgt direct + double elev = elevReader.getElevationFromHgt(lat, lon); + System.out.println("-----> elv for hgt " + lat + ", " + lon + " = " + elev); + + } + } +}