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]"); } } 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; } }