264 lines
7.5 KiB
Java
264 lines
7.5 KiB
Java
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;
|
|
}
|
|
}
|