From 254aff19b8370a973ca79771185ec58283a1a9c3 Mon Sep 17 00:00:00 2001 From: afischerdev Date: Fri, 27 Oct 2023 16:21:58 +0200 Subject: [PATCH 01/15] added 1sec srtm use --- .../java/btools/mapcreator/PosUnifier.java | 84 ++++++++++++++----- 1 file changed, 64 insertions(+), 20 deletions(-) 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 bccbf7001..10e660554 100644 --- a/brouter-map-creator/src/main/java/btools/mapcreator/PosUnifier.java +++ b/brouter-map-creator/src/main/java/btools/mapcreator/PosUnifier.java @@ -23,6 +23,9 @@ * @author ab */ public class PosUnifier extends MapCreatorBase { + + public static final boolean UseLidarRd5FileName = false; + private DiffCoderDataOutputStream nodesOutStream; private DiffCoderDataOutputStream borderNodesOut; private File nodeTilesOut; @@ -33,6 +36,7 @@ public class PosUnifier extends MapCreatorBase { private int lastSrtmLatIdx; private SrtmRaster lastSrtmRaster; private String srtmdir; + private String srtmfallbackdir; private CompactLongSet borderNids; @@ -46,25 +50,40 @@ public static void main(String[] args) throws Exception { 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; + SrtmRaster srtm = null; + /* + // check hgt direct + srtm = posu.hgtForNode(n.ilon, n.ilat); + if (srtm != null) { + selev = srtm.getElevation(n.ilon, n.ilat); + } else { + System.out.println("hgtForNode no data"); + } + posu.resetSrtm(); + System.out.println("-----> selv for hgt " + lat + ", " + lon + " = " + selev + " = " + (selev / 4.)); + srtm = null; + 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.)); + System.out.println("-----> selv for bef " + lat + ", " + lon + " = " + selev + " = " + (selev / 4.)); return; - } else if (args.length != 5) { - System.out.println("usage: java PosUnifier "); + } else if (args.length != 5 && args.length != 6) { + System.out.println("usage: java PosUnifier [srtm-fallback-data-dir]"); + System.out.println("or java PosUnifier "); return; } - new PosUnifier().process(new File(args[0]), new File(args[1]), new File(args[2]), new File(args[3]), args[4]); + new PosUnifier().process(new File(args[0]), new File(args[1]), new File(args[2]), new File(args[3]), args[4], (args.length == 6 ? args[5] : null)); } - public void process(File nodeTilesIn, File nodeTilesOut, File bordernidsinfile, File bordernodesoutfile, String srtmdir) throws Exception { + public void process(File nodeTilesIn, File nodeTilesOut, File bordernidsinfile, File bordernodesoutfile, String srtmdir, String srtmfallbackdir) throws Exception { this.nodeTilesOut = nodeTilesOut; this.srtmdir = srtmdir; + this.srtmfallbackdir = srtmfallbackdir; // read border nids set DataInputStream dis = createInStream(bordernidsinfile); @@ -178,15 +197,17 @@ private SrtmRaster srtmForNode(int ilon, int ilat) throws Exception { 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); + String filename; + if (UseLidarRd5FileName) { + filename = genFilenameRd5(ilon, ilat); + } else { + filename = genFilenameXY(srtmLonIdx, srtmLatIdx); + } 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 RasterCoder().decodeRaster(isc); @@ -194,19 +215,24 @@ private SrtmRaster srtmForNode(int ilon, int ilat) throws Exception { } catch (Exception e) { System.out.println("**** ERROR reading " + f + " ****"); } + System.out.println("*** reading: " + f + " " + lastSrtmRaster.ncols); srtmmap.put(filename, lastSrtmRaster); return lastSrtmRaster; } - - f = new File(new File(srtmdir), filename + ".zip"); - // System.out.println("reading: " + f + " ilon=" + ilon + " ilat=" + ilat); - if (f.exists()) { - try { - lastSrtmRaster = new SrtmData(f).getRaster(); + if (srtmfallbackdir != null) { + f = new File(new File(srtmfallbackdir), filename + ".bef"); + if (f.exists()) { + try { + InputStream isc = new BufferedInputStream(new FileInputStream(f)); + //lastSrtmRaster = new StatRasterCoder().decodeRaster(isc); + lastSrtmRaster = new RasterCoder().decodeRaster(isc); + isc.close(); + } catch (Exception e) { + System.out.println("**** ERROR reading " + f + " ****"); + } + System.out.println("*** reading: " + f + " " + lastSrtmRaster.cellsize); srtmmap.put(filename, lastSrtmRaster); return lastSrtmRaster; - } catch (Exception e) { - System.out.println("**** ERROR reading " + f + " ****"); } } srtmmap.put(filename, lastSrtmRaster); @@ -214,6 +240,24 @@ private SrtmRaster srtmForNode(int ilon, int ilat) throws Exception { return lastSrtmRaster; } + static String genFilenameXY(int srtmLonIdx, int srtmLatIdx) { + String slonidx = "0" + srtmLonIdx; + String slatidx = "0" + srtmLatIdx; + return "srtm_" + slonidx.substring(slonidx.length() - 2) + "_" + slatidx.substring(slatidx.length() - 2); + } + + static String genFilenameRd5(int ilon, int ilat) { + int lonDegree = ilon / 1000000; + int latDegree = ilat / 1000000; + int lonMod5 = lonDegree % 5; + int latMod5 = latDegree % 5; + lonDegree = lonDegree - 180 - lonMod5; + latDegree = latDegree - 90 - latMod5; + return String.format("srtm_%s_%s", lonDegree < 0 ? "W" + (-lonDegree) : "E" + lonDegree, + latDegree < 0 ? "S" + (-latDegree) : "N" + latDegree); + } + + private SrtmRaster hgtForNode(int ilon, int ilat) throws Exception { double lon = (ilon - 180000000) / 1000000.; double lat = (ilat - 90000000) / 1000000.; @@ -222,13 +266,13 @@ private SrtmRaster hgtForNode(int ilon, int ilat) throws Exception { // don't block lastSrtmRaster SrtmRaster srtm = srtmmap.get(filename); if (srtm == null) { - File f = new File(new File(srtmdir), filename + ".hgt"); + File f = new File(new File(srtmdir), filename + ".zip"); if (f.exists()) { srtm = new ConvertLidarTile().getRaster(f, lon, lat); srtmmap.put(filename, srtm); return srtm; } - f = new File(new File(srtmdir), filename + ".zip"); + f = new File(new File(srtmdir), filename + ".hgt"); if (f.exists()) { srtm = new ConvertLidarTile().getRaster(f, lon, lat); srtmmap.put(filename, srtm); From c7786f03ecb410d059daa683b1a04722896ab5e0 Mon Sep 17 00:00:00 2001 From: afischerdev Date: Fri, 27 Oct 2023 16:22:52 +0200 Subject: [PATCH 02/15] added 1sec bef generation --- .../btools/mapcreator/ConvertLidarTile.java | 276 +++++++++++++++--- 1 file changed, 236 insertions(+), 40 deletions(-) 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 2a00e96be..cab1d0b7f 100644 --- a/brouter-map-creator/src/main/java/btools/mapcreator/ConvertLidarTile.java +++ b/brouter-map-creator/src/main/java/btools/mapcreator/ConvertLidarTile.java @@ -9,13 +9,15 @@ import java.io.FileOutputStream; import java.io.InputStream; import java.io.OutputStream; +import java.util.ArrayList; import java.util.Arrays; +import java.util.List; import java.util.zip.ZipEntry; import java.util.zip.ZipInputStream; public class ConvertLidarTile { - private static int NROWS; - private static int NCOLS; + + public static final boolean DEBUG = false; public static final short NODATA2 = -32767; // hgt-formats nodata public static final short NODATA = Short.MIN_VALUE; @@ -25,17 +27,89 @@ public class ConvertLidarTile { 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) + private static final int SRTM3_ROW_LENGTH = 1200; // number of elevation values per line + private static final int SRTM1_ROW_LENGTH = 3600; + private static final boolean SRTM_NO_ZERO = true; + + private int NROWS; + private int NCOLS; + private int ROW_LENGTH; + private short[] imagePixels; + + + public static void main(String[] args) throws Exception { + if (args.length == 3 || args.length == 4 || args.length == 5) { + String filename90 = args[0]; + if ("all".equals(filename90)) { + //if (DEBUG) + System.out.println("lidar convert all "); + new ConvertLidarTile().doConvertAll(args[1], args[2], (args.length > 3 ? args[3] : null), (args.length == 5 ? args[4] : null)); + return; + } + // old filenames only + String filename30 = filename90 + ".bef"; //filename90.substring(0, filename90.length() - 3) + "bef"; + + int srtmLonIdx = Integer.parseInt(filename90.substring(5, 7).toLowerCase()); + int srtmLatIdx = Integer.parseInt(filename90.substring(8, 10).toLowerCase()); + + int ilon_base = (srtmLonIdx - 1) * 5 - 180; + int ilat_base = 150 - srtmLatIdx * 5 - 90; + int row_length = SRTM3_ROW_LENGTH; + String fallbackdir = null; + if (args.length > 3) { + row_length = (Integer.parseInt(args[3]) == 1 ? SRTM1_ROW_LENGTH : SRTM3_ROW_LENGTH); + fallbackdir = (args.length == 5 ? args[4] : null); + } + //if (DEBUG) + System.out.println("lidar convert " + ilon_base + " " + ilat_base + " from " + srtmLonIdx + " " + srtmLatIdx + " f: " + filename90 + " rowl " + row_length); + + new ConvertLidarTile().doConvert(args[1], ilon_base, ilat_base, args[2] + "/" + filename30, row_length, fallbackdir); + } else { + System.out.println("usage: java [arc seconds (1 or 3,default=3)] [hgt-fallback-data-dir]"); + System.out.println("or java all [arc seconds (1 or 3, default=3)] [hgt-fallback-data-dir]"); + return; + } + } + + private void doConvertAll(String hgtdata, String outdir, String rlen, String hgtfallbackdata) throws Exception { + int row_length = SRTM3_ROW_LENGTH; + if (rlen != null) { + row_length = (Integer.parseInt(rlen) == 1 ? SRTM1_ROW_LENGTH : SRTM3_ROW_LENGTH); + } + String filename30; + for (int ilon_base = -180; ilon_base < 180; ilon_base += 5) { + for (int ilat_base = 85; ilat_base > -90; ilat_base -= 5) { + if (PosUnifier.UseLidarRd5FileName) { + filename30 = genFilenameRd5(ilon_base, ilat_base); + } else { + filename30 = genFilenameOld(ilon_base, ilat_base); + } + if (DEBUG) + System.out.println("lidar convert all: " + filename30); + doConvert(hgtdata, ilon_base, ilat_base, outdir + "/" + filename30, row_length, hgtfallbackdata); + } + } + } - static short[] imagePixels; + static String genFilenameOld(int ilon_base, int ilat_base) { + int srtmLonIdx = ((ilon_base + 180) / 5) + 1; + int srtmLatIdx = (60 - ilat_base) / 5; + return String.format("srtm_%02d_%02d.bef", srtmLonIdx, srtmLatIdx); + } + + static String genFilenameRd5(int ilon_base, int ilat_base) { + return String.format("srtm_%s_%s.bef", ilon_base < 0 ? "W" + (-ilon_base) : "E" + ilon_base, + ilat_base < 0 ? "S" + (-ilat_base) : "N" + ilat_base); + } - private static void readHgtZip(String filename, int rowOffset, int colOffset) throws Exception { + private void readHgtZip(String filename, int rowOffset, int colOffset, int row_length) throws Exception { ZipInputStream zis = new ZipInputStream(new BufferedInputStream(new FileInputStream(filename))); try { for (; ; ) { ZipEntry ze = zis.getNextEntry(); if (ze == null) break; if (ze.getName().toLowerCase().endsWith(HGT_FILE_EXT)) { - readHgtFromStream(zis, rowOffset, colOffset, HGT_3ASEC_ROWS); + readHgtFromStream(zis, rowOffset, colOffset, row_length, 1); return; } } @@ -44,14 +118,14 @@ private static void readHgtZip(String filename, int rowOffset, int colOffset) th } } - private static void readHgtFromStream(InputStream is, int rowOffset, int colOffset, int rowLength) + private void readHgtFromStream(InputStream is, int rowOffset, int colOffset, int rowLength, int scale) throws Exception { DataInputStream dis = new DataInputStream(new BufferedInputStream(is)); for (int ir = 0; ir < rowLength; ir++) { - int row = rowOffset + ir; + int row = rowOffset + ir * scale; for (int ic = 0; ic < rowLength; ic++) { - int col = colOffset + ic; + int col = colOffset + ic * scale; int i1 = dis.read(); // msb first! int i0 = dis.read(); @@ -64,19 +138,52 @@ private static void readHgtFromStream(InputStream is, int rowOffset, int colOffs if (val == NODATA2) { val = NODATA; } - - setPixel(row, col, val); + if (scale == 3) { + setPixel(row, col, val); + setPixel(row + 1, col, val); + setPixel(row + 2, col, val); + setPixel(row, col + 1, val); + setPixel(row + 1, col + 1, val); + setPixel(row + 2, col + 1, val); + setPixel(row, col + 2, val); + setPixel(row + 1, col + 2, val); + setPixel(row + 2, col + 2, val); + } else { + setPixel(row, col, val); + } } } } - private static void setPixel(int row, int col, short val) { + private void readFallbackFile(File file, int rowOffset, int colOffset, int row_length) + throws Exception { + int rowLength; + int scale; + if (file.length() > HGT_3ASEC_FILE_SIZE) { + rowLength = HGT_1ASEC_ROWS; + scale = 1; + } else { + rowLength = HGT_3ASEC_ROWS; + scale = 3; + } + if (DEBUG) + System.out.println("read fallback: " + file + " " + rowLength); + + FileInputStream fis = new FileInputStream(file); + try { + readHgtFromStream(fis, rowOffset, colOffset, rowLength, scale); + } finally { + fis.close(); + } + } + + private void setPixel(int row, int col, short val) { if (row >= 0 && row < NROWS && col >= 0 && col < NCOLS) { imagePixels[row * NCOLS + col] = val; } } - private static short getPixel(int row, int col) { + private short getPixel(int row, int col) { if (row >= 0 && row < NROWS && col >= 0 && col < NCOLS) { return imagePixels[row * NCOLS + col]; } @@ -84,36 +191,137 @@ private static short getPixel(int row, int col) { } - public static void doConvert(String inputDir, int lonDegreeStart, int latDegreeStart, String outputFile) throws Exception { + public void doConvert(String inputDir, int lonDegreeStart, int latDegreeStart, String outputFile, int row_length, String hgtfallbackdata) throws Exception { int extraBorder = 0; - NROWS = 5 * 1200 + 1 + 2 * extraBorder; - NCOLS = 5 * 1200 + 1 + 2 * extraBorder; + List foundList = new ArrayList<>(); + List notfoundList = new ArrayList<>(); - imagePixels = new short[NROWS * NCOLS]; // 650 MB ! + boolean found = false; - // prefill as NODATA - for (int row = 0; row < NROWS; row++) { - for (int col = 0; col < NCOLS; col++) { - imagePixels[row * NCOLS + col] = NODATA; + if (row_length == SRTM1_ROW_LENGTH) { + // check for sources w/o border + for (int latIdx = 0; latIdx < 5; latIdx++) { + int latDegree = latDegreeStart + latIdx; + + for (int lonIdx = 0; lonIdx < 5; lonIdx++) { + int lonDegree = lonDegreeStart + lonIdx; + + String filename = inputDir + "/" + formatLat(latDegree) + formatLon(lonDegree) + ".zip"; + File f = new File(filename); + if (f.exists() && f.length() > 0) { + found = true; + break; + } + } } + } else { + // ignore when srtm3 + found = true; + } + if (found) { // init when found + NROWS = 5 * row_length + 1 + 2 * extraBorder; + NCOLS = 5 * row_length + 1 + 2 * extraBorder; + imagePixels = new short[NROWS * NCOLS]; // 650 MB ! + + // prefill as NODATA + Arrays.fill(imagePixels, NODATA); + } else { + if (DEBUG) + System.out.println("none 1sec data: " + lonDegreeStart + " " + latDegreeStart); + return; } for (int latIdx = -1; latIdx <= 5; latIdx++) { int latDegree = latDegreeStart + latIdx; - int rowOffset = extraBorder + (4 - latIdx) * 1200; + int rowOffset = extraBorder + (4 - latIdx) * row_length; for (int lonIdx = -1; lonIdx <= 5; lonIdx++) { int lonDegree = lonDegreeStart + lonIdx; - int colOffset = extraBorder + lonIdx * 1200; + int colOffset = extraBorder + lonIdx * row_length; String filename = inputDir + "/" + formatLat(latDegree) + formatLon(lonDegree) + ".zip"; File f = new File(filename); if (f.exists() && f.length() > 0) { - System.out.println("exist: " + filename); - readHgtZip(filename, rowOffset, colOffset); + if (DEBUG) + System.out.println("exist: " + filename); + readHgtZip(filename, rowOffset, colOffset, row_length + 1); } else { - System.out.println("none : " + filename); + if (hgtfallbackdata != null) { + String filenamehgt = hgtfallbackdata + "/" + formatLat(latDegree) + formatLon(lonDegree) + ".hgt"; + f = new File(filenamehgt); + if (f.exists() && f.length() > 0) { + readFallbackFile(f, rowOffset, colOffset, row_length + 1); + /*if (imagePixels == null) { + imagePixels = new short[NROWS * NCOLS]; + Arrays.fill(imagePixels, NODATA); + found = true; + } + */ + /* + int rowLength; + int arcspace; + if (f.length() > HGT_3ASEC_FILE_SIZE) { + rowLength = HGT_1ASEC_ROWS; + arcspace = 1; + } else { + rowLength = HGT_3ASEC_ROWS; + arcspace = 3; + } + if (DEBUG) + System.out.println("read fallback: " + f + " " + rowLength); + + FileInputStream fis = new FileInputStream(f); + DataInputStream dis = new DataInputStream(new BufferedInputStream(fis)); + for (int ir = 0; ir < rowLength; ir++) { + int row = rowOffset + ir * arcspace; + + for (int ic = 0; ic < rowLength; ic++) { + int col = colOffset + ic * arcspace; + + int i1 = dis.read(); // msb first! + int i0 = dis.read(); + + if (i0 == -1 || i1 == -1) + throw new RuntimeException("unexpected end of file reading hgt entry!"); + + short val = (short) ((i1 << 8) | i0); + + if (val == NODATA2) { + val = NODATA; + } + if (arcspace == 3) { + setPixel(row, col, val); + setPixel(row+1, col, val); + setPixel(row+2, col, val); + setPixel(row, col+1, val); + setPixel(row+1, col+1, val); + setPixel(row+2, col+1, val); + setPixel(row, col+2, val); + setPixel(row+1, col+2, val); + setPixel(row+2, col+2, val); + } else { + setPixel(row, col, val); + } + } + } + fis.close(); + */ + } else { + if (DEBUG) + System.out.println("none : " + filename); + } + } + + } + } + } + + // post fill zero + if (SRTM_NO_ZERO) { + for (int row = 0; row < NROWS; row++) { + for (int col = 0; col < NCOLS; col++) { + if (imagePixels[row * NCOLS + col] == 0) imagePixels[row * NCOLS + col] = NODATA; } } } @@ -126,7 +334,7 @@ public static void doConvert(String inputDir, int lonDegreeStart, int latDegreeS raster.ncols = NCOLS; raster.halfcol = halfCol5; raster.noDataValue = NODATA; - raster.cellsize = 1 / 1200.; + raster.cellsize = 1. / row_length; raster.xllcorner = lonDegreeStart - (0.5 + extraBorder) * raster.cellsize; raster.yllcorner = latDegreeStart - (0.5 + extraBorder) * raster.cellsize; raster.eval_array = imagePixels; @@ -156,6 +364,7 @@ public static void doConvert(String inputDir, int lonDegreeStart, int latDegreeS } } } + imagePixels = null; } private static String formatLon(int lon) { @@ -181,19 +390,6 @@ private static String formatLat(int lat) { return s + n.substring(n.length() - 2); } - public static void main(String[] args) throws Exception { - String filename90 = args[0]; - String filename30 = filename90.substring(0, filename90.length() - 3) + "bef"; - - int srtmLonIdx = Integer.parseInt(filename90.substring(5, 7).toLowerCase()); - int srtmLatIdx = Integer.parseInt(filename90.substring(8, 10).toLowerCase()); - - int ilon_base = (srtmLonIdx - 1) * 5 - 180; - int ilat_base = 150 - srtmLatIdx * 5 - 90; - - doConvert(args[1], ilon_base, ilat_base, filename30); - } - public SrtmRaster getRaster(File f, double lon, double lat) throws Exception { long fileSize; InputStream inputStream; @@ -231,7 +427,7 @@ public SrtmRaster getRaster(File f, double lon, double lat) throws Exception { // prefill as NODATA Arrays.fill(imagePixels, NODATA); - readHgtFromStream(inputStream, 0, 0, rowLength); + readHgtFromStream(inputStream, 0, 0, rowLength, 1); inputStream.close(); SrtmRaster raster = new SrtmRaster(); From 859b401c8b9e1a9ebe206989a27093f1e0b6e553 Mon Sep 17 00:00:00 2001 From: afischerdev Date: Fri, 27 Oct 2023 16:25:44 +0200 Subject: [PATCH 03/15] unable raster weighting --- .../src/main/java/btools/mapcreator/RasterCoder.java | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/brouter-map-creator/src/main/java/btools/mapcreator/RasterCoder.java b/brouter-map-creator/src/main/java/btools/mapcreator/RasterCoder.java index 233e3ea88..72b028cc0 100644 --- a/brouter-map-creator/src/main/java/btools/mapcreator/RasterCoder.java +++ b/brouter-map-creator/src/main/java/btools/mapcreator/RasterCoder.java @@ -45,7 +45,7 @@ public SrtmRaster decodeRaster(InputStream is) throws IOException { _decodeRaster(raster, is); - raster.usingWeights = raster.ncols > 6001; + raster.usingWeights = false; // raster.ncols > 6001; long t1 = System.currentTimeMillis(); System.out.println("finished decoding in " + (t1 - t0) + " ms ncols=" + raster.ncols + " nrows=" + raster.nrows); From 69489c2b7e3063e156e9d22cfc0633686e0347ae Mon Sep 17 00:00:00 2001 From: afischerdev Date: Fri, 27 Oct 2023 16:26:21 +0200 Subject: [PATCH 04/15] updated test --- .../src/test/java/btools/mapcreator/MapcreatorTest.java | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/brouter-map-creator/src/test/java/btools/mapcreator/MapcreatorTest.java b/brouter-map-creator/src/test/java/btools/mapcreator/MapcreatorTest.java index beb4cbd33..b509bffa2 100644 --- a/brouter-map-creator/src/test/java/btools/mapcreator/MapcreatorTest.java +++ b/brouter-map-creator/src/test/java/btools/mapcreator/MapcreatorTest.java @@ -42,7 +42,7 @@ public void mapcreatorTest() throws Exception { File unodes55 = new File(tmpdir, "unodes55"); File bordernodes = new File(tmpdir, "bordernodes.dat"); unodes55.mkdir(); - new PosUnifier().process(nodes55, unodes55, borderFile, bordernodes, workingDir.getAbsolutePath()); + new PosUnifier().process(nodes55, unodes55, borderFile, bordernodes, workingDir.getAbsolutePath(), null); // run WayLinker File segments = new File(tmpdir, "segments"); From 2be7e0c19c8a6329f28bb8d4b754bc529e7f8c6d Mon Sep 17 00:00:00 2001 From: afischerdev Date: Mon, 6 Nov 2023 12:04:47 +0100 Subject: [PATCH 05/15] rename SrmtRaster --- .../btools/mapcreator/ElevationRaster.java | 264 ++++++++++++++++++ 1 file changed, 264 insertions(+) create mode 100644 brouter-map-creator/src/main/java/btools/mapcreator/ElevationRaster.java diff --git a/brouter-map-creator/src/main/java/btools/mapcreator/ElevationRaster.java b/brouter-map-creator/src/main/java/btools/mapcreator/ElevationRaster.java new file mode 100644 index 000000000..478811ce1 --- /dev/null +++ b/brouter-map-creator/src/main/java/btools/mapcreator/ElevationRaster.java @@ -0,0 +1,264 @@ +package btools.mapcreator; + +import btools.util.ReducedMedianFilter; + +/** + * Container for a elevation raster + it's meta-data + * + * @author ab + */ +public class ElevationRaster { + public int ncols; + public int nrows; + public boolean halfcol; + public double xllcorner; + public double yllcorner; + public double cellsize; + public short[] eval_array; + public short noDataValue; + + public boolean usingWeights = false; + + private boolean missingData = false; + + public short getElevation(int ilon, int ilat) { + double lon = ilon / 1000000. - 180.; + double lat = ilat / 1000000. - 90.; + + if (usingWeights) { + return getElevationFromShiftWeights(lon, lat); + } + + // no weights calculated, use 2d linear interpolation + double dcol = (lon - xllcorner) / cellsize - 0.5; + double drow = (lat - yllcorner) / cellsize - 0.5; + int row = (int) drow; + int col = (int) dcol; + if (col < 0) col = 0; + if (col >= ncols - 1) col = ncols - 2; + if (row < 0) row = 0; + if (row >= nrows - 1) row = nrows - 2; + double wrow = drow - row; + double wcol = dcol - col; + missingData = false; + +// System.out.println( "wrow=" + wrow + " wcol=" + wcol + " row=" + row + " col=" + col ); + double eval = (1. - wrow) * (1. - wcol) * get(row, col) + + (wrow) * (1. - wcol) * get(row + 1, col) + + (1. - wrow) * (wcol) * get(row, col + 1) + + (wrow) * (wcol) * get(row + 1, col + 1); +// System.out.println( "eval=" + eval ); + return missingData ? Short.MIN_VALUE : (short) (eval * 4); + } + + private short get(int r, int c) { + short e = eval_array[(nrows - 1 - r) * ncols + c]; + if (e == Short.MIN_VALUE) missingData = true; + return e; + } + + private short getElevationFromShiftWeights(double lon, double lat) { + // calc lat-idx and -weight + double alat = lat < 0. ? -lat : lat; + alat /= 5.; + int latIdx = (int) alat; + double wlat = alat - latIdx; + + double dcol = (lon - xllcorner) / cellsize; + double drow = (lat - yllcorner) / cellsize; + int row = (int) drow; + int col = (int) dcol; + + double dgx = (dcol - col) * gridSteps; + double dgy = (drow - row) * gridSteps; + +// System.out.println( "wrow=" + wrow + " wcol=" + wcol + " row=" + row + " col=" + col ); + + int gx = (int) (dgx); + int gy = (int) (dgy); + + double wx = dgx - gx; + double wy = dgy - gy; + + double w00 = (1. - wx) * (1. - wy); + double w01 = (1. - wx) * (wy); + double w10 = (wx) * (1. - wy); + double w11 = (wx) * (wy); + + Weights[][] w0 = getWeights(latIdx); + Weights[][] w1 = getWeights(latIdx + 1); + + missingData = false; + + double m0 = w00 * getElevation(w0[gx][gy], row, col) + + w01 * getElevation(w0[gx][gy + 1], row, col) + + w10 * getElevation(w0[gx + 1][gy], row, col) + + w11 * getElevation(w0[gx + 1][gy + 1], row, col); + double m1 = w00 * getElevation(w1[gx][gy], row, col) + + w01 * getElevation(w1[gx][gy + 1], row, col) + + w10 * getElevation(w1[gx + 1][gy], row, col) + + w11 * getElevation(w1[gx + 1][gy + 1], row, col); + + if (missingData) return Short.MIN_VALUE; + double m = (1. - wlat) * m0 + wlat * m1; + return (short) (m * 2); + } + + private ReducedMedianFilter rmf = new ReducedMedianFilter(256); + + private double getElevation(Weights w, int row, int col) { + if (missingData) { + return 0.; + } + int nx = w.nx; + int ny = w.ny; + int mx = nx / 2; // mean pixels + int my = ny / 2; + + // System.out.println( "nx="+ nx + " ny=" + ny ); + + rmf.reset(); + + for (int ix = 0; ix < nx; ix++) { + for (int iy = 0; iy < ny; iy++) { + short val = get(row + iy - my, col + ix - mx); + rmf.addSample(w.getWeight(ix, iy), val); + } + } + return missingData ? 0. : rmf.calcEdgeReducedMedian(filterCenterFraction); + } + + + private static class Weights { + int nx; + int ny; + double[] weights; + long total = 0; + + Weights(int nx, int ny) { + this.nx = nx; + this.ny = ny; + weights = new double[nx * ny]; + } + + void inc(int ix, int iy) { + weights[iy * nx + ix] += 1.; + total++; + } + + void normalize(boolean verbose) { + for (int iy = 0; iy < ny; iy++) { + StringBuilder sb = verbose ? new StringBuilder() : null; + for (int ix = 0; ix < nx; ix++) { + weights[iy * nx + ix] /= total; + if (sb != null) { + int iweight = (int) (1000 * weights[iy * nx + ix] + 0.5); + String sval = " " + iweight; + sb.append(sval.substring(sval.length() - 4)); + } + } + if (sb != null) { + System.out.println(sb); + System.out.println(); + } + } + } + + double getWeight(int ix, int iy) { + return weights[iy * nx + ix]; + } + } + + private static int gridSteps = 10; + private static Weights[][][] allShiftWeights = new Weights[17][][]; + + private static double filterCenterFraction = 0.2; + private static double filterDiscRadius = 4.999; // in pixels + + static { + String sRadius = System.getProperty("filterDiscRadius"); + if (sRadius != null && sRadius.length() > 0) { + filterDiscRadius = Integer.parseInt(sRadius); + System.out.println("using filterDiscRadius = " + filterDiscRadius); + } + String sFraction = System.getProperty("filterCenterFraction"); + if (sFraction != null && sFraction.length() > 0) { + filterCenterFraction = Integer.parseInt(sFraction) / 100.; + System.out.println("using filterCenterFraction = " + filterCenterFraction); + } + } + + + // calculate interpolation weights from the overlap of a probe disc of given radius at given latitude + // ( latIndex = 0 -> 0 deg, latIndex = 16 -> 80 degree) + + private static Weights[][] getWeights(int latIndex) { + int idx = latIndex < 16 ? latIndex : 16; + + Weights[][] res = allShiftWeights[idx]; + if (res == null) { + res = calcWeights(idx); + allShiftWeights[idx] = res; + } + return res; + } + + private static Weights[][] calcWeights(int latIndex) { + double coslat = Math.cos(latIndex * 5. / 57.3); + + // radius in pixel units + double ry = filterDiscRadius; + double rx = ry / coslat; + + // gridsize is 2*radius + 1 cell + int nx = ((int) rx) * 2 + 3; + int ny = ((int) ry) * 2 + 3; + + System.out.println("nx=" + nx + " ny=" + ny); + + int mx = nx / 2; // mean pixels + int my = ny / 2; + + // create a matrix for the relative intergrid-position + + Weights[][] shiftWeights = new Weights[gridSteps + 1][]; + + // loop the intergrid-position + for (int gx = 0; gx <= gridSteps; gx++) { + shiftWeights[gx] = new Weights[gridSteps + 1]; + double x0 = mx + ((double) gx) / gridSteps; + + for (int gy = 0; gy <= gridSteps; gy++) { + double y0 = my + ((double) gy) / gridSteps; + + // create the weight-matrix + Weights weights = new Weights(nx, ny); + shiftWeights[gx][gy] = weights; + + double sampleStep = 0.001; + + for (double x = -1. + sampleStep / 2.; x < 1.; x += sampleStep) { + double mx2 = 1. - x * x; + + int x_idx = (int) (x0 + x * rx); + + for (double y = -1. + sampleStep / 2.; y < 1.; y += sampleStep) { + if (y * y > mx2) { + continue; + } + // we are in the ellipse, see what pixel we are on + int y_idx = (int) (y0 + y * ry); + weights.inc(x_idx, y_idx); + } + } + weights.normalize(true); + } + } + return shiftWeights; + } + + @Override + public String toString() { + return ncols + "," + nrows + "," + halfcol + "," + xllcorner + "," + yllcorner + "," + cellsize + "," + noDataValue + "," + usingWeights; + } +} From 50eb62361daca00b9838e3e3effccdc79542ea6f Mon Sep 17 00:00:00 2001 From: afischerdev Date: Mon, 6 Nov 2023 12:06:43 +0100 Subject: [PATCH 06/15] rename RasterCoder --- .../mapcreator/ElevationRasterCoder.java | 117 ++++++++++++++++++ 1 file changed, 117 insertions(+) create mode 100644 brouter-map-creator/src/main/java/btools/mapcreator/ElevationRasterCoder.java diff --git a/brouter-map-creator/src/main/java/btools/mapcreator/ElevationRasterCoder.java b/brouter-map-creator/src/main/java/btools/mapcreator/ElevationRasterCoder.java new file mode 100644 index 000000000..02a841dd2 --- /dev/null +++ b/brouter-map-creator/src/main/java/btools/mapcreator/ElevationRasterCoder.java @@ -0,0 +1,117 @@ +package btools.mapcreator; + +import java.io.DataInputStream; +import java.io.DataOutputStream; +import java.io.IOException; +import java.io.InputStream; +import java.io.OutputStream; + +import btools.util.MixCoderDataInputStream; +import btools.util.MixCoderDataOutputStream; + +// +// Encode/decode a raster +// + +public class ElevationRasterCoder { + public void encodeRaster(ElevationRaster raster, OutputStream os) throws IOException { + DataOutputStream dos = new DataOutputStream(os); + + long t0 = System.currentTimeMillis(); + + dos.writeInt(raster.ncols); + dos.writeInt(raster.nrows); + dos.writeBoolean(raster.halfcol); + dos.writeDouble(raster.xllcorner); + dos.writeDouble(raster.yllcorner); + dos.writeDouble(raster.cellsize); + dos.writeShort(raster.noDataValue); + + _encodeRaster(raster, os); + long t1 = System.currentTimeMillis(); + + System.out.println("finished encoding in " + (t1 - t0) + " ms"); + } + + public ElevationRaster decodeRaster(InputStream is) throws IOException { + DataInputStream dis = new DataInputStream(is); + + long t0 = System.currentTimeMillis(); + + ElevationRaster raster = new ElevationRaster(); + raster.ncols = dis.readInt(); + raster.nrows = dis.readInt(); + raster.halfcol = dis.readBoolean(); + raster.xllcorner = dis.readDouble(); + raster.yllcorner = dis.readDouble(); + raster.cellsize = dis.readDouble(); + raster.noDataValue = dis.readShort(); + raster.eval_array = new short[raster.ncols * raster.nrows]; + + _decodeRaster(raster, is); + + raster.usingWeights = false; // raster.ncols > 6001; + + long t1 = System.currentTimeMillis(); + System.out.println("finished decoding in " + (t1 - t0) + " ms ncols=" + raster.ncols + " nrows=" + raster.nrows); + return raster; + } + + + private void _encodeRaster(ElevationRaster raster, OutputStream os) throws IOException { + MixCoderDataOutputStream mco = new MixCoderDataOutputStream(os); + int nrows = raster.nrows; + int ncols = raster.ncols; + short[] pixels = raster.eval_array; + int colstep = raster.halfcol ? 2 : 1; + + for (int row = 0; row < nrows; row++) { + short lastval = Short.MIN_VALUE; // nodata + for (int col = 0; col < ncols; col += colstep) { + short val = pixels[row * ncols + col]; + if (val == -32766) { + val = lastval; // replace remaining (border) skips + } else { + lastval = val; + } + + // remap nodata + int code = val == Short.MIN_VALUE ? -1 : (val < 0 ? val - 1 : val); + mco.writeMixed(code); + } + } + mco.flush(); + } + + private void _decodeRaster(ElevationRaster raster, InputStream is) throws IOException { + MixCoderDataInputStream mci = new MixCoderDataInputStream(is); + int nrows = raster.nrows; + int ncols = raster.ncols; + short[] pixels = raster.eval_array; + int colstep = raster.halfcol ? 2 : 1; + + for (int row = 0; row < nrows; row++) { + for (int col = 0; col < ncols; col += colstep) { + int code = mci.readMixed(); + + // remap nodata + int v30 = code == -1 ? Short.MIN_VALUE : (code < 0 ? code + 1 : code); + if (raster.usingWeights && v30 > -32766) { + v30 *= 2; + } + pixels[row * ncols + col] = (short) (v30); + } + if (raster.halfcol) { + for (int col = 1; col < ncols - 1; col += colstep) { + int l = (int) pixels[row * ncols + col - 1]; + int r = (int) pixels[row * ncols + col + 1]; + short v30 = Short.MIN_VALUE; // nodata + if (l > -32766 && r > -32766) { + v30 = (short) ((l + r) / 2); + } + pixels[row * ncols + col] = v30; + } + } + } + } +} From 5198559c7795d140237282de073b18243b564ef4 Mon Sep 17 00:00:00 2001 From: afischerdev Date: Mon, 6 Nov 2023 12:48:42 +0100 Subject: [PATCH 07/15] moved converter to one file --- .../ElevationRasterTileConverter.java | 546 ++++++++++++++++++ 1 file changed, 546 insertions(+) create mode 100644 brouter-map-creator/src/main/java/btools/mapcreator/ElevationRasterTileConverter.java diff --git a/brouter-map-creator/src/main/java/btools/mapcreator/ElevationRasterTileConverter.java b/brouter-map-creator/src/main/java/btools/mapcreator/ElevationRasterTileConverter.java new file mode 100644 index 000000000..9b251fc2b --- /dev/null +++ b/brouter-map-creator/src/main/java/btools/mapcreator/ElevationRasterTileConverter.java @@ -0,0 +1,546 @@ +package btools.mapcreator; + +import java.io.BufferedInputStream; +import java.io.BufferedOutputStream; +import java.io.BufferedReader; +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.InputStreamReader; +import java.io.OutputStream; +import java.util.Arrays; +import java.util.Locale; +import java.util.StringTokenizer; +import java.util.zip.ZipEntry; +import java.util.zip.ZipInputStream; + +public class ElevationRasterTileConverter { + + public static final boolean DEBUG = false; + + public static final short NODATA2 = -32767; // hgt-formats nodata + public static final short NODATA = Short.MIN_VALUE; + + 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) + private static final int SRTM3_ROW_LENGTH = 1200; // number of elevation values per line + private static final int SRTM1_ROW_LENGTH = 3600; + private static final boolean SRTM_NO_ZERO = true; + + private int NROWS; + private int NCOLS; + private int ROW_LENGTH; + private short[] imagePixels; + + /** + * This generates elevation raster files with a 5x5 degree scope + * The output can be for 1sec (18000x18000 points) + * or for 3sec (6000x6000 points) + * When using 1sec input files a not found area can be called from 3sec pool + * The input can be 1x1 degree 1sec/3sec hgt files (also packed as zip) + * or 5x5 degree 3sec asc files (delivered as zip) + * Arguments for single file generation: + * ElevationRasterTileConverter [arc seconds (1 or 3,default=3)] [hgt-fallback-data-dir] + * Samples + * $ ... ElevationRasterTileConverter srtm_34_-1 ./srtm/hgt3sec ./srtm/srtm3_bef + * $ ... ElevationRasterTileConverter srtm_34_-1 ./srtm/hgt1sec ./srtm/srtm1_bef 1 + * $ ... ElevationRasterTileConverter srtm_34_-1 ./srtm/hgt1sec ./srtm/srtm1_bef 1 ./srtm/hgt3sec + *

+ * Arguments for multi file generation (world wide): + * $ ... ElevationRasterTileConverter all ./srtm/hgt3sec ./srtm/srtm3_bef + * $ ... ElevationRasterTileConverter all ./srtm/hgt1sec ./srtm/srtm1_bef 1 + * $ ... ElevationRasterTileConverter all ./srtm/hgt1sec ./srtm/srtm1_bef 1 ./srtm/hgt3sec + * + * @param args + * @throws Exception + */ + public static void main(String[] args) throws Exception { + if (args.length == 3 || args.length == 4 || args.length == 5) { + String filename90 = args[0]; + if ("all".equals(filename90)) { + //if (DEBUG) + System.out.println("raster convert all "); + new ElevationRasterTileConverter().doConvertAll(args[1], args[2], (args.length > 3 ? args[3] : null), (args.length == 5 ? args[4] : null)); + return; + } + // old filenames only + String filename30 = filename90 + ".bef"; //filename90.substring(0, filename90.length() - 3) + "bef"; + + int srtmLonIdx = Integer.parseInt(filename90.substring(5, 7).toLowerCase()); + int srtmLatIdx = Integer.parseInt(filename90.substring(8, 10).toLowerCase()); + + int ilon_base = (srtmLonIdx - 1) * 5 - 180; + int ilat_base = 150 - srtmLatIdx * 5 - 90; + int row_length = SRTM3_ROW_LENGTH; + String fallbackdir = null; + if (args.length > 3) { + row_length = (Integer.parseInt(args[3]) == 1 ? SRTM1_ROW_LENGTH : SRTM3_ROW_LENGTH); + fallbackdir = (args.length == 5 ? args[4] : null); + } + //if (DEBUG) + System.out.println("raster convert " + ilon_base + " " + ilat_base + " from " + srtmLonIdx + " " + srtmLatIdx + " f: " + filename90 + " rowl " + row_length); + + new ElevationRasterTileConverter().doConvert(args[1], ilon_base, ilat_base, args[2] + "/" + filename30, row_length, fallbackdir); + } else { + System.out.println("usage: java [arc seconds (1 or 3,default=3)] [hgt-fallback-data-dir]"); + System.out.println("or java all [arc seconds (1 or 3, default=3)] [hgt-fallback-data-dir]"); + return; + } + } + + private void doConvertAll(String hgtdata, String outdir, String rlen, String hgtfallbackdata) throws Exception { + int row_length = SRTM3_ROW_LENGTH; + if (rlen != null) { + row_length = (Integer.parseInt(rlen) == 1 ? SRTM1_ROW_LENGTH : SRTM3_ROW_LENGTH); + } + String filename30; + for (int ilon_base = -180; ilon_base < 180; ilon_base += 5) { + for (int ilat_base = 85; ilat_base > -90; ilat_base -= 5) { + if (PosUnifier.UseRasterRd5FileName) { + filename30 = genFilenameRd5(ilon_base, ilat_base); + } else { + filename30 = genFilenameOld(ilon_base, ilat_base); + } + if (DEBUG) + System.out.println("lidar convert all: " + filename30); + doConvert(hgtdata, ilon_base, ilat_base, outdir + "/" + filename30, row_length, hgtfallbackdata); + } + } + } + + static String genFilenameOld(int ilon_base, int ilat_base) { + int srtmLonIdx = ((ilon_base + 180) / 5) + 1; + int srtmLatIdx = (60 - ilat_base) / 5; + return String.format(Locale.US, "srtm_%02d_%02d.bef", srtmLonIdx, srtmLatIdx); + } + + static String genFilenameRd5(int ilon_base, int ilat_base) { + return String.format("srtm_%s_%s.bef", ilon_base < 0 ? "W" + (-ilon_base) : "E" + ilon_base, + ilat_base < 0 ? "S" + (-ilat_base) : "N" + ilat_base); + } + + private void readHgtZip(String filename, int rowOffset, int colOffset, int row_length, int scale) throws Exception { + ZipInputStream zis = new ZipInputStream(new BufferedInputStream(new FileInputStream(filename))); + try { + for (; ; ) { + ZipEntry ze = zis.getNextEntry(); + if (ze == null) break; + if (ze.getName().toLowerCase().endsWith(HGT_FILE_EXT)) { + readHgtFromStream(zis, rowOffset, colOffset, row_length, scale); + return; + } + } + } finally { + zis.close(); + } + } + + private void readHgtFromStream(InputStream is, int rowOffset, int colOffset, int rowLength, int scale) + throws Exception { + DataInputStream dis = new DataInputStream(new BufferedInputStream(is)); + for (int ir = 0; ir < rowLength; ir++) { + int row = rowOffset + ir * scale; + + for (int ic = 0; ic < rowLength; ic++) { + int col = colOffset + ic * scale; + + int i1 = dis.read(); // msb first! + int i0 = dis.read(); + + if (i0 == -1 || i1 == -1) + throw new RuntimeException("unexpected end of file reading hgt entry!"); + + short val = (short) ((i1 << 8) | i0); + + if (val == NODATA2) { + val = NODATA; + } + if (scale == 3) { + setPixel(row, col, val); + setPixel(row + 1, col, val); + setPixel(row + 2, col, val); + setPixel(row, col + 1, val); + setPixel(row + 1, col + 1, val); + setPixel(row + 2, col + 1, val); + setPixel(row, col + 2, val); + setPixel(row + 1, col + 2, val); + setPixel(row + 2, col + 2, val); + } else { + setPixel(row, col, val); + } + } + } + } + + private void readHgtFile(File file, int rowOffset, int colOffset, int row_length, int scale) + throws Exception { + + if (DEBUG) + System.out.println("read: " + file + " " + row_length); + + FileInputStream fis = new FileInputStream(file); + try { + readHgtFromStream(fis, rowOffset, colOffset, row_length, scale); + } finally { + fis.close(); + } + } + + /* + private void readFallbackFile(File file, int rowOffset, int colOffset, int row_length) + throws Exception { + int rowLength; + int scale; + if (file.length() > HGT_3ASEC_FILE_SIZE) { + rowLength = HGT_1ASEC_ROWS; + scale = 1; + } else { + rowLength = HGT_3ASEC_ROWS; + scale = 3; + } + if (DEBUG) + System.out.println("read fallback: " + file + " " + rowLength); + + FileInputStream fis = new FileInputStream(file); + try { + readHgtFromStream(fis, rowOffset, colOffset, rowLength, scale); + } finally { + fis.close(); + } + } + */ + + private void readAscZip(File file, ElevationRaster raster) throws Exception { + + ZipInputStream zis = new ZipInputStream(new BufferedInputStream(new FileInputStream(file))); + try { + for (; ; ) { + ZipEntry ze = zis.getNextEntry(); + if (ze.getName().endsWith(".asc")) { + readAscFromStream(zis, raster); + return; + } + } + } finally { + zis.close(); + } + } + + private String secondToken(String s) { + StringTokenizer tk = new StringTokenizer(s, " "); + tk.nextToken(); + return tk.nextToken(); + } + + private void readAscFromStream(InputStream is, ElevationRaster raster) throws Exception { + BufferedReader br = new BufferedReader(new InputStreamReader(is)); + int linenr = 0; + for (; ; ) { + linenr++; + if (linenr <= 6) { + String line = br.readLine(); + if (linenr == 1) + raster.ncols = Integer.parseInt(secondToken(line)); + else if (linenr == 2) + raster.nrows = Integer.parseInt(secondToken(line)); + else if (linenr == 3) + raster.xllcorner = Double.parseDouble(secondToken(line)); + else if (linenr == 4) + raster.yllcorner = Double.parseDouble(secondToken(line)); + else if (linenr == 5) + raster.cellsize = Double.parseDouble(secondToken(line)); + else if (linenr == 6) { + // nodata ignored here ( < -250 assumed nodata... ) + // raster.noDataValue = Short.parseShort( secondToken( line ) ); + raster.eval_array = new short[raster.ncols * raster.nrows]; + } + } else { + int row = 0; + int col = 0; + int n = 0; + boolean negative = false; + for (; ; ) { + int c = br.read(); + if (c < 0) + break; + if (c == ' ') { + if (negative) + n = -n; + short val = n < -250 ? Short.MIN_VALUE : (short) (n); + + raster.eval_array[row * raster.ncols + col] = val; + if (++col == raster.ncols) { + col = 0; + ++row; + } + n = 0; + negative = false; + } else if (c >= '0' && c <= '9') { + n = 10 * n + (c - '0'); + } else if (c == '-') { + negative = true; + } + } + break; + } + } + br.close(); + } + + + private void setPixel(int row, int col, short val) { + if (row >= 0 && row < NROWS && col >= 0 && col < NCOLS) { + imagePixels[row * NCOLS + col] = val; + } + } + + private short getPixel(int row, int col) { + if (row >= 0 && row < NROWS && col >= 0 && col < NCOLS) { + return imagePixels[row * NCOLS + col]; + } + return NODATA; + } + + + public void doConvert(String inputDir, int lonDegreeStart, int latDegreeStart, String outputFile, int row_length, String hgtfallbackdata) throws Exception { + int extraBorder = 0; + + //List foundList = new ArrayList<>(); + //List notfoundList = new ArrayList<>(); + + boolean hgtfound = false; + boolean ascfound = false; + String filename = null; + //if (row_length == SRTM1_ROW_LENGTH) + { + // check for sources w/o border + for (int latIdx = 0; latIdx < 5; latIdx++) { + int latDegree = latDegreeStart + latIdx; + + for (int lonIdx = 0; lonIdx < 5; lonIdx++) { + int lonDegree = lonDegreeStart + lonIdx; + + filename = inputDir + "/" + formatLat(latDegree) + formatLon(lonDegree) + ".zip"; + File f = new File(filename); + if (f.exists() && f.length() > 0) { + hgtfound = true; + break; + } + filename = filename.substring(0, filename.length() - 4) + ".hgt"; + f = new File(filename); + if (f.exists() && f.length() > 0) { + hgtfound = true; + break; + } + } + } + if (!hgtfound) { + filename = inputDir + "/" + genFilenameOld(lonDegreeStart, latDegreeStart).substring(0, 10) + ".zip"; + File f = new File(filename); + if (f.exists() && f.length() > 0) { + ascfound = true; + } + } + } + if (hgtfound) { // init when found + NROWS = 5 * row_length + 1 + 2 * extraBorder; + NCOLS = 5 * row_length + 1 + 2 * extraBorder; + imagePixels = new short[NROWS * NCOLS]; // 650 MB ! + + // prefill as NODATA + Arrays.fill(imagePixels, NODATA); + } else if (!ascfound) { + if (DEBUG) + System.out.println("none data: " + lonDegreeStart + " " + latDegreeStart); + return; + } + + if (hgtfound) { + for (int latIdx = -1; latIdx <= 5; latIdx++) { + int latDegree = latDegreeStart + latIdx; + int rowOffset = extraBorder + (4 - latIdx) * row_length; + + for (int lonIdx = -1; lonIdx <= 5; lonIdx++) { + int lonDegree = lonDegreeStart + lonIdx; + int colOffset = extraBorder + lonIdx * row_length; + + filename = inputDir + "/" + formatLat(latDegree) + formatLon(lonDegree) + ".zip"; + File f = new File(filename); + if (f.exists() && f.length() > 0) { + if (DEBUG) + System.out.println("exist: " + filename); + readHgtZip(filename, rowOffset, colOffset, row_length + 1, 1); + continue; + } + filename = filename.substring(0, filename.length() - 4) + ".hgt"; + f = new File(filename); + if (f.exists() && f.length() > 0) { + if (DEBUG) + System.out.println("exist: " + filename); + readHgtFile(f, rowOffset, colOffset, row_length + 1, 1); + continue; + } else { + if (hgtfallbackdata != null) { + filename = hgtfallbackdata + "/" + formatLat(latDegree) + formatLon(lonDegree) + ".hgt"; + f = new File(filename); + if (f.exists() && f.length() > 0) { + readHgtFile(f, rowOffset, colOffset, SRTM3_ROW_LENGTH + 1, 3); + continue; + } + filename = filename.substring(0, filename.length() - 4) + ".zip"; + f = new File(filename); + if (f.exists() && f.length() > 0) { + readHgtZip(filename, rowOffset, colOffset, SRTM3_ROW_LENGTH + 1, 3); + } else { + if (DEBUG) + System.out.println("none : " + filename); + } + } + + } + } + } + + // post fill zero + if (SRTM_NO_ZERO) { + for (int row = 0; row < NROWS; row++) { + for (int col = 0; col < NCOLS; col++) { + if (imagePixels[row * NCOLS + col] == 0) imagePixels[row * NCOLS + col] = NODATA; + } + } + } + } + + boolean halfCol5 = false; // no halfcol tiles in lidar data (?) + + + ElevationRaster raster = new ElevationRaster(); + if (hgtfound) { + raster.nrows = NROWS; + raster.ncols = NCOLS; + raster.halfcol = halfCol5; + raster.noDataValue = NODATA; + raster.cellsize = 1. / row_length; + raster.xllcorner = lonDegreeStart - (0.5 + extraBorder) * raster.cellsize; + raster.yllcorner = latDegreeStart - (0.5 + extraBorder) * raster.cellsize; + raster.eval_array = imagePixels; + } + + if (ascfound) { + File f = new File(filename); + readAscZip(f, raster); + } + + // encode the raster + OutputStream os = new BufferedOutputStream(new FileOutputStream(outputFile)); + new ElevationRasterCoder().encodeRaster(raster, os); + os.close(); + + // decode the raster + InputStream is = new BufferedInputStream(new FileInputStream(outputFile)); + ElevationRaster raster2 = new ElevationRasterCoder().decodeRaster(is); + is.close(); + + short[] pix2 = raster2.eval_array; + if (pix2.length != raster.eval_array.length) + throw new RuntimeException("length mismatch!"); + + // compare decoding result + for (int row = 0; row < raster.nrows; row++) { + int colstep = halfCol5 ? 2 : 1; + for (int col = 0; col < raster.ncols; col += colstep) { + int idx = row * raster.ncols + col; + short p2 = pix2[idx]; + if (p2 != raster.eval_array[idx]) { + throw new RuntimeException("content mismatch: p2=" + p2 + " p1=" + raster.eval_array[idx]); + } + } + } + imagePixels = null; + } + + private static String formatLon(int lon) { + if (lon >= 180) + lon -= 180; // TODO: w180 oder E180 ? + + String s = "E"; + if (lon < 0) { + lon = -lon; + s = "W"; + } + String n = "000" + lon; + return s + n.substring(n.length() - 3); + } + + private static String formatLat(int lat) { + String s = "N"; + if (lat < 0) { + lat = -lat; + s = "S"; + } + String n = "00" + lat; + return s + n.substring(n.length() - 2); + } + + + public ElevationRaster 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))); + 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; + } + } + } else { + fileSize = f.length(); + inputStream = new FileInputStream(f); + } + + 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, 1); + inputStream.close(); + + ElevationRaster raster = new ElevationRaster(); + raster.nrows = NROWS; + raster.ncols = NCOLS; + raster.halfcol = false; // assume full resolution + raster.noDataValue = NODATA; + 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; + } + + +} From 36d692da84d96d8454fa2c7876ea3fd63452c65c Mon Sep 17 00:00:00 2001 From: afischerdev Date: Mon, 6 Nov 2023 12:50:50 +0100 Subject: [PATCH 08/15] set new raster calls --- .../java/btools/mapcreator/PosUnifier.java | 40 +++++++++---------- 1 file changed, 20 insertions(+), 20 deletions(-) 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 10e660554..e9a43fe3d 100644 --- a/brouter-map-creator/src/main/java/btools/mapcreator/PosUnifier.java +++ b/brouter-map-creator/src/main/java/btools/mapcreator/PosUnifier.java @@ -7,7 +7,6 @@ import java.io.FileInputStream; import java.io.InputStream; import java.util.HashMap; -import java.util.Locale; import java.util.Map; import btools.util.CompactLongSet; @@ -24,17 +23,17 @@ */ public class PosUnifier extends MapCreatorBase { - public static final boolean UseLidarRd5FileName = false; + public static final boolean UseRasterRd5FileName = false; private DiffCoderDataOutputStream nodesOutStream; private DiffCoderDataOutputStream borderNodesOut; private File nodeTilesOut; private CompactLongSet[] positionSets; - private Map srtmmap; + private Map srtmmap; private int lastSrtmLonIdx; private int lastSrtmLatIdx; - private SrtmRaster lastSrtmRaster; + private ElevationRaster lastSrtmRaster; private String srtmdir; private String srtmfallbackdir; @@ -51,7 +50,7 @@ public static void main(String[] args) throws Exception { NodeData n = new NodeData(1, lon, lat); short selev = Short.MIN_VALUE; - SrtmRaster srtm = null; + ElevationRaster srtm = null; /* // check hgt direct srtm = posu.hgtForNode(n.ilon, n.ilat); @@ -60,7 +59,7 @@ public static void main(String[] args) throws Exception { } else { System.out.println("hgtForNode no data"); } - posu.resetSrtm(); + posu.resetElevationRaster(); System.out.println("-----> selv for hgt " + lat + ", " + lon + " = " + selev + " = " + (selev / 4.)); srtm = null; selev = Short.MIN_VALUE; @@ -69,7 +68,7 @@ public static void main(String[] args) throws Exception { srtm = posu.srtmForNode(n.ilon, n.ilat); } if (srtm != null) selev = srtm.getElevation(n.ilon, n.ilat); - posu.resetSrtm(); + posu.resetElevationRaster(); System.out.println("-----> selv for bef " + lat + ", " + lon + " = " + selev + " = " + (selev / 4.)); return; } else if (args.length != 5 && args.length != 6) { @@ -107,7 +106,7 @@ public void process(File nodeTilesIn, File nodeTilesOut, File bordernidsinfile, @Override public void nodeFileStart(File nodefile) throws Exception { - resetSrtm(); + resetElevationRaster(); nodesOutStream = createOutStream(fileFromTemplate(nodefile, nodeTilesOut, "u5d")); @@ -125,7 +124,7 @@ public void nextNode(NodeData n) throws Exception { srtm = srtmForNode(n.ilon, n.ilat); } */ - SrtmRaster srtm = srtmForNode(n.ilon, n.ilat); + ElevationRaster srtm = srtmForNode(n.ilon, n.ilat); if (srtm != null) n.selev = srtm.getElevation(n.ilon, n.ilat); findUniquePos(n); @@ -139,7 +138,7 @@ public void nextNode(NodeData n) throws Exception { @Override public void nodeFileEnd(File nodeFile) throws Exception { nodesOutStream.close(); - resetSrtm(); + resetElevationRaster(); } private boolean checkAdd(int lon, int lat) { @@ -187,7 +186,7 @@ private void _findUniquePos(NodeData n) { * get the srtm data set for a position srtm coords are * srtm__ where srtmLon = 180 + lon, srtmLat = 60 - lat */ - private SrtmRaster srtmForNode(int ilon, int ilat) throws Exception { + private ElevationRaster srtmForNode(int ilon, int ilat) throws Exception { int srtmLonIdx = (ilon + 5000000) / 5000000; int srtmLatIdx = (654999999 - ilat) / 5000000 - 100; // ugly negative rounding... @@ -198,7 +197,7 @@ private SrtmRaster srtmForNode(int ilon, int ilat) throws Exception { lastSrtmLatIdx = srtmLatIdx; String filename; - if (UseLidarRd5FileName) { + if (UseRasterRd5FileName) { filename = genFilenameRd5(ilon, ilat); } else { filename = genFilenameXY(srtmLonIdx, srtmLatIdx); @@ -210,7 +209,7 @@ private SrtmRaster srtmForNode(int ilon, int ilat) throws Exception { if (f.exists()) { try { InputStream isc = new BufferedInputStream(new FileInputStream(f)); - lastSrtmRaster = new RasterCoder().decodeRaster(isc); + lastSrtmRaster = new ElevationRasterCoder().decodeRaster(isc); isc.close(); } catch (Exception e) { System.out.println("**** ERROR reading " + f + " ****"); @@ -225,7 +224,7 @@ private SrtmRaster srtmForNode(int ilon, int ilat) throws Exception { try { InputStream isc = new BufferedInputStream(new FileInputStream(f)); //lastSrtmRaster = new StatRasterCoder().decodeRaster(isc); - lastSrtmRaster = new RasterCoder().decodeRaster(isc); + lastSrtmRaster = new ElevationRasterCoder().decodeRaster(isc); isc.close(); } catch (Exception e) { System.out.println("**** ERROR reading " + f + " ****"); @@ -258,23 +257,23 @@ static String genFilenameRd5(int ilon, int ilat) { } - private SrtmRaster hgtForNode(int ilon, int ilat) throws Exception { + private ElevationRaster hgtForNode(int ilon, int ilat) throws Exception { double lon = (ilon - 180000000) / 1000000.; double lat = (ilat - 90000000) / 1000000.; String filename = buildHgtFilename(lat, lon); // don't block lastSrtmRaster - SrtmRaster srtm = srtmmap.get(filename); + ElevationRaster srtm = srtmmap.get(filename); if (srtm == null) { File f = new File(new File(srtmdir), filename + ".zip"); if (f.exists()) { - srtm = new ConvertLidarTile().getRaster(f, lon, lat); + srtm = new ElevationRasterTileConverter().getRaster(f, lon, lat); srtmmap.put(filename, srtm); return srtm; } f = new File(new File(srtmdir), filename + ".hgt"); if (f.exists()) { - srtm = new ConvertLidarTile().getRaster(f, lon, lat); + srtm = new ElevationRasterTileConverter().getRaster(f, lon, lat); srtmmap.put(filename, srtm); return srtm; } @@ -282,6 +281,7 @@ private SrtmRaster hgtForNode(int ilon, int ilat) throws Exception { return srtm; } + private String buildHgtFilename(double llat, double llon) { int lat = (int) llat; int lon = (int) llon; @@ -297,10 +297,10 @@ private String buildHgtFilename(double llat, double llon) { lon = -lon + 1; } - return String.format(Locale.US, "%s%02d%s%03d", latPref, lat, lonPref, lon); + return String.format("%s%02d%s%03d", latPref, lat, lonPref, lon); } - private void resetSrtm() { + private void resetElevationRaster() { srtmmap = new HashMap<>(); lastSrtmLonIdx = -1; lastSrtmLatIdx = -1; From fbad69474644853fb86f6d5072a54333d069323c Mon Sep 17 00:00:00 2001 From: afischerdev Date: Mon, 6 Nov 2023 12:55:50 +0100 Subject: [PATCH 09/15] removed unused files --- .../btools/mapcreator/ConvertLidarTile.java | 446 ------------------ .../btools/mapcreator/ConvertSrtmTile.java | 260 ---------- .../btools/mapcreator/ConvertUrlList.java | 53 --- .../java/btools/mapcreator/RasterCoder.java | 112 ----- .../main/java/btools/mapcreator/SrtmData.java | 166 ------- .../java/btools/mapcreator/SrtmRaster.java | 264 ----------- 6 files changed, 1301 deletions(-) delete mode 100644 brouter-map-creator/src/main/java/btools/mapcreator/ConvertLidarTile.java delete mode 100644 brouter-map-creator/src/main/java/btools/mapcreator/ConvertSrtmTile.java delete mode 100644 brouter-map-creator/src/main/java/btools/mapcreator/ConvertUrlList.java delete mode 100644 brouter-map-creator/src/main/java/btools/mapcreator/RasterCoder.java delete mode 100644 brouter-map-creator/src/main/java/btools/mapcreator/SrtmData.java delete mode 100644 brouter-map-creator/src/main/java/btools/mapcreator/SrtmRaster.java diff --git a/brouter-map-creator/src/main/java/btools/mapcreator/ConvertLidarTile.java b/brouter-map-creator/src/main/java/btools/mapcreator/ConvertLidarTile.java deleted file mode 100644 index cab1d0b7f..000000000 --- a/brouter-map-creator/src/main/java/btools/mapcreator/ConvertLidarTile.java +++ /dev/null @@ -1,446 +0,0 @@ -package btools.mapcreator; - -import java.io.BufferedInputStream; -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.ArrayList; -import java.util.Arrays; -import java.util.List; -import java.util.zip.ZipEntry; -import java.util.zip.ZipInputStream; - -public class ConvertLidarTile { - - public static final boolean DEBUG = false; - - public static final short NODATA2 = -32767; // hgt-formats nodata - public static final short NODATA = Short.MIN_VALUE; - - 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) - private static final int SRTM3_ROW_LENGTH = 1200; // number of elevation values per line - private static final int SRTM1_ROW_LENGTH = 3600; - private static final boolean SRTM_NO_ZERO = true; - - private int NROWS; - private int NCOLS; - private int ROW_LENGTH; - private short[] imagePixels; - - - public static void main(String[] args) throws Exception { - if (args.length == 3 || args.length == 4 || args.length == 5) { - String filename90 = args[0]; - if ("all".equals(filename90)) { - //if (DEBUG) - System.out.println("lidar convert all "); - new ConvertLidarTile().doConvertAll(args[1], args[2], (args.length > 3 ? args[3] : null), (args.length == 5 ? args[4] : null)); - return; - } - // old filenames only - String filename30 = filename90 + ".bef"; //filename90.substring(0, filename90.length() - 3) + "bef"; - - int srtmLonIdx = Integer.parseInt(filename90.substring(5, 7).toLowerCase()); - int srtmLatIdx = Integer.parseInt(filename90.substring(8, 10).toLowerCase()); - - int ilon_base = (srtmLonIdx - 1) * 5 - 180; - int ilat_base = 150 - srtmLatIdx * 5 - 90; - int row_length = SRTM3_ROW_LENGTH; - String fallbackdir = null; - if (args.length > 3) { - row_length = (Integer.parseInt(args[3]) == 1 ? SRTM1_ROW_LENGTH : SRTM3_ROW_LENGTH); - fallbackdir = (args.length == 5 ? args[4] : null); - } - //if (DEBUG) - System.out.println("lidar convert " + ilon_base + " " + ilat_base + " from " + srtmLonIdx + " " + srtmLatIdx + " f: " + filename90 + " rowl " + row_length); - - new ConvertLidarTile().doConvert(args[1], ilon_base, ilat_base, args[2] + "/" + filename30, row_length, fallbackdir); - } else { - System.out.println("usage: java [arc seconds (1 or 3,default=3)] [hgt-fallback-data-dir]"); - System.out.println("or java all [arc seconds (1 or 3, default=3)] [hgt-fallback-data-dir]"); - return; - } - } - - private void doConvertAll(String hgtdata, String outdir, String rlen, String hgtfallbackdata) throws Exception { - int row_length = SRTM3_ROW_LENGTH; - if (rlen != null) { - row_length = (Integer.parseInt(rlen) == 1 ? SRTM1_ROW_LENGTH : SRTM3_ROW_LENGTH); - } - String filename30; - for (int ilon_base = -180; ilon_base < 180; ilon_base += 5) { - for (int ilat_base = 85; ilat_base > -90; ilat_base -= 5) { - if (PosUnifier.UseLidarRd5FileName) { - filename30 = genFilenameRd5(ilon_base, ilat_base); - } else { - filename30 = genFilenameOld(ilon_base, ilat_base); - } - if (DEBUG) - System.out.println("lidar convert all: " + filename30); - doConvert(hgtdata, ilon_base, ilat_base, outdir + "/" + filename30, row_length, hgtfallbackdata); - } - } - } - - static String genFilenameOld(int ilon_base, int ilat_base) { - int srtmLonIdx = ((ilon_base + 180) / 5) + 1; - int srtmLatIdx = (60 - ilat_base) / 5; - return String.format("srtm_%02d_%02d.bef", srtmLonIdx, srtmLatIdx); - } - - static String genFilenameRd5(int ilon_base, int ilat_base) { - return String.format("srtm_%s_%s.bef", ilon_base < 0 ? "W" + (-ilon_base) : "E" + ilon_base, - ilat_base < 0 ? "S" + (-ilat_base) : "N" + ilat_base); - } - - private void readHgtZip(String filename, int rowOffset, int colOffset, int row_length) throws Exception { - ZipInputStream zis = new ZipInputStream(new BufferedInputStream(new FileInputStream(filename))); - try { - for (; ; ) { - ZipEntry ze = zis.getNextEntry(); - if (ze == null) break; - if (ze.getName().toLowerCase().endsWith(HGT_FILE_EXT)) { - readHgtFromStream(zis, rowOffset, colOffset, row_length, 1); - return; - } - } - } finally { - zis.close(); - } - } - - private void readHgtFromStream(InputStream is, int rowOffset, int colOffset, int rowLength, int scale) - throws Exception { - DataInputStream dis = new DataInputStream(new BufferedInputStream(is)); - for (int ir = 0; ir < rowLength; ir++) { - int row = rowOffset + ir * scale; - - for (int ic = 0; ic < rowLength; ic++) { - int col = colOffset + ic * scale; - - int i1 = dis.read(); // msb first! - int i0 = dis.read(); - - if (i0 == -1 || i1 == -1) - throw new RuntimeException("unexpected end of file reading hgt entry!"); - - short val = (short) ((i1 << 8) | i0); - - if (val == NODATA2) { - val = NODATA; - } - if (scale == 3) { - setPixel(row, col, val); - setPixel(row + 1, col, val); - setPixel(row + 2, col, val); - setPixel(row, col + 1, val); - setPixel(row + 1, col + 1, val); - setPixel(row + 2, col + 1, val); - setPixel(row, col + 2, val); - setPixel(row + 1, col + 2, val); - setPixel(row + 2, col + 2, val); - } else { - setPixel(row, col, val); - } - } - } - } - - private void readFallbackFile(File file, int rowOffset, int colOffset, int row_length) - throws Exception { - int rowLength; - int scale; - if (file.length() > HGT_3ASEC_FILE_SIZE) { - rowLength = HGT_1ASEC_ROWS; - scale = 1; - } else { - rowLength = HGT_3ASEC_ROWS; - scale = 3; - } - if (DEBUG) - System.out.println("read fallback: " + file + " " + rowLength); - - FileInputStream fis = new FileInputStream(file); - try { - readHgtFromStream(fis, rowOffset, colOffset, rowLength, scale); - } finally { - fis.close(); - } - } - - private void setPixel(int row, int col, short val) { - if (row >= 0 && row < NROWS && col >= 0 && col < NCOLS) { - imagePixels[row * NCOLS + col] = val; - } - } - - private short getPixel(int row, int col) { - if (row >= 0 && row < NROWS && col >= 0 && col < NCOLS) { - return imagePixels[row * NCOLS + col]; - } - return NODATA; - } - - - public void doConvert(String inputDir, int lonDegreeStart, int latDegreeStart, String outputFile, int row_length, String hgtfallbackdata) throws Exception { - int extraBorder = 0; - - List foundList = new ArrayList<>(); - List notfoundList = new ArrayList<>(); - - boolean found = false; - - if (row_length == SRTM1_ROW_LENGTH) { - // check for sources w/o border - for (int latIdx = 0; latIdx < 5; latIdx++) { - int latDegree = latDegreeStart + latIdx; - - for (int lonIdx = 0; lonIdx < 5; lonIdx++) { - int lonDegree = lonDegreeStart + lonIdx; - - String filename = inputDir + "/" + formatLat(latDegree) + formatLon(lonDegree) + ".zip"; - File f = new File(filename); - if (f.exists() && f.length() > 0) { - found = true; - break; - } - } - } - } else { - // ignore when srtm3 - found = true; - } - if (found) { // init when found - NROWS = 5 * row_length + 1 + 2 * extraBorder; - NCOLS = 5 * row_length + 1 + 2 * extraBorder; - imagePixels = new short[NROWS * NCOLS]; // 650 MB ! - - // prefill as NODATA - Arrays.fill(imagePixels, NODATA); - } else { - if (DEBUG) - System.out.println("none 1sec data: " + lonDegreeStart + " " + latDegreeStart); - return; - } - - for (int latIdx = -1; latIdx <= 5; latIdx++) { - int latDegree = latDegreeStart + latIdx; - int rowOffset = extraBorder + (4 - latIdx) * row_length; - - for (int lonIdx = -1; lonIdx <= 5; lonIdx++) { - int lonDegree = lonDegreeStart + lonIdx; - int colOffset = extraBorder + lonIdx * row_length; - - String filename = inputDir + "/" + formatLat(latDegree) + formatLon(lonDegree) + ".zip"; - File f = new File(filename); - if (f.exists() && f.length() > 0) { - if (DEBUG) - System.out.println("exist: " + filename); - readHgtZip(filename, rowOffset, colOffset, row_length + 1); - } else { - if (hgtfallbackdata != null) { - String filenamehgt = hgtfallbackdata + "/" + formatLat(latDegree) + formatLon(lonDegree) + ".hgt"; - f = new File(filenamehgt); - if (f.exists() && f.length() > 0) { - readFallbackFile(f, rowOffset, colOffset, row_length + 1); - /*if (imagePixels == null) { - imagePixels = new short[NROWS * NCOLS]; - Arrays.fill(imagePixels, NODATA); - found = true; - } - */ - /* - int rowLength; - int arcspace; - if (f.length() > HGT_3ASEC_FILE_SIZE) { - rowLength = HGT_1ASEC_ROWS; - arcspace = 1; - } else { - rowLength = HGT_3ASEC_ROWS; - arcspace = 3; - } - if (DEBUG) - System.out.println("read fallback: " + f + " " + rowLength); - - FileInputStream fis = new FileInputStream(f); - DataInputStream dis = new DataInputStream(new BufferedInputStream(fis)); - for (int ir = 0; ir < rowLength; ir++) { - int row = rowOffset + ir * arcspace; - - for (int ic = 0; ic < rowLength; ic++) { - int col = colOffset + ic * arcspace; - - int i1 = dis.read(); // msb first! - int i0 = dis.read(); - - if (i0 == -1 || i1 == -1) - throw new RuntimeException("unexpected end of file reading hgt entry!"); - - short val = (short) ((i1 << 8) | i0); - - if (val == NODATA2) { - val = NODATA; - } - if (arcspace == 3) { - setPixel(row, col, val); - setPixel(row+1, col, val); - setPixel(row+2, col, val); - setPixel(row, col+1, val); - setPixel(row+1, col+1, val); - setPixel(row+2, col+1, val); - setPixel(row, col+2, val); - setPixel(row+1, col+2, val); - setPixel(row+2, col+2, val); - } else { - setPixel(row, col, val); - } - } - } - fis.close(); - */ - } else { - if (DEBUG) - System.out.println("none : " + filename); - } - } - - } - } - } - - // post fill zero - if (SRTM_NO_ZERO) { - for (int row = 0; row < NROWS; row++) { - for (int col = 0; col < NCOLS; col++) { - if (imagePixels[row * NCOLS + col] == 0) imagePixels[row * NCOLS + col] = NODATA; - } - } - } - - 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. / row_length; - raster.xllcorner = lonDegreeStart - (0.5 + extraBorder) * raster.cellsize; - raster.yllcorner = latDegreeStart - (0.5 + extraBorder) * raster.cellsize; - raster.eval_array = imagePixels; - - // encode the raster - OutputStream os = new BufferedOutputStream(new FileOutputStream(outputFile)); - new RasterCoder().encodeRaster(raster, os); - os.close(); - - // decode the raster - InputStream is = new BufferedInputStream(new FileInputStream(outputFile)); - SrtmRaster raster2 = new RasterCoder().decodeRaster(is); - is.close(); - - short[] pix2 = raster2.eval_array; - if (pix2.length != imagePixels.length) - throw new RuntimeException("length mismatch!"); - - // compare decoding result - for (int row = 0; row < NROWS; row++) { - int colstep = halfCol5 ? 2 : 1; - for (int col = 0; col < NCOLS; col += colstep) { - int idx = row * NCOLS + col; - short p2 = pix2[idx]; - if (p2 != imagePixels[idx]) { - throw new RuntimeException("content mismatch: p2=" + p2 + " p1=" + imagePixels[idx]); - } - } - } - imagePixels = null; - } - - private static String formatLon(int lon) { - if (lon >= 180) - lon -= 180; // TODO: w180 oder E180 ? - - String s = "E"; - if (lon < 0) { - lon = -lon; - s = "W"; - } - String n = "000" + lon; - return s + n.substring(n.length() - 3); - } - - private static String formatLat(int lat) { - String s = "N"; - if (lat < 0) { - lat = -lat; - s = "S"; - } - String n = "00" + lat; - return s + n.substring(n.length() - 2); - } - - 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))); - 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; - } - } - } else { - fileSize = f.length(); - inputStream = new FileInputStream(f); - } - - 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, 1); - inputStream.close(); - - SrtmRaster raster = new SrtmRaster(); - raster.nrows = NROWS; - raster.ncols = NCOLS; - raster.halfcol = false; // assume full resolution - raster.noDataValue = NODATA; - 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/ConvertSrtmTile.java b/brouter-map-creator/src/main/java/btools/mapcreator/ConvertSrtmTile.java deleted file mode 100644 index e18c66451..000000000 --- a/brouter-map-creator/src/main/java/btools/mapcreator/ConvertSrtmTile.java +++ /dev/null @@ -1,260 +0,0 @@ -package btools.mapcreator; - -import java.io.*; -import java.util.zip.*; - -public class ConvertSrtmTile { - public static int NROWS; - public static int NCOLS; - - public static final short SKIPDATA = -32766; // >50 degree skipped pixel - public static final short NODATA2 = -32767; // bil-formats nodata - public static final short NODATA = Short.MIN_VALUE; - - static short[] imagePixels; - - public static int[] diffs = new int[100]; - - private static void readBilZip(String filename, int rowOffset, int colOffset, boolean halfCols) throws Exception { - ZipInputStream zis = new ZipInputStream(new BufferedInputStream(new FileInputStream(filename))); - try { - for (; ; ) { - ZipEntry ze = zis.getNextEntry(); - if (ze.getName().endsWith(".bil")) { - readBilFromStream(zis, rowOffset, colOffset, halfCols); - return; - } - } - } finally { - zis.close(); - } - } - - private static void readBilFromStream(InputStream is, int rowOffset, int colOffset, boolean halfCols) - throws Exception { - DataInputStream dis = new DataInputStream(new BufferedInputStream(is)); - for (int ir = 0; ir < 3601; ir++) { - int row = rowOffset + ir; - - for (int ic = 0; ic < 3601; ic++) { - int col = colOffset + ic; - - if ((ic % 2) == 1 && halfCols) { - if (getPixel(row, col) == NODATA) { - setPixel(row, col, SKIPDATA); - } - continue; - } - - int i0 = dis.read(); - int i1 = dis.read(); - - if (i0 == -1 || i1 == -1) - throw new RuntimeException("unexcepted end of file reading bil entry!"); - - short val = (short) ((i1 << 8) | i0); - - if (val == NODATA2) { - val = NODATA; - } - - setPixel(row, col, val); - } - } - } - - - private static void setPixel(int row, int col, short val) { - if (row >= 0 && row < NROWS && col >= 0 && col < NCOLS) { - imagePixels[row * NCOLS + col] = val; - } - } - - private static short getPixel(int row, int col) { - if (row >= 0 && row < NROWS && col >= 0 && col < NCOLS) { - return imagePixels[row * NCOLS + col]; - } - return NODATA; - } - - - public static void doConvert(String inputDir, String v1Dir, int lonDegreeStart, int latDegreeStart, String outputFile, SrtmRaster raster90) throws Exception { - int extraBorder = 10; - int datacells = 0; - int mismatches = 0; - - NROWS = 5 * 3600 + 1 + 2 * extraBorder; - NCOLS = 5 * 3600 + 1 + 2 * extraBorder; - - 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; - } - } - - for (int latIdx = -1; latIdx <= 5; latIdx++) { - int latDegree = latDegreeStart + latIdx; - int rowOffset = extraBorder + (4 - latIdx) * 3600; - - for (int lonIdx = -1; lonIdx <= 5; lonIdx++) { - int lonDegree = lonDegreeStart + lonIdx; - int colOffset = extraBorder + lonIdx * 3600; - - String filename = inputDir + "/" + formatLat(latDegree) + "_" + formatLon(lonDegree) + "_1arc_v3_bil.zip"; - File f = new File(filename); - if (f.exists() && f.length() > 0) { - System.out.println("exist: " + filename); - boolean halfCol = latDegree >= 50 || latDegree < -50; - readBilZip(filename, rowOffset, colOffset, halfCol); - } else { - System.out.println("none : " + filename); - } - } - } - - boolean halfCol5 = latDegreeStart >= 50 || latDegreeStart < -50; - - for (int row90 = 0; row90 < 6001; row90++) { - int crow = 3 * row90 + extraBorder; // center row of 3x3 - for (int col90 = 0; col90 < 6001; col90++) { - int ccol = 3 * col90 + extraBorder; // center col of 3x3 - - // evaluate 3x3 area - if (raster90 != null && (!halfCol5 || (col90 % 2) == 0)) { - short v90 = raster90.eval_array[row90 * 6001 + col90]; - - int sum = 0; - int nodatas = 0; - int datas = 0; - int colstep = halfCol5 ? 2 : 1; - for (int row = crow - 1; row <= crow + 1; row++) { - for (int col = ccol - colstep; col <= ccol + colstep; col += colstep) { - short v30 = imagePixels[row * NCOLS + col]; - if (v30 == NODATA) { - nodatas++; - } else if (v30 != SKIPDATA) { - sum += v30; - datas++; - } - } - } - boolean doReplace = nodatas > 0 || v90 == NODATA || datas < 7; - if (!doReplace) { - datacells++; - int diff = sum - datas * v90; - if (diff < -4 || diff > 4) { - doReplace = true; - mismatches++; - } - - if (diff > -50 && diff < 50 && (row90 % 1200) != 0 && (col90 % 1200) != 0) { - diffs[diff + 50]++; - } - } - if (doReplace) { - for (int row = crow - 1; row <= crow + 1; row++) { - for (int col = ccol - colstep; col <= ccol + colstep; col += colstep) { - imagePixels[row * NCOLS + col] = v90; - } - } - } - } - } - } - - SrtmRaster raster = new SrtmRaster(); - raster.nrows = NROWS; - raster.ncols = NCOLS; - raster.halfcol = halfCol5; - raster.noDataValue = NODATA; - raster.cellsize = 1 / 3600.; - raster.xllcorner = lonDegreeStart - (0.5 + extraBorder) * raster.cellsize; - raster.yllcorner = latDegreeStart - (0.5 + extraBorder) * raster.cellsize; - raster.eval_array = imagePixels; - - // encode the raster - OutputStream os = new BufferedOutputStream(new FileOutputStream(outputFile)); - new RasterCoder().encodeRaster(raster, os); - os.close(); - - // decode the raster - InputStream is = new BufferedInputStream(new FileInputStream(outputFile)); - SrtmRaster raster2 = new RasterCoder().decodeRaster(is); - is.close(); - - short[] pix2 = raster2.eval_array; - if (pix2.length != imagePixels.length) - throw new RuntimeException("length mismatch!"); - - // compare decoding result - for (int row = 0; row < NROWS; row++) { - int colstep = halfCol5 ? 2 : 1; - for (int col = 0; col < NCOLS; col += colstep) { - int idx = row * NCOLS + col; - if (imagePixels[idx] == SKIPDATA) { - continue; - } - short p2 = pix2[idx]; - if (p2 > SKIPDATA) { - p2 /= 2; - } - if (p2 != imagePixels[idx]) { - throw new RuntimeException("content mismatch!"); - } - } - } - - for (int i = 1; i < 100; i++) System.out.println("diff[" + (i - 50) + "] = " + diffs[i]); - System.out.println("datacells=" + datacells + " mismatch%=" + (100. * mismatches) / datacells); - btools.util.MixCoderDataOutputStream.stats(); - // test( raster ); - // raster.calcWeights( 50. ); - // test( raster ); - // 39828330 &lon=3115280&layer=OpenStreetMap - } - - private static void test(SrtmRaster raster) { - int lat0 = 39828330; - int lon0 = 3115280; - - for (int iy = -9; iy <= 9; iy++) { - StringBuilder sb = new StringBuilder(); - for (int ix = -9; ix <= 9; ix++) { - int lat = lat0 + 90000000 - 100 * iy; - int lon = lon0 + 180000000 + 100 * ix; - int ival = (int) (raster.getElevation(lon, lat) / 4.); - String sval = " " + ival; - sb.append(sval.substring(sval.length() - 4)); - } - System.out.println(sb); - System.out.println(); - } - } - - private static String formatLon(int lon) { - if (lon >= 180) - lon -= 180; // TODO: w180 oder E180 ? - - String s = "e"; - if (lon < 0) { - lon = -lon; - s = "w"; - } - String n = "000" + lon; - return s + n.substring(n.length() - 3); - } - - private static String formatLat(int lat) { - String s = "n"; - if (lat < 0) { - lat = -lat; - s = "s"; - } - String n = "00" + lat; - return s + n.substring(n.length() - 2); - } - -} diff --git a/brouter-map-creator/src/main/java/btools/mapcreator/ConvertUrlList.java b/brouter-map-creator/src/main/java/btools/mapcreator/ConvertUrlList.java deleted file mode 100644 index ab5b8dbef..000000000 --- a/brouter-map-creator/src/main/java/btools/mapcreator/ConvertUrlList.java +++ /dev/null @@ -1,53 +0,0 @@ -package btools.mapcreator; - -import java.io.BufferedReader; -import java.io.File; -import java.io.FileReader; - -public class ConvertUrlList { - public static final short NODATA = -32767; - - public static void main(String[] args) throws Exception { - BufferedReader br = new BufferedReader(new FileReader(args[0])); - - for (; ; ) { - String line = br.readLine(); - if (line == null) { - break; - } - int idx1 = line.indexOf("srtm_"); - if (idx1 < 0) { - continue; - } - - String filename90 = line.substring(idx1); - String filename30 = filename90.substring(0, filename90.length() - 3) + "bef"; - - if (new File(filename30).exists()) { - continue; - } - - // int srtmLonIdx = (ilon+5000000)/5000000; -> ilon = (srtmLonIdx-1)*5 - // int srtmLatIdx = (154999999-ilat)/5000000; -> ilat = 155 - srtmLatIdx*5 - - int srtmLonIdx = Integer.parseInt(filename90.substring(5, 7).toLowerCase()); - int srtmLatIdx = Integer.parseInt(filename90.substring(8, 10).toLowerCase()); - - int ilon_base = (srtmLonIdx - 1) * 5 - 180; - int ilat_base = 150 - srtmLatIdx * 5 - 90; - - SrtmRaster raster90 = null; - - File file90 = new File(new File(args[1]), filename90); - if (file90.exists()) { - System.out.println("reading " + file90); - raster90 = new SrtmData(file90).getRaster(); - } - - ConvertSrtmTile.doConvert(args[2], args[3], ilon_base, ilat_base, filename30, raster90); - } - br.close(); - } - - -} diff --git a/brouter-map-creator/src/main/java/btools/mapcreator/RasterCoder.java b/brouter-map-creator/src/main/java/btools/mapcreator/RasterCoder.java deleted file mode 100644 index 72b028cc0..000000000 --- a/brouter-map-creator/src/main/java/btools/mapcreator/RasterCoder.java +++ /dev/null @@ -1,112 +0,0 @@ -package btools.mapcreator; - -import java.io.*; - -import btools.util.*; - -// -// Encode/decode a raster -// - -public class RasterCoder { - public void encodeRaster(SrtmRaster raster, OutputStream os) throws IOException { - DataOutputStream dos = new DataOutputStream(os); - - long t0 = System.currentTimeMillis(); - - dos.writeInt(raster.ncols); - dos.writeInt(raster.nrows); - dos.writeBoolean(raster.halfcol); - dos.writeDouble(raster.xllcorner); - dos.writeDouble(raster.yllcorner); - dos.writeDouble(raster.cellsize); - dos.writeShort(raster.noDataValue); - - _encodeRaster(raster, os); - long t1 = System.currentTimeMillis(); - - System.out.println("finished encoding in " + (t1 - t0) + " ms"); - } - - public SrtmRaster decodeRaster(InputStream is) throws IOException { - DataInputStream dis = new DataInputStream(is); - - long t0 = System.currentTimeMillis(); - - SrtmRaster raster = new SrtmRaster(); - raster.ncols = dis.readInt(); - raster.nrows = dis.readInt(); - raster.halfcol = dis.readBoolean(); - raster.xllcorner = dis.readDouble(); - raster.yllcorner = dis.readDouble(); - raster.cellsize = dis.readDouble(); - raster.noDataValue = dis.readShort(); - raster.eval_array = new short[raster.ncols * raster.nrows]; - - _decodeRaster(raster, is); - - raster.usingWeights = false; // raster.ncols > 6001; - - long t1 = System.currentTimeMillis(); - System.out.println("finished decoding in " + (t1 - t0) + " ms ncols=" + raster.ncols + " nrows=" + raster.nrows); - return raster; - } - - - private void _encodeRaster(SrtmRaster raster, OutputStream os) throws IOException { - MixCoderDataOutputStream mco = new MixCoderDataOutputStream(os); - int nrows = raster.nrows; - int ncols = raster.ncols; - short[] pixels = raster.eval_array; - int colstep = raster.halfcol ? 2 : 1; - - for (int row = 0; row < nrows; row++) { - short lastval = Short.MIN_VALUE; // nodata - for (int col = 0; col < ncols; col += colstep) { - short val = pixels[row * ncols + col]; - if (val == -32766) { - val = lastval; // replace remaining (border) skips - } else { - lastval = val; - } - - // remap nodata - int code = val == Short.MIN_VALUE ? -1 : (val < 0 ? val - 1 : val); - mco.writeMixed(code); - } - } - mco.flush(); - } - - private void _decodeRaster(SrtmRaster raster, InputStream is) throws IOException { - MixCoderDataInputStream mci = new MixCoderDataInputStream(is); - int nrows = raster.nrows; - int ncols = raster.ncols; - short[] pixels = raster.eval_array; - int colstep = raster.halfcol ? 2 : 1; - - for (int row = 0; row < nrows; row++) { - for (int col = 0; col < ncols; col += colstep) { - int code = mci.readMixed(); - - // remap nodata - int v30 = code == -1 ? Short.MIN_VALUE : (code < 0 ? code + 1 : code); - if (raster.usingWeights && v30 > -32766) { - v30 *= 2; - } - pixels[row * ncols + col] = (short) (v30); - } - if (raster.halfcol) { - for (int col = 1; col < ncols - 1; col += colstep) { - int l = (int) pixels[row * ncols + col - 1]; - int r = (int) pixels[row * ncols + col + 1]; - short v30 = Short.MIN_VALUE; // nodata - if (l > -32766 && r > -32766) { - v30 = (short) ((l + r) / 2); - } - pixels[row * ncols + col] = v30; - } - } - } - } -} diff --git a/brouter-map-creator/src/main/java/btools/mapcreator/SrtmData.java b/brouter-map-creator/src/main/java/btools/mapcreator/SrtmData.java deleted file mode 100644 index b101cd889..000000000 --- a/brouter-map-creator/src/main/java/btools/mapcreator/SrtmData.java +++ /dev/null @@ -1,166 +0,0 @@ -package btools.mapcreator; - -/** - * This is a wrapper for a 5*5 degree srtm file in ascii/zip-format - *

- * - filter out unused nodes according to the way file - * - enhance with SRTM elevation data - * - split further in smaller (5*5 degree) tiles - * - * @author ab - */ - -import java.io.BufferedInputStream; -import java.io.BufferedOutputStream; -import java.io.BufferedReader; -import java.io.File; -import java.io.FileInputStream; -import java.io.FileOutputStream; -import java.io.InputStream; -import java.io.InputStreamReader; -import java.io.OutputStream; -import java.util.StringTokenizer; -import java.util.zip.ZipEntry; -import java.util.zip.ZipInputStream; - -public class SrtmData { - private SrtmRaster raster; - - public SrtmData(File file) throws Exception { - raster = new SrtmRaster(); - - ZipInputStream zis = new ZipInputStream(new BufferedInputStream(new FileInputStream(file))); - try { - for (; ; ) { - ZipEntry ze = zis.getNextEntry(); - if (ze.getName().endsWith(".asc")) { - readFromStream(zis); - return; - } - } - } finally { - zis.close(); - } - } - - public SrtmRaster getRaster() { - return raster; - } - - private String secondToken(String s) { - StringTokenizer tk = new StringTokenizer(s, " "); - tk.nextToken(); - return tk.nextToken(); - } - - public void readFromStream(InputStream is) throws Exception { - BufferedReader br = new BufferedReader(new InputStreamReader(is)); - int linenr = 0; - for (; ; ) { - linenr++; - if (linenr <= 6) { - String line = br.readLine(); - if (linenr == 1) - raster.ncols = Integer.parseInt(secondToken(line)); - else if (linenr == 2) - raster.nrows = Integer.parseInt(secondToken(line)); - else if (linenr == 3) - raster.xllcorner = Double.parseDouble(secondToken(line)); - else if (linenr == 4) - raster.yllcorner = Double.parseDouble(secondToken(line)); - else if (linenr == 5) - raster.cellsize = Double.parseDouble(secondToken(line)); - else if (linenr == 6) { - // nodata ignored here ( < -250 assumed nodata... ) - // raster.noDataValue = Short.parseShort( secondToken( line ) ); - raster.eval_array = new short[raster.ncols * raster.nrows]; - } - } else { - int row = 0; - int col = 0; - int n = 0; - boolean negative = false; - for (; ; ) { - int c = br.read(); - if (c < 0) - break; - if (c == ' ') { - if (negative) - n = -n; - short val = n < -250 ? Short.MIN_VALUE : (short) (n); - - raster.eval_array[row * raster.ncols + col] = val; - if (++col == raster.ncols) { - col = 0; - ++row; - } - n = 0; - negative = false; - } else if (c >= '0' && c <= '9') { - n = 10 * n + (c - '0'); - } else if (c == '-') { - negative = true; - } - } - break; - } - } - br.close(); - } - - public static void main(String[] args) throws Exception { - String fromDir = args[0]; - String toDir = args[1]; - - File[] files = new File(fromDir).listFiles(); - for (File f : files) { - if (!f.getName().endsWith(".zip")) { - continue; - } - System.out.println("*** reading: " + f); - long t0 = System.currentTimeMillis(); - SrtmRaster raster = new SrtmData(f).getRaster(); - long t1 = System.currentTimeMillis(); - String name = f.getName(); - - long zipTime = t1 - t0; - - File fbef = new File(new File(toDir), name.substring(0, name.length() - 3) + "bef"); - System.out.println("recoding: " + f + " to " + fbef); - OutputStream osbef = new BufferedOutputStream(new FileOutputStream(fbef)); - new RasterCoder().encodeRaster(raster, osbef); - osbef.close(); - - System.out.println("*** re-reading: " + fbef); - - long t2 = System.currentTimeMillis(); - InputStream isc = new BufferedInputStream(new FileInputStream(fbef)); - SrtmRaster raster2 = new RasterCoder().decodeRaster(isc); - isc.close(); - long t3 = System.currentTimeMillis(); - - long befTime = t3 - t2; - - System.out.println("*** zip-time: " + zipTime + "*** bef-time: " + befTime); - - String s1 = raster.toString(); - String s2 = raster2.toString(); - - if (!s1.equals(s2)) { - throw new IllegalArgumentException("missmatch: " + s1 + "<--->" + s2); - } - - int cols = raster.ncols; - int rows = raster.nrows; - for (int c = 0; c < cols; c++) { - for (int r = 0; r < rows; r++) { - int idx = r * cols + c; - - if (raster.eval_array[idx] != raster2.eval_array[idx]) { - throw new IllegalArgumentException("missmatch: at " + c + "," + r + ": " + raster.eval_array[idx] + "<--->" + raster2.eval_array[idx]); - } - } - } - } - } -} diff --git a/brouter-map-creator/src/main/java/btools/mapcreator/SrtmRaster.java b/brouter-map-creator/src/main/java/btools/mapcreator/SrtmRaster.java deleted file mode 100644 index 3a307bfe1..000000000 --- a/brouter-map-creator/src/main/java/btools/mapcreator/SrtmRaster.java +++ /dev/null @@ -1,264 +0,0 @@ -package btools.mapcreator; - -import btools.util.ReducedMedianFilter; - -/** - * Container for a srtm-raster + it's meta-data - * - * @author ab - */ -public class SrtmRaster { - public int ncols; - public int nrows; - public boolean halfcol; - public double xllcorner; - public double yllcorner; - public double cellsize; - public short[] eval_array; - public short noDataValue; - - public boolean usingWeights = false; - - private boolean missingData = false; - - public short getElevation(int ilon, int ilat) { - double lon = ilon / 1000000. - 180.; - double lat = ilat / 1000000. - 90.; - - if (usingWeights) { - return getElevationFromShiftWeights(lon, lat); - } - - // no weights calculated, use 2d linear interpolation - double dcol = (lon - xllcorner) / cellsize - 0.5; - double drow = (lat - yllcorner) / cellsize - 0.5; - int row = (int) drow; - int col = (int) dcol; - if (col < 0) col = 0; - if (col >= ncols - 1) col = ncols - 2; - if (row < 0) row = 0; - if (row >= nrows - 1) row = nrows - 2; - double wrow = drow - row; - double wcol = dcol - col; - missingData = false; - -// System.out.println( "wrow=" + wrow + " wcol=" + wcol + " row=" + row + " col=" + col ); - double eval = (1. - wrow) * (1. - wcol) * get(row, col) - + (wrow) * (1. - wcol) * get(row + 1, col) - + (1. - wrow) * (wcol) * get(row, col + 1) - + (wrow) * (wcol) * get(row + 1, col + 1); -// System.out.println( "eval=" + eval ); - return missingData ? Short.MIN_VALUE : (short) (eval * 4); - } - - private short get(int r, int c) { - short e = eval_array[(nrows - 1 - r) * ncols + c]; - if (e == Short.MIN_VALUE) missingData = true; - return e; - } - - private short getElevationFromShiftWeights(double lon, double lat) { - // calc lat-idx and -weight - double alat = lat < 0. ? -lat : lat; - alat /= 5.; - int latIdx = (int) alat; - double wlat = alat - latIdx; - - double dcol = (lon - xllcorner) / cellsize; - double drow = (lat - yllcorner) / cellsize; - int row = (int) drow; - int col = (int) dcol; - - double dgx = (dcol - col) * gridSteps; - double dgy = (drow - row) * gridSteps; - -// System.out.println( "wrow=" + wrow + " wcol=" + wcol + " row=" + row + " col=" + col ); - - int gx = (int) (dgx); - int gy = (int) (dgy); - - double wx = dgx - gx; - double wy = dgy - gy; - - double w00 = (1. - wx) * (1. - wy); - double w01 = (1. - wx) * (wy); - double w10 = (wx) * (1. - wy); - double w11 = (wx) * (wy); - - Weights[][] w0 = getWeights(latIdx); - Weights[][] w1 = getWeights(latIdx + 1); - - missingData = false; - - double m0 = w00 * getElevation(w0[gx][gy], row, col) - + w01 * getElevation(w0[gx][gy + 1], row, col) - + w10 * getElevation(w0[gx + 1][gy], row, col) - + w11 * getElevation(w0[gx + 1][gy + 1], row, col); - double m1 = w00 * getElevation(w1[gx][gy], row, col) - + w01 * getElevation(w1[gx][gy + 1], row, col) - + w10 * getElevation(w1[gx + 1][gy], row, col) - + w11 * getElevation(w1[gx + 1][gy + 1], row, col); - - if (missingData) return Short.MIN_VALUE; - double m = (1. - wlat) * m0 + wlat * m1; - return (short) (m * 2); - } - - private ReducedMedianFilter rmf = new ReducedMedianFilter(256); - - private double getElevation(Weights w, int row, int col) { - if (missingData) { - return 0.; - } - int nx = w.nx; - int ny = w.ny; - int mx = nx / 2; // mean pixels - int my = ny / 2; - - // System.out.println( "nx="+ nx + " ny=" + ny ); - - rmf.reset(); - - for (int ix = 0; ix < nx; ix++) { - for (int iy = 0; iy < ny; iy++) { - short val = get(row + iy - my, col + ix - mx); - rmf.addSample(w.getWeight(ix, iy), val); - } - } - return missingData ? 0. : rmf.calcEdgeReducedMedian(filterCenterFraction); - } - - - private static class Weights { - int nx; - int ny; - double[] weights; - long total = 0; - - Weights(int nx, int ny) { - this.nx = nx; - this.ny = ny; - weights = new double[nx * ny]; - } - - void inc(int ix, int iy) { - weights[iy * nx + ix] += 1.; - total++; - } - - void normalize(boolean verbose) { - for (int iy = 0; iy < ny; iy++) { - StringBuilder sb = verbose ? new StringBuilder() : null; - for (int ix = 0; ix < nx; ix++) { - weights[iy * nx + ix] /= total; - if (sb != null) { - int iweight = (int) (1000 * weights[iy * nx + ix] + 0.5); - String sval = " " + iweight; - sb.append(sval.substring(sval.length() - 4)); - } - } - if (sb != null) { - System.out.println(sb); - System.out.println(); - } - } - } - - double getWeight(int ix, int iy) { - return weights[iy * nx + ix]; - } - } - - private static int gridSteps = 10; - private static Weights[][][] allShiftWeights = new Weights[17][][]; - - private static double filterCenterFraction = 0.2; - private static double filterDiscRadius = 4.999; // in pixels - - static { - String sRadius = System.getProperty("filterDiscRadius"); - if (sRadius != null && sRadius.length() > 0) { - filterDiscRadius = Integer.parseInt(sRadius); - System.out.println("using filterDiscRadius = " + filterDiscRadius); - } - String sFraction = System.getProperty("filterCenterFraction"); - if (sFraction != null && sFraction.length() > 0) { - filterCenterFraction = Integer.parseInt(sFraction) / 100.; - System.out.println("using filterCenterFraction = " + filterCenterFraction); - } - } - - - // calculate interpolation weights from the overlap of a probe disc of given radius at given latitude - // ( latIndex = 0 -> 0 deg, latIndex = 16 -> 80 degree) - - private static Weights[][] getWeights(int latIndex) { - int idx = latIndex < 16 ? latIndex : 16; - - Weights[][] res = allShiftWeights[idx]; - if (res == null) { - res = calcWeights(idx); - allShiftWeights[idx] = res; - } - return res; - } - - private static Weights[][] calcWeights(int latIndex) { - double coslat = Math.cos(latIndex * 5. / 57.3); - - // radius in pixel units - double ry = filterDiscRadius; - double rx = ry / coslat; - - // gridsize is 2*radius + 1 cell - int nx = ((int) rx) * 2 + 3; - int ny = ((int) ry) * 2 + 3; - - System.out.println("nx=" + nx + " ny=" + ny); - - int mx = nx / 2; // mean pixels - int my = ny / 2; - - // create a matrix for the relative intergrid-position - - Weights[][] shiftWeights = new Weights[gridSteps + 1][]; - - // loop the intergrid-position - for (int gx = 0; gx <= gridSteps; gx++) { - shiftWeights[gx] = new Weights[gridSteps + 1]; - double x0 = mx + ((double) gx) / gridSteps; - - for (int gy = 0; gy <= gridSteps; gy++) { - double y0 = my + ((double) gy) / gridSteps; - - // create the weight-matrix - Weights weights = new Weights(nx, ny); - shiftWeights[gx][gy] = weights; - - double sampleStep = 0.001; - - for (double x = -1. + sampleStep / 2.; x < 1.; x += sampleStep) { - double mx2 = 1. - x * x; - - int x_idx = (int) (x0 + x * rx); - - for (double y = -1. + sampleStep / 2.; y < 1.; y += sampleStep) { - if (y * y > mx2) { - continue; - } - // we are in the ellipse, see what pixel we are on - int y_idx = (int) (y0 + y * ry); - weights.inc(x_idx, y_idx); - } - } - weights.normalize(true); - } - } - return shiftWeights; - } - - @Override - public String toString() { - return ncols + "," + nrows + "," + halfcol + "," + xllcorner + "," + yllcorner + "," + cellsize + "," + noDataValue + "," + usingWeights; - } -} From cae367025f966c0ea93f2074b5415bebba3c1fa2 Mon Sep 17 00:00:00 2001 From: afischerdev Date: Mon, 6 Nov 2023 17:48:10 +0100 Subject: [PATCH 10/15] added an elevation raster generation part --- docs/developers/build_segments.md | 36 ++++++++++++++++++++++++++++++- 1 file changed, 35 insertions(+), 1 deletion(-) diff --git a/docs/developers/build_segments.md b/docs/developers/build_segments.md index 513fab5a0..79bd2967b 100644 --- a/docs/developers/build_segments.md +++ b/docs/developers/build_segments.md @@ -26,7 +26,7 @@ official BRouter segments files are the ones provided by [CGIAR](https://cgiarcsi.community/data/srtm-90m-digital-elevation-database-v4-1/). If you are working with rather small geographical extracts, you can download tiles manually using [this -interface](http://srtm.csi.cgiar.org/SELECTION/inputCoord.asp) (use the +interface](https://srtm.csi.cgiar.org/srtmdata/) (use the "ArcInfo ASCII" format), instead of having to ask for an access for bulk download of data. There is no need to unzip the downloaded files, the `process_pbf_planet.sh` script expects a folder with all the ZIP files inside @@ -48,3 +48,37 @@ and set the `PLANET_FILE` variable to point to it. _Note:_ It is possible that you encounter an error complaining about not being able to run `bash^M` on Linux/Mac OS. You can fix this one by running `sed -i -e 's/\r$//' process_pbf_planet.sh`. + + +## Run a generation for elevation data tiles + +To match the 5x5 OSM data grid (*.rd5) files from BRouter, there are elevation +data in a 5x5 degree format (*.bef). At the moment (end of 2023) the naming of +this elevation tiles follows the konvention used by srtm.csi.cgiar.org: srtm_x_y + +As the srtm files are only available between 60N and 60S the filenames above 60N +contains negative values. e.g. srtm_34_-1 as a tile above srtm_34_00. + +Please see OSM wiki for more info on [srtm](https://wiki.openstreetmap.org/wiki/SRTM). + +The converter generates bef tiles from `hgt` files, `zipped hgt` files and `zipped 'ESRI' asc` files. + +Converter call with arguments for a single tile generation: + +`ElevationRasterTileConverter [arc seconds (1 or 3,default=3)] [hgt-fallback-data-dir] +Samples: +$ ... ElevationRasterTileConverter srtm_34_-1 ./srtm/hgt3sec ./srtm/srtm3_bef +$ ... ElevationRasterTileConverter srtm_34_-1 ./srtm/hgt1sec ./srtm/srtm1_bef 1 +$ ... ElevationRasterTileConverter srtm_34_-1 ./srtm/hgt1sec ./srtm/srtm1_bef 1 ./srtm/hgt3sec +` +Arguments for multi file generation (world wide): + +`$ ... ElevationRasterTileConverter all ./srtm/hgt3sec ./srtm/srtm3_bef +$ ... ElevationRasterTileConverter all ./srtm/hgt1sec ./srtm/srtm1_bef 1 +$ ... ElevationRasterTileConverter all ./srtm/hgt1sec ./srtm/srtm1_bef 1 ./srtm/hgt3sec +` + +To use 1sec and 3sec bef tiles at rd5 generation time you need an extra parameter to the fallback folder. +E.g. +`$ ... PosUnifier nodes55 unodes55 bordernids.dat bordernodes.dat ../srtm/srtm1_bef ../srtm/srtm3_bef +` From 477c675d466c9bf3d0c774bdc04f1795afb695f6 Mon Sep 17 00:00:00 2001 From: afischerdev Date: Tue, 28 Nov 2023 15:04:14 +0100 Subject: [PATCH 11/15] prepare elevation type in rd5 --- .../src/main/java/btools/mapcreator/PosUnifier.java | 11 +++++++++-- .../src/main/java/btools/mapcreator/WayLinker.java | 11 ++++++++++- 2 files changed, 19 insertions(+), 3 deletions(-) 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 e9a43fe3d..cea8cb5ad 100644 --- a/brouter-map-creator/src/main/java/btools/mapcreator/PosUnifier.java +++ b/brouter-map-creator/src/main/java/btools/mapcreator/PosUnifier.java @@ -28,6 +28,7 @@ public class PosUnifier extends MapCreatorBase { private DiffCoderDataOutputStream nodesOutStream; private DiffCoderDataOutputStream borderNodesOut; private File nodeTilesOut; + private File outNodeFile; private CompactLongSet[] positionSets; private Map srtmmap; @@ -107,8 +108,8 @@ public void process(File nodeTilesIn, File nodeTilesOut, File bordernidsinfile, @Override public void nodeFileStart(File nodefile) throws Exception { resetElevationRaster(); - - nodesOutStream = createOutStream(fileFromTemplate(nodefile, nodeTilesOut, "u5d")); + outNodeFile = fileFromTemplate(nodefile, nodeTilesOut, "u5d"); + nodesOutStream = createOutStream(outNodeFile); positionSets = new CompactLongSet[2500]; } @@ -138,6 +139,12 @@ public void nextNode(NodeData n) throws Exception { @Override public void nodeFileEnd(File nodeFile) throws Exception { nodesOutStream.close(); + if (outNodeFile != null) { + if (lastSrtmRaster != null) { + String newName = outNodeFile.getAbsolutePath() + (lastSrtmRaster.nrows > 6001 ? "_1": "_3"); + outNodeFile.renameTo(new File(newName)); + } + } resetElevationRaster(); } diff --git a/brouter-map-creator/src/main/java/btools/mapcreator/WayLinker.java b/brouter-map-creator/src/main/java/btools/mapcreator/WayLinker.java index 91c70a969..0a01b9328 100644 --- a/brouter-map-creator/src/main/java/btools/mapcreator/WayLinker.java +++ b/brouter-map-creator/src/main/java/btools/mapcreator/WayLinker.java @@ -51,6 +51,7 @@ public class WayLinker extends MapCreatorBase implements Runnable { private short lookupMinorVersion; private long creationTimeStamp; + private byte elevationType; private BExpressionContextWay expctxWay; @@ -217,7 +218,14 @@ public boolean wayFileStart(File wayfile) throws Exception { File trafficFile = fileFromTemplate(wayfile, trafficTilesIn, "trf"); // process corresponding node-file, if any - File nodeFile = fileFromTemplate(wayfile, nodeTilesIn, "u5d"); + elevationType = 3; + File nodeFile = fileFromTemplate(wayfile, nodeTilesIn, "u5d_1"); + if (nodeFile.exists()) { + elevationType = 1; + } else { + nodeFile = fileFromTemplate(wayfile, nodeTilesIn, "u5d_3"); + if (!nodeFile.exists()) nodeFile = fileFromTemplate(wayfile, nodeTilesIn, "u5d"); + } if (nodeFile.exists()) { reset(); @@ -539,6 +547,7 @@ public void wayFileEnd(File wayfile) throws Exception { for (int i55 = 0; i55 < 25; i55++) { os.writeInt(fileHeaderCrcs[i55]); } + os.writeByte(elevationType); os.close(); From 16d019c1d09e2e060f860c497a938bf62c86b694 Mon Sep 17 00:00:00 2001 From: afischerdev Date: Tue, 28 Nov 2023 15:14:13 +0100 Subject: [PATCH 12/15] read elevation type from rd5 --- .../main/java/btools/mapaccess/NodesCache.java | 16 ++++++++++++++++ .../src/main/java/btools/mapaccess/OsmFile.java | 2 ++ .../main/java/btools/mapaccess/PhysicalFile.java | 8 ++++++++ 3 files changed, 26 insertions(+) diff --git a/brouter-mapaccess/src/main/java/btools/mapaccess/NodesCache.java b/brouter-mapaccess/src/main/java/btools/mapaccess/NodesCache.java index f2a63160d..e16b8d264 100644 --- a/brouter-mapaccess/src/main/java/btools/mapaccess/NodesCache.java +++ b/brouter-mapaccess/src/main/java/btools/mapaccess/NodesCache.java @@ -385,4 +385,20 @@ public void close() { } } } + + public int getElevationType(int ilon, int ilat) { + int lonDegree = ilon / 1000000; + int latDegree = ilat / 1000000; + OsmFile[] fileRow = fileRows[latDegree]; + int ndegrees = fileRow == null ? 0 : fileRow.length; + for (int i = 0; i < ndegrees; i++) { + if (fileRow[i].lonDegree == lonDegree) { + OsmFile osmf = fileRow[i]; + if (osmf != null) return osmf.elevationType; + break; + } + } + return 3; + } + } diff --git a/brouter-mapaccess/src/main/java/btools/mapaccess/OsmFile.java b/brouter-mapaccess/src/main/java/btools/mapaccess/OsmFile.java index 0f773cf20..67eff299f 100644 --- a/brouter-mapaccess/src/main/java/btools/mapaccess/OsmFile.java +++ b/brouter-mapaccess/src/main/java/btools/mapaccess/OsmFile.java @@ -33,6 +33,7 @@ final class OsmFile { private int cellsize; private int ncaches; private int indexsize; + protected byte elevationType = 3; public OsmFile(PhysicalFile rafile, int lonDegree, int latDegree, DataBuffers dataBuffers) throws IOException { this.lonDegree = lonDegree; @@ -43,6 +44,7 @@ public OsmFile(PhysicalFile rafile, int lonDegree, int latDegree, DataBuffers da if (rafile != null) { divisor = rafile.divisor; + elevationType = rafile.elevationType; cellsize = 1000000 / divisor; ncaches = divisor * divisor; diff --git a/brouter-mapaccess/src/main/java/btools/mapaccess/PhysicalFile.java b/brouter-mapaccess/src/main/java/btools/mapaccess/PhysicalFile.java index 2c2a18620..e61f46965 100644 --- a/brouter-mapaccess/src/main/java/btools/mapaccess/PhysicalFile.java +++ b/brouter-mapaccess/src/main/java/btools/mapaccess/PhysicalFile.java @@ -24,6 +24,7 @@ final public class PhysicalFile { String fileName; public int divisor = 80; + public byte elevationType = 3; public static void main(String[] args) { MicroCache.debug = true; @@ -113,6 +114,10 @@ public PhysicalFile(File f, DataBuffers dataBuffers, int lookupVersion, int look if (len == pos) return; // old format o.k. + if ((len-pos) > extraLen) { + extraLen++; + } + if (len < pos + extraLen) { // > is o.k. for future extensions! throw new IOException("file of size " + len + " too short, should be " + (pos + extraLen)); } @@ -134,5 +139,8 @@ public PhysicalFile(File f, DataBuffers dataBuffers, int lookupVersion, int look for (int i = 0; i < 25; i++) { fileHeaderCrcs[i] = dis.readInt(); } + try { + elevationType = dis.readByte(); + } catch (Exception e) {} } } From 068a5ff714e6924ac0acc3a9c116e42c0b9dc009 Mon Sep 17 00:00:00 2001 From: afischerdev Date: Tue, 28 Nov 2023 15:17:51 +0100 Subject: [PATCH 13/15] use elevation type for filter value --- .../java/btools/router/RoutingEngine.java | 21 +++++++++++++++++-- 1 file changed, 19 insertions(+), 2 deletions(-) diff --git a/brouter-core/src/main/java/btools/router/RoutingEngine.java b/brouter-core/src/main/java/btools/router/RoutingEngine.java index baff3dca4..dd7a2af4b 100644 --- a/brouter-core/src/main/java/btools/router/RoutingEngine.java +++ b/brouter-core/src/main/java/btools/router/RoutingEngine.java @@ -20,6 +20,7 @@ import btools.mapaccess.OsmLinkHolder; import btools.mapaccess.OsmNode; import btools.mapaccess.OsmNodePairSet; +import btools.mapaccess.OsmPos; import btools.util.CompactLongMap; import btools.util.SortedHeap; import btools.util.StackSampler; @@ -855,11 +856,12 @@ private void recalcTrack(OsmTrack t) { if (ele_last != Short.MIN_VALUE) { ehb = ehb + (ele_last - ele) * eleFactor; } + double filter = elevationFilter(n); if (ehb > 0) { ascend += ehb; ehb = 0; - } else if (ehb < -10) { - ehb = -10; + } else if (ehb < filter) { + ehb = filter; } } @@ -896,6 +898,21 @@ private void recalcTrack(OsmTrack t) { logInfo("filtered ascend = " + t.ascend); } + /** + * find the elevation type for position + * to determine the filter value + * + * @param n the point + * @return the filter value for 1sec / 3sec elevation source + */ + double elevationFilter(OsmPos n) { + if (nodesCache != null) { + int r = nodesCache.getElevationType(n.getILon(), n.getILat()); + if (r == 1) return -5.; + } + return -10.; + } + // geometric position matching finding the nearest routable way-section private void matchWaypointsToNodes(List unmatchedWaypoints) { resetCache(false); From 526bb53b703a3a187825dbcc2994e6fca74f2a4d Mon Sep 17 00:00:00 2001 From: afischerdev Date: Sat, 2 Mar 2024 10:44:38 +0100 Subject: [PATCH 14/15] 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 000000000..50b8d8140 --- /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 000000000..f161d7dc4 --- /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); + + } + } +} From dcc9719ba61136019799e0d49dc18dc374d62944 Mon Sep 17 00:00:00 2001 From: afischerdev Date: Sun, 3 Mar 2024 17:52:23 +0100 Subject: [PATCH 15/15] added create elev image with diff colors --- .../CreateElevationRasterImage.java | 89 ++++++++++++++++--- 1 file changed, 79 insertions(+), 10 deletions(-) diff --git a/brouter-map-creator/src/main/java/btools/mapcreator/CreateElevationRasterImage.java b/brouter-map-creator/src/main/java/btools/mapcreator/CreateElevationRasterImage.java index 50b8d8140..935965f03 100644 --- a/brouter-map-creator/src/main/java/btools/mapcreator/CreateElevationRasterImage.java +++ b/brouter-map-creator/src/main/java/btools/mapcreator/CreateElevationRasterImage.java @@ -4,9 +4,11 @@ import java.awt.image.BufferedImage; import java.awt.image.DataBufferInt; import java.io.BufferedInputStream; +import java.io.BufferedReader; import java.io.File; import java.io.FileInputStream; import java.io.FileOutputStream; +import java.io.FileReader; import java.io.InputStream; import java.util.HashMap; import java.util.Map; @@ -25,26 +27,40 @@ public class CreateElevationRasterImage { Map srtmmap; int lastSrtmLonIdx; int lastSrtmLatIdx; + short maxElev = Short.MIN_VALUE; + short minElev = Short.MAX_VALUE; 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 { + + private void createImage(double lon, double lat, String dir, String imageName, int maxX, int maxY, int downscale, String format, String colors) throws Exception { srtmdir = dir; + if (colors != null) { + loadColors(colors); + } if (format.equals("hgt")) { createImageFromHgt(lon, lat, dir, imageName, maxX, maxY); return; } + if (!format.equals("bef")) { + System.out.println("wrong format (bef|hgt)"); + return; + } srtmmap = new HashMap<>(); lastSrtmLonIdx = -1; lastSrtmLatIdx = -1; lastSrtmRaster = null; NodeData n = new NodeData(1, lon, lat); ElevationRaster srtm = srtmForNode(n.ilon, n.ilat); + if (srtm == null) { + System.out.println("no data"); + return; + } 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)); + //System.out.println("srtm pos " + Math.round(pos[0]) + " " + Math.round(pos[1])); short[] raster = srtm.eval_array; int rasterX = srtm.ncols; int rasterY = srtm.nrows; @@ -53,11 +69,12 @@ private void createImage(double lon, double lat, String dir, String imageName, i 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); + //short e = getElevationXY(srtm, pos[0] + (sizeY - y) * downscale, pos[1] + (x * downscale)); + short e = get(srtm, (int) Math.round(pos[0]) + (sizeY - y), x + (int) Math.round(pos[1])); + if (e != Short.MIN_VALUE && e < minElev) minElev = e; + if (e != Short.MIN_VALUE && e > maxElev) maxElev = e; if (e == Short.MIN_VALUE) { imgraster[sizeY * y + x] = 0xffff; @@ -67,6 +84,8 @@ private void createImage(double lon, double lat, String dir, String imageName, i } } } + System.out.println("srtm target " + sizeX + " " + sizeY + " (" + rasterX + " " + rasterY + ")" + " min " + minElev + " max " + maxElev); + if (DEBUG) { String out = "short "; for (int i = 0; i < 100; i++) { @@ -97,9 +116,13 @@ private void createImage(double lon, double lat, String dir, String imageName, i 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); + if (data == null) { + System.out.println("no data"); + return; + } + 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]; @@ -107,15 +130,19 @@ private void createImageFromHgt(double lon, double lat, String dir, String image 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 && e < minElev) minElev = e; + if (e != HgtReader.HGT_VOID && e > maxElev) maxElev = e; + 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.); + imgraster[sizeY * y + x] = getColorForHeight((short) (e)); } } } + System.out.println("hgt size " + rowlen + " x " + rowlen + " min " + minElev + " max " + maxElev); if (DEBUG) { String out = "short "; for (int i = 0; i < 100; i++) { @@ -143,6 +170,47 @@ private void createImageFromHgt(double lon, double lat, String dir, String image } + private void loadColors(String colors) { + if (DEBUG) System.out.println("colors=" + colors); + File colFile = new File(colors); + if (colFile.exists()) { + BufferedReader reader = null; + colorMap = new TreeMap<>(); + try { + reader = new BufferedReader(new FileReader(colors)); + String line = reader.readLine(); + + while (line != null) { + if (DEBUG) System.out.println(line); + String[] sa = line.split(","); + if (!line.startsWith("#") && sa.length == 4) { + short e = Short.valueOf(sa[0].trim()); + short r = Short.valueOf(sa[1].trim()); + short g = Short.valueOf(sa[2].trim()); + short b = Short.valueOf(sa[3].trim()); + colorMap.put(e, new Color(r, g, b)); + } + // read next line + line = reader.readLine(); + } + + } catch (Exception e) { + e.printStackTrace(); + colorMap = null; + } finally { + if (reader != null) { + try { + reader.close(); + } catch (Exception e) { + e.printStackTrace(); + } + } + } + } else { + System.out.println("color file " + colors + " not found"); + } + } + public double[] getElevationPos(ElevationRaster srtm, int ilon, int ilat) { double lon = ilon / 1000000. - 180.; double lat = ilat / 1000000. - 90.; @@ -255,12 +323,13 @@ private ElevationRaster srtmForNode(int ilon, int ilat) throws Exception { public static void main(String[] args) throws Exception { if (args.length < 6) { - System.out.println("usage: java CreateLidarImage [type]"); + System.out.println("usage: java CreateLidarImage [type] [color_file]"); System.out.println("\nwhere: type = [bef|hgt] downscale = [1|2|4|..]"); return; } - String format = args.length == 8 ? args[7] : "bef"; + String format = args.length >= 8 ? args[7] : "bef"; + String colors = args.length == 9 ? args[8] : null; 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); + Integer.parseInt(args[4]), Integer.parseInt(args[5]), Integer.parseInt(args[6]), format, colors); } }