removed unused files
This commit is contained in:
afischerdev
parent
36d692da84
commit
fbad694746
6 changed files with 0 additions and 1301 deletions
|
|
@ -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 <srtm-filename> <hgt-data-dir> <srtm-output-dir> [arc seconds (1 or 3,default=3)] [hgt-fallback-data-dir]");
|
||||
System.out.println("or java all <hgt-data-dir> <srtm-output-dir> [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<String> foundList = new ArrayList<>();
|
||||
List<String> 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;
|
||||
}
|
||||
|
||||
}
|
||||
|
|
@ -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);
|
||||
}
|
||||
|
||||
}
|
||||
|
|
@ -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();
|
||||
}
|
||||
|
||||
|
||||
}
|
||||
|
|
@ -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;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -1,166 +0,0 @@
|
|||
package btools.mapcreator;
|
||||
|
||||
/**
|
||||
* This is a wrapper for a 5*5 degree srtm file in ascii/zip-format
|
||||
* <p>
|
||||
* - 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]);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -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;
|
||||
}
|
||||
}
|
||||
Loading…
Reference in a new issue