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 0000000..9b251fc --- /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; + } + + +}