Fabio/brouter
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions 1

Merge pull request #553 from zod/elev-source-cleanup

Browse source 
Cleanup hgt reader
This commit is contained in:
afischerdev 2023-05-17 10:48:35 +02:00 committed by GitHub
commit 11a9843f41
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
2 changed files with 47 additions and 336 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

102
brouter-map-creator/src/main/java/btools/mapcreator/ConvertLidarTile.java
View file

@ -5,22 +5,26 @@ import java.io.BufferedOutputStream;
import java.io.DataInputStream;
import java.io.File;
import java.io.FileInputStream;
import java.io.FileNotFoundException;
import java.io.FileOutputStream;
import java.io.InputStream;
import java.io.OutputStream;
import java.util.Arrays;
import java.util.zip.ZipEntry;
import java.util.zip.ZipInputStream;
public class ConvertLidarTile {
public static int NROWS;
public static int NCOLS;
public static int ROW_LENGTH;
private static int NROWS;
private static int NCOLS;
public static final short NODATA2 = -32767; // hgt-formats nodata
public static final short NODATA = Short.MIN_VALUE;
public static final int SRTM3_ROW_LENGTH = 1200; // number of elevation values per line
public static final int SRTM1_ROW_LENGTH = 3600; //-- New file resolution is 3601x3601
private static final String HGT_FILE_EXT = ".hgt";
private static final int HGT_BORDER_OVERLAP = 1;
private static final int HGT_3ASEC_ROWS = 1201; // 3 arc second resolution (90m)
private static final int HGT_3ASEC_FILE_SIZE = HGT_3ASEC_ROWS * HGT_3ASEC_ROWS * Short.BYTES;
private static final int HGT_1ASEC_ROWS = 3601; // 1 arc second resolution (30m)
static short[] imagePixels;
@ -30,8 +34,8 @@ public class ConvertLidarTile {
for (; ; ) {
ZipEntry ze = zis.getNextEntry();
if (ze == null) break;
if (ze.getName().endsWith(".hgt")) {
readHgtFromStream(zis, rowOffset, colOffset, 1200);
if (ze.getName().toLowerCase().endsWith(HGT_FILE_EXT)) {
readHgtFromStream(zis, rowOffset, colOffset, HGT_3ASEC_ROWS);
return;
}
}
@ -40,20 +44,20 @@ public class ConvertLidarTile {
}
}
private static void readHgtFromStream(InputStream is, int rowOffset, int colOffset, int row_length)
private static void readHgtFromStream(InputStream is, int rowOffset, int colOffset, int rowLength)
throws Exception {
DataInputStream dis = new DataInputStream(new BufferedInputStream(is));
for (int ir = 0; ir < row_length + 1; ir++) {
for (int ir = 0; ir < rowLength; ir++) {
int row = rowOffset + ir;
for (int ic = 0; ic < row_length + 1; ic++) {
for (int ic = 0; ic < rowLength; ic++) {
int col = colOffset + ic;
int i1 = dis.read(); // msb first!
int i0 = dis.read();
if (i0 == -1 || i1 == -1)
throw new RuntimeException("unexcepted end of file reading hgt entry!");
throw new RuntimeException("unexpected end of file reading hgt entry!");
short val = (short) ((i1 << 8) | i0);
@ -191,80 +195,56 @@ public class ConvertLidarTile {
}
public SrtmRaster getRaster(File f, double lon, double lat) throws Exception {
long fileSize;
InputStream inputStream;
if (f.getName().toLowerCase().endsWith(".zip")) {
ZipInputStream zis = new ZipInputStream(new BufferedInputStream(new FileInputStream(f)));
try {
for (; ; ) {
ZipEntry ze = zis.getNextEntry();
if (ze == null) break;
if (ze.getName().toLowerCase().endsWith(".hgt")) {
if (ze.getSize() > ((SRTM3_ROW_LENGTH + 1) * (SRTM3_ROW_LENGTH + 1) * 2)) {
ROW_LENGTH = SRTM1_ROW_LENGTH;
} else {
ROW_LENGTH = SRTM3_ROW_LENGTH;
if (ze == null) {
throw new FileNotFoundException(f.getName() + " doesn't contain a " + HGT_FILE_EXT + " file.");
}
// stay a 1x1 raster
NROWS = ROW_LENGTH + 1;
NCOLS = ROW_LENGTH + 1;
imagePixels = new short[NROWS * NCOLS]; // 650 MB !
// prefill as NODATA
for (int row = 0; row < NROWS; row++) {
for (int col = 0; col < NCOLS; col++) {
imagePixels[row * NCOLS + col] = NODATA;
}
}
readHgtFromStream(zis, 0, 0, ROW_LENGTH);
if (ze.getName().toLowerCase().endsWith(HGT_FILE_EXT)) {
fileSize = ze.getSize();
inputStream = zis;
break;
}
}
} finally {
zis.close();
}
} else {
if (f.length() > ((SRTM3_ROW_LENGTH + 1) * (SRTM3_ROW_LENGTH + 1) * 2)) {
ROW_LENGTH = SRTM1_ROW_LENGTH;
} else {
ROW_LENGTH = SRTM3_ROW_LENGTH;
fileSize = f.length();
inputStream = new FileInputStream(f);
}
// stay a 1x1 raster
NROWS = ROW_LENGTH + 1;
NCOLS = ROW_LENGTH + 1;
imagePixels = new short[NROWS * NCOLS]; // 650 MB !
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
for (int row = 0; row < NROWS; row++) {
for (int col = 0; col < NCOLS; col++) {
imagePixels[row * NCOLS + col] = NODATA;
}
}
FileInputStream dis = new FileInputStream(f);
try {
readHgtFromStream(dis, 0, 0, ROW_LENGTH);
} finally {
dis.close();
}
}
System.out.println("read file " + f + " rl " + ROW_LENGTH);
boolean halfCol5 = false; // no halfcol tiles in lidar data (?)
Arrays.fill(imagePixels, NODATA);
readHgtFromStream(inputStream, 0, 0, rowLength);
inputStream.close();
SrtmRaster raster = new SrtmRaster();
raster.nrows = NROWS;
raster.ncols = NCOLS;
raster.halfcol = halfCol5;
raster.halfcol = false; // assume full resolution
raster.noDataValue = NODATA;
raster.cellsize = 1. / (double) ROW_LENGTH;
raster.cellsize = 1. / (double) (rowLength - HGT_BORDER_OVERLAP);
raster.xllcorner = (int) (lon < 0 ? lon - 1 : lon); //onDegreeStart - raster.cellsize;
raster.yllcorner = (int) (lat < 0 ? lat - 1 : lat); //latDegreeStart - raster.cellsize;
raster.eval_array = imagePixels;
return raster;
}
}

269
brouter-map-creator/src/main/java/btools/mapcreator/HgtReader.java
View file

@ -1,269 +0,0 @@
// License: GPL. For details, see LICENSE file.
package btools.mapcreator;
import java.io.BufferedInputStream;
import java.io.File;
import java.io.FileInputStream;
import java.io.FileNotFoundException;
import java.io.InputStream;
import java.nio.ByteBuffer;
import java.nio.ByteOrder;
import java.nio.ShortBuffer;
import java.nio.channels.FileChannel;
import java.util.HashMap;
import java.util.Map;
import java.util.zip.ZipEntry;
import java.util.zip.ZipInputStream;
/**
* adapted from https://github.com/JOSM/josm-plugins/blob/master/ElevationProfile/src/org/openstreetmap/josm/plugins/elevation/HgtReader.java
* <p>
* Class HgtReader reads data from SRTM HGT files. Currently this class is restricted to a resolution of 3 arc seconds.
* <p>
* SRTM data files are available at the <a href="http://dds.cr.usgs.gov/srtm/version2_1/SRTM3">NASA SRTM site</a>
*
* @author Oliver Wieland &lt;oliver.wieland@online.de&gt;
*/
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<String, ShortBuffer> cache = new HashMap<>();
public HgtReader(String folder) {
srtmFolder = folder;
}
public static double getElevationFromHgt(double lat, double lon) {
try {
String file = getHgtFileName(lat, lon);
if (DEBUG) System.out.println("SRTM buffer " + file + " for " + lat + " " + lon);
// given area in cache?
if (!cache.containsKey(file)) {
// fill initial cache value. If no file is found, then
// we use it as a marker to indicate 'file has been searched
// but is not there'
cache.put(file, null);
// Try all resource directories
//for (String location : Main.pref.getAllPossiblePreferenceDirs())
{
String fullPath = new File(srtmFolder, file + HGT_EXT).getPath();
File f = new File(fullPath);
if (f.exists()) {
// found something: read HGT file...
ShortBuffer data = readHgtFile(fullPath);
// ... and store result in cache
cache.put(file, data);
//break;
} else {
fullPath = new File(srtmFolder, file + ZIP_EXT).getPath();
f = new File(fullPath);
if (f.exists()) {
ZipInputStream zis = new ZipInputStream(new BufferedInputStream(new FileInputStream(f)));
try {
for (; ; ) {
ZipEntry ze = zis.getNextEntry();
if (ze == null) break;
if (ze.getName().toLowerCase().endsWith(HGT_EXT)) {
// System.out.println("read zip " + ze.getName());
ShortBuffer data = readHgtStream(zis);
// ... and store result in cache
cache.put(file, data);
break;
}
zis.closeEntry();
}
} finally {
zis.close();
}
}
}
}
}
// read elevation value
return readElevation(lat, lon);
} catch (FileNotFoundException e) {
System.err.println("SRTM Get elevation from HGT " + lat + ", " + lon + " failed: => " + e.getMessage());
// no problem... file not there
return NO_ELEVATION;
} catch (Exception ioe) {
// oops...
ioe.printStackTrace(System.err);
// fallback
return NO_ELEVATION;
}
}
@SuppressWarnings("resource")
private static ShortBuffer readHgtFile(String file) throws Exception {
if (file == null) throw new Exception("no srmt file " + file);
FileChannel fc = null;
ShortBuffer sb = null;
try {
// Eclipse complains here about resource leak on 'fc' - even with 'finally' clause???
fc = new FileInputStream(file).getChannel();
// choose the right endianness
ByteBuffer bb = ByteBuffer.allocateDirect((int) fc.size());
while (bb.remaining() > 0) fc.read(bb);
bb.flip();
//sb = bb.order(ByteOrder.LITTLE_ENDIAN).asShortBuffer();
sb = bb.order(ByteOrder.BIG_ENDIAN).asShortBuffer();
} finally {
if (fc != null) fc.close();
}
return sb;
}
// @SuppressWarnings("resource")
private static ShortBuffer readHgtStream(InputStream zis) throws Exception {
if (zis == null) throw new Exception("no srmt stream ");
ShortBuffer sb = null;
try {
// choose the right endianness
byte[] bytes = zis.readAllBytes();
ByteBuffer bb = ByteBuffer.allocate(bytes.length);
bb.put(bytes, 0, bytes.length);
//while (bb.remaining() > 0) zis.read(bb, 0, size);
//ByteBuffer bb = ByteBuffer.allocate(zis.available());
//Channels.newChannel(zis).read(bb);
bb.flip();
//sb = bb.order(ByteOrder.LITTLE_ENDIAN).asShortBuffer();
sb = bb.order(ByteOrder.BIG_ENDIAN).asShortBuffer();
} finally {
}
return sb;
}
/**
* Reads the elevation value for the given coordinate.
* <p>
* See also <a href="http://gis.stackexchange.com/questions/43743/how-to-extract-elevation-from-hgt-file">stackexchange.com</a>
*
* @param lat, lon the coordinate to get the elevation data for
* @return the elevation value or <code>Double.NaN</code>, if no value is present
*/
public static double readElevation(double lat, double lon) {
String tag = getHgtFileName(lat, lon);
ShortBuffer sb = cache.get(tag);
if (sb == null) {
return NO_ELEVATION;
}
if (DEBUG) System.out.println("SRTM buffer size " + sb.capacity() + " limit " + sb.limit());
try {
int rowLength = HGT3_ROW_LENGTH;
int resolution = HGT3_RES;
if (sb.capacity() > (HGT3_ROW_LENGTH * HGT3_ROW_LENGTH)) {
rowLength = HGT1_ROW_LENGTH;
resolution = HGT1_RES;
}
// see http://gis.stackexchange.com/questions/43743/how-to-extract-elevation-from-hgt-file
double fLat = frac(lat) * SECONDS_PER_MINUTE;
double fLon = frac(lon) * SECONDS_PER_MINUTE;
// compute offset within HGT file
int row = (int) Math.round((fLat) * SECONDS_PER_MINUTE / resolution);
int col = (int) Math.round((fLon) * SECONDS_PER_MINUTE / resolution);
if (lon < 0) col = rowLength - col - 1;
if (lat > 0) row = rowLength - row - 1;
//row = rowLength - row;
int cell = (rowLength * (row)) + col;
//int cell = ((rowLength * (latitude)) + longitude);
if (DEBUG)
System.out.println("Read SRTM elevation data from row/col/cell " + row + "," + col + ", " + cell + ", " + sb.limit());
// valid position in buffer?
if (cell < sb.limit()) {
short ele = sb.get(cell);
// check for data voids
if (ele == HGT_VOID) {
return NO_ELEVATION;
} else {
return ele;
}
} else {
return NO_ELEVATION;
}
} catch (Exception e) {
System.err.println("error at " + lon + " " + lat + " ");
e.printStackTrace();
}
return NO_ELEVATION;
}
/**
* Gets the associated HGT file name for the given way point. Usually the
* format is <tt>[N|S]nn[W|E]mmm.hgt</tt> where <i>nn</i> is the integral latitude
* without decimals and <i>mmm</i> 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();
}
}
}