- mkgmap-dev - The mkgmap lists

[PATCH] Creating contour lines from DEM data
by Christian Gawron 05 Jul '09

05 Jul '09

Dear all, the attached java classes together with the patch enables mkgmap to create contour lines directly from digital elevation model (DEM) data. This may be useful for those who don't want to ceate contour lines with other tools and store them in huge(!) files. It introduces the following new options: --contours If set, create contours. --dem-type Valid values are CGIAR (CGIAR geotiff), ASTER (ASTER geotiff) or SRTM (.hgt files). Default is CGIAR. --dem-path Path to the DEM data files. The default is type-dependant. --dem-increment Verical distance between the contour lines (default is 10m). The algorithm used is a modified version of that described in "An Adaptive Grid Contouring Algorithm" by Downing and Zoraster (see http://mapcontext.com/autocarto/proceedings/auto-carto-5/pdf/an-adaptive-gr…) and uses bicubic interpolation of the DEM. Please feel free to send me improvements, bug reports etc. Best wishes Christian Index: src/uk/me/parabola/mkgmap/Option.java =================================================================== --- src/uk/me/parabola/mkgmap/Option.java (Revision 1077) +++ src/uk/me/parabola/mkgmap/Option.java (Arbeitskopie) @@ -34,7 +34,7 @@ } } - protected Option(String option, String value) { + public Option(String option, String value) { this.option = option; this.value = value; } Index: src/uk/me/parabola/mkgmap/CommandArgsReader.java =================================================================== --- src/uk/me/parabola/mkgmap/CommandArgsReader.java (Revision 1077) +++ src/uk/me/parabola/mkgmap/CommandArgsReader.java (Arbeitskopie) @@ -129,7 +129,7 @@ * @param option The option name. * @param value Its value. */ - private void addOption(String option, String value) { + public void addOption(String option, String value) { CommandOption opt = new CommandOption(option, value); addOption(opt); } @@ -181,13 +181,21 @@ return arglist.iterator(); } + public ArgList getArgList() { + return arglist; + } + + public EnhancedProperties getArgs() { + return args; + } + /** * Read a config file that contains more options. When the number of * options becomes large it is more convenient to place them in a file. * * @param filename The filename to obtain options from. */ - private void readConfigFile(String filename) { + public void readConfigFile(String filename) { Options opts = new Options(new OptionProcessor() { public void processOption(Option opt) { log.debug("incoming opt", opt.getOption(), opt.getValue()); @@ -206,14 +214,14 @@ * the argument to be processed in order. Options can be intersperced with * filenames. The options take effect where they appear. */ - interface ArgType { + public interface ArgType { public abstract void processArg(); } /** * A filename. */ - class Filename implements ArgType { + public class Filename implements ArgType { private final String name; private boolean useFilenameAsMapname = true; @@ -270,7 +278,7 @@ /** * An option argument. A key value pair. */ - class CommandOption implements ArgType { + public class CommandOption implements ArgType { private final Option option; private CommandOption(Option option) { @@ -297,7 +305,7 @@ /** * The arguments are held in this list. */ - class ArgList implements Iterable<CommandArgsReader.ArgType> { + public class ArgList implements Iterable<CommandArgsReader.ArgType> { private final List<ArgType> alist; private int filenameCount; Index: src/uk/me/parabola/mkgmap/reader/MapperBasedMapDataSource.java =================================================================== --- src/uk/me/parabola/mkgmap/reader/MapperBasedMapDataSource.java (Revision 1077) +++ src/uk/me/parabola/mkgmap/reader/MapperBasedMapDataSource.java (Arbeitskopie) @@ -33,6 +33,7 @@ import uk.me.parabola.mkgmap.general.MapPoint; import uk.me.parabola.mkgmap.general.MapShape; import uk.me.parabola.mkgmap.general.RoadNetwork; +import uk.me.parabola.mkgmap.reader.dem.DEM; import uk.me.parabola.util.Configurable; import uk.me.parabola.util.EnhancedProperties; @@ -146,5 +147,8 @@ // the overview section. mapper.addShape(background); } + if (getConfig().getProperty("contours", false)) { + DEM.createContours(mapper, getConfig()); + } } } /* * Copyright (C) 2009 Christian Gawron * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License version 2 as * published by the Free Software Foundation. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * * Author: Christian Gawron * Create date: 03-Jul-2009 */ package uk.me.parabola.mkgmap.reader.dem; import java.io.*; import java.nio.channels.FileChannel; import java.nio.MappedByteBuffer; import java.util.Vector; import java.util.Locale; import com.sun.media.jai.codec.*; import javax.media.jai.*; import java.awt.image.renderable.ParameterBlock; import java.awt.image.Raster; import java.awt.image.DataBuffer; import uk.me.parabola.imgfmt.Utils; import uk.me.parabola.imgfmt.app.Area; import uk.me.parabola.imgfmt.app.Coord; import uk.me.parabola.mkgmap.general.MapDetails; import uk.me.parabola.mkgmap.general.MapPoint; import uk.me.parabola.mkgmap.general.MapLine; import uk.me.parabola.util.EnhancedProperties; import static java.nio.channels.FileChannel.MapMode.READ_ONLY; /** * Create contour lines using an algorithm similar to that described in * <a href="http://mapcontext.com/autocarto/proceedings/auto-carto-5/pdf/an-adaptive-gr…">An Adaptive Grid Contouring Algorithm</a> by Downing and Zoraster. */ public abstract class DEM { final static double epsilon = 1e-9; final static double delta = 1.5; final static int maxPoints=200000; final static double minDist = 15; final static double maxDist = 21; final static double step = 0.01; final static double semiMajorAxis = 6378137.0; final static double inverseFlattening = 298.257223563; static int M=1200; static int N=M; static double res = 1.0/N; short values[] = null; int lat; int lon; abstract double ele(int x, int y); abstract void read(int minLon, int minLat, int maxLon, int maxLat); abstract void serializeCopyRight(Writer out) throws IOException; public static void createContours(MapDetails mapper, EnhancedProperties config) { Area bounds = mapper.getBounds(); double minLat = Utils.toDegrees(bounds.getMinLat()); double minLon = Utils.toDegrees(bounds.getMinLong()); double maxLat = Utils.toDegrees(bounds.getMaxLat()); double maxLon = Utils.toDegrees(bounds.getMaxLong()); System.out.printf("bounds: %f %f %f %f\n", minLat, minLon, maxLat, maxLon); DEM data; String demType = config.getProperty("dem-type", "CGIAR"); String dataPath; if (demType.equals("ASTER")) { dataPath = config.getProperty("dem-path", "ASTER"); data = new ASTERGeoTiff(dataPath, minLat, minLon, maxLat, maxLon); } else if (demType.equals("CGIAR")) { dataPath = config.getProperty("dem-path", "CGIAR"); data = new CGIARGeoTiff(dataPath, minLat, minLon, maxLat, maxLon); } else { dataPath = config.getProperty("dem-path", "SRTM"); data = new HGTDEM(dataPath, minLat, minLon, maxLat, maxLon); } Isolines lines = data.new Isolines(data, minLat, minLon, maxLat, maxLon); int increment = config.getProperty("dem-increment", 10); lines.addLevels(0, increment); for (Isolines.Isoline line : lines.isolines) { MapLine contour = new MapLine(); if ((line.level % 200) == 0) { contour.setType(0x22); // major contour contour.setMinResolution(18); } else if ((line.level % 50) == 0) { contour.setType(0x21); // major contour contour.setMinResolution(20); } else { contour.setType(0x20); // major contour contour.setMinResolution(22); } contour.setName(String.format("%.2f", line.level * 3.2808399)); contour.setPoints(line.points); mapper.addLine(contour); } } public static class HGTDEM extends DEM { MappedByteBuffer buffer = null; public HGTDEM(String dataPath, double minLat, double minLon, double maxLat, double maxLon) { this.lat = (int) minLat; this.lon = (int) minLon; if (maxLat > lat+1 || maxLon > lon+1) throw new RuntimeException("Area too large (must not span more than one SRTM file)"); String northSouth = lat < 0 ? "S" : "N"; String eastWest = lon > 0 ? "E" : "W"; String fileName = String.format("%s/%s%02d%s%03d.hgt", dataPath, northSouth, lat < 0 ? -lat : lat, eastWest, lon < 0 ? -lon : lon); try { FileInputStream is = new FileInputStream(fileName); buffer = is.getChannel().map(READ_ONLY, 0, 2*(M+1)*(M+1)); } catch (Exception e) { throw new RuntimeException(e); } } public void read(int minLon, int minLat, int maxLon, int maxLat) { } public double ele(int x, int y) { return buffer.getShort(2*((M-y)*(M+1)+x))+delta; } public void serializeCopyRight(Writer out) throws IOException { out.write(" <copyright>\n"); out.write(" Contour lines generated from DEM data by NASA\n"); out.write(" </copyright>\n"); } } public static class CGIARGeoTiff extends DEM { Raster raster; String fileName; int minLat, minLon, maxLat, maxLon; PlanarImage image; public CGIARGeoTiff(String dataPath, double minLat, double minLon, double maxLat, double maxLon) { this.lat = ((int) (minLat/5))*5; this.lon = ((int) (minLon/5))*5; if (maxLat > lat+5 || maxLon > lon+5) throw new RuntimeException("Area too large (must not span more than one CGIAR GeoTIFF)"); int tileX, tileY; tileX = (180 + lon)/5 + 1; tileY = (60 - lat)/5; this.fileName = String.format("%s/srtm_%02d_%02d.tif", dataPath, tileX, tileY); init(); } public void serializeCopyRight(Writer out) throws IOException { out.write(" <copyright>\n"); out.write(" Contour lines generated from improved SRTM data by CIAT-CSI (see http://srtm.csi.cgiar.org)\n"); out.write(" </copyright>\n"); } public void read(int minLon, int minLat, int maxLon, int maxLat) { this.minLon = minLon; this.minLat = minLat; this.maxLon = maxLon; this.maxLat = maxLat; raster = image.getData(new java.awt.Rectangle(minLon, 6000-maxLat-1, maxLon-minLon+1, maxLat-minLat+1)); System.out.printf("read: %d %d %d %d\n", minLon, 6000-maxLat-1, maxLon-minLon+1, maxLat-minLat+1); } void init() { System.out.printf("CGIAR GeoTIFF: %s\n", fileName); N = 6000; M = 6000; res = 5.0/M; try { SeekableStream s = new FileSeekableStream(fileName); ParameterBlock pb = new ParameterBlock(); pb.add(s); TIFFDecodeParam param = new TIFFDecodeParam(); pb.add(param); RenderedOp op = JAI.create("tiff", pb); image = op.createInstance(); System.out.printf("Image: %d %d %d %d\n", image.getWidth(), image.getHeight(), image.getNumXTiles(), image.getNumYTiles()); } catch (Exception e) { throw new RuntimeException(e); } } public double ele(int x, int y) { try { int elevation = raster.getPixel(x, 6000-y-1, (int[])null)[0]; return elevation+delta; } catch (ArrayIndexOutOfBoundsException ex) { System.out.printf("ele: (%d, %d) (%d, %d, %d, %d) %s\n", x, 6000-y-1, raster.getMinX(), raster.getMinY(), raster.getWidth(), raster.getHeight(), ex.toString()); throw ex; } } } public static class ASTERGeoTiff extends DEM { Raster raster; String fileName; int minLat, minLon, maxLat, maxLon; PlanarImage image; public ASTERGeoTiff(String dataPath, double minLat, double minLon, double maxLat, double maxLon) { this.lat = (int) minLat; this.lon = (int) minLon; if (maxLat > lat+1 || maxLon > lon+1) throw new RuntimeException("Area too large (must not span more than one ASTER GeoTIFF)"); String northSouth = lat < 0 ? "S" : "N"; String eastWest = lon > 0 ? "E" : "W"; fileName = String.format("%s/ASTGTM_%s%02d%s%03d_dem.tif", dataPath, northSouth, lat < 0 ? -lat : lat, eastWest, lon < 0 ? -lon : lon); init(); } public void serializeCopyRight(Writer out) throws IOException { out.write(" <copyright>\n"); out.write(" Contour lines generated from DGM data by ASTER (see https://wist.echo.nasa.gov/~wist/api/imswelcome)\n"); out.write(" </copyright>\n"); } public void read(int minLon, int minLat, int maxLon, int maxLat) { this.minLon = minLon; this.minLat = minLat; this.maxLon = maxLon; this.maxLat = maxLat; raster = image.getData(new java.awt.Rectangle(minLon, 3601-maxLat-1, maxLon-minLon+1, maxLat-minLat+1)); System.out.printf("read: %d %d %d %d\n", minLon, 3601-maxLat-1, maxLon-minLon+1, maxLat-minLat+1); } void init() { System.out.printf("ASTER GeoTIFF: %s\n", fileName); N = 3600; M = 3600; res = 1.0/M; try { SeekableStream s = new FileSeekableStream(fileName); ParameterBlock pb = new ParameterBlock(); pb.add(s); TIFFDecodeParam param = new TIFFDecodeParam(); pb.add(param); RenderedOp op = JAI.create("tiff", pb); image = op.createInstance(); System.out.printf("Image: %d %d %d %d\n", image.getWidth(), image.getHeight(), image.getNumXTiles(), image.getNumYTiles()); } catch (Exception e) { throw new RuntimeException(e); } } public double ele(int x, int y) { try { int elevation = raster.getPixel(x, 3601-y-1, (int[])null)[0]; return elevation+delta; } catch (ArrayIndexOutOfBoundsException ex) { System.out.printf("ele: (%d, %d) (%d, %d, %d, %d) %s\n", x, 3601-y-1, raster.getMinX(), raster.getMinY(), raster.getWidth(), raster.getHeight(), ex.toString()); throw ex; } } } private static void debug(String format, Object ... args) { //System.out.printf(format + "\n", args); } int lastXi = -1; int lastYi = -1; final static int bcInv[][]= { { 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, { 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0}, { -3, 0, 0, 3, 0, 0, 0, 0, -2, 0, 0, -1, 0, 0, 0, 0}, { 2, 0, 0, -2, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0}, { 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0}, { 0, 0, 0, 0, -3, 0, 0, 3, 0, 0, 0, 0, -2, 0, 0, -1}, { 0, 0, 0, 0, 2, 0, 0, -2, 0, 0, 0, 0, 1, 0, 0, 1}, { -3, 3, 0, 0, -2, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, { 0, 0, 0, 0, 0, 0, 0, 0, -3, 3, 0, 0, -2, -1, 0, 0}, { 9, -9, 9, -9, 6, 3, -3, -6, 6, -6, -3, 3, 4, 2, 1, 2}, { -6, 6, -6, 6, -4, -2, 2, 4, -3, 3, 3, -3, -2, -1, -1, -2}, { 2, -2, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, { 0, 0, 0, 0, 0, 0, 0, 0, 2, -2, 0, 0, 1, 1, 0, 0}, { -6, 6, -6, 6, -3, -3, 3, 3, -4, 4, 2, -2, -2, -2, -1, -1}, { 4, -4, 4, -4, 2, 2, -2, -2, 2, -2, -2, 2, 1, 1, 1, 1} }; static int lastId=1000000000; static double lastX = 0; static double lastY = 0; static int edge[] = new int[2]; static int off0[][] = {{ 0, 0}, { 0, 0}, { 0, 1}, { 1, 1}}; static int off1[][] = {{ 1, 0}, { 0, 1}, { 1, 1}, { 1, 0}}; static int brd[] = { 1, 2, 4, 8 }; static int inv[] = { 4, 8, 1, 2 }; static int rev[] = { 2, 3, 0, 1 }; static int mov[][] = {{ 0, -1}, {-1, 0}, { 0, 1}, { 1, 0}}; double bc[][] = new double[4][4]; double bc_y[] = new double[4]; double bc_y1[] = new double[4]; double bc_y2[] = new double[4]; double bc_y12[] = new double[4]; double bc_Coeff[] = new double[16]; double bc_x[] = new double[16]; private void recalculateCoefficients(int xi, int yi) { double v00, vp0, v0p, vpp; double vm0, v0m, vpm, vmp, vmm, vmP, vPm; double vP0, v0P, vPp, vpP, vPP; v00 = ele(xi, yi); v0p = ele(xi, yi+1); vpp = ele(xi+1, yi+1); vp0 = ele(xi+1, yi); vm0 = ele(xi-1, yi); v0m = ele(xi, yi-1); vmp = ele(xi-1, yi+1); vpm = ele(xi+1, yi-1); vmm = ele(xi-1, yi-1); vmP = ele(xi+2, yi-1); vPm = ele(xi-1, yi+2); vP0 = ele(xi+2, yi); v0P = ele(xi, yi+2); vPp = ele(xi+2, yi+1); vpP = ele(xi+1, yi+2); vPP = ele(xi+2, yi+2); bc_y[0] = v00; bc_y[1] = vp0; bc_y[2] = vpp; bc_y[3] = v0p; bc_y1[0] = (vp0-vm0)/2; bc_y1[1] = (vP0-v00)/2; bc_y1[2] = (vPp-v0p)/2; bc_y1[3] = (vpp-vmp)/2; bc_y2[0] = (v0p-v0m)/2; bc_y2[1] = (vpp-vpm)/2; bc_y2[2] = (vpP-vp0)/2; bc_y2[3] = (v0P-v00)/2; bc_y12[0] = (vpp - vpm - vmp + vmm) / 4; bc_y12[0] = (vPp - vPm - v0p + v0m) / 4; bc_y12[2] = (vPP - vP0 - v0P + v00) / 4; bc_y12[0] = (vpP - vp0 - vmP + vm0) / 4; int j, i; double s; for (i=0; i<4; i++) { bc_x[i]=bc_y[i]; bc_x[i+4]=bc_y1[i]; bc_x[i+8]=bc_y2[i]; bc_x[i+12]=bc_y12[i]; } for (i=0; i<16; i++) { s = 0; for (int k=0; k<16; k++) s += bcInv[i][k]*bc_x[k]; bc_Coeff[i] = s; } int l = 0; for (i=0; i<4; i++) for (j=0; j<4; j++) bc[i][j] = bc_Coeff[l++]; } public double gradient(double lat, double lon, double[] grad) { grad[0] = 0; grad[1] = 0; double x = (lon-this.lon)/res; double y = (lat-this.lat)/res; int xi = (int) x; int yi = (int) y; if (lastXi != xi || lastYi != yi) { debug("new Cell for interpolation: %d %d", xi, yi); recalculateCoefficients(xi, yi); lastXi = xi; lastYi = yi; } double t = x - xi; double u = y - yi; if (xi < 0 || xi > N+1 || yi < 0 || yi > N+1) throw new IndexOutOfBoundsException(String.format("(%f, %f)->(%d, %d)", lat, lon, xi, yi)); double val = 0; for (int i=3; i>=0; i--) { val = t*val + ((bc[i][3]*u + bc[i][2])*u + bc[i][1])*u + bc[i][0]; grad[0] = u*grad[0] + (3*bc[3][i]*t + 2*bc[2][i])*t + bc[1][i]; grad[1] = t*grad[1] + (3*bc[i][3]*t + 2*bc[i][2])*t + bc[i][1]; } return val; } public double elevation(double lat, double lon) { double x = (lon-this.lon)/res; double y = (lat-this.lat)/res; int xi = (int) x; int yi = (int) y; if (lastXi != xi || lastYi != yi) { debug("new Cell for interpolation: %d %d", xi, yi); recalculateCoefficients(xi, yi); lastXi = xi; lastYi = yi; } double t = x - xi; double u = y - yi; if (xi < 0 || xi > N+1 || yi < 0 || yi > N+1) throw new IndexOutOfBoundsException(String.format("(%f, %f)->(%d, %d)", lat, lon, xi, yi)); double val = 0; for (int i=3; i>=0; i--) { val = t*val + ((bc[i][3]*u + bc[i][2])*u + bc[i][1])*u + bc[i][0]; } return val; } public double elevation(int x, int y) { if (x < 0 || x > N || y < 0 || y > N) throw new IndexOutOfBoundsException(String.format("elevation: %d %d", x, y)); return ele(x, y); } class Isolines { DEM data; int minX; int maxX; int minY; int maxY; double min; double max; Vector<Isoline> isolines = new Vector<Isoline>(); class Isoline { int id; Vector<Coord> points; double level; boolean isClosed; private Isoline(double level) { this.level = level; isClosed = false; id = lastId++; points = new Vector<Coord>(); } private class Edge implements Brent.Function { double x0, y0, x1, y1; Edge(double x0, double y0, double x1, double y1) { this.x0 = x0; this.y0 = y0; this.x1 = x1; this.y1 = y1; } public double eval(double d) { double f = data.elevation(x0 + d * (x1-x0), y0 + d * (y1-y0)) - level; //System.out.printf("evalEdge: %f %f\n", d, f); return f; } } private class FN implements Brent.Function { double x0, y0; double dx, dy; public void setParameter(double x0, double y0, double dx, double dy) { this.x0 = x0; this.y0 = y0; this.dx = dx; this.dy = dy; } public double eval(double t) { double f = data.elevation(y0+t*dy, x0+t*dx) - level; return f; } } private FN fn = new FN(); double grad[] = new double[2]; double px[] = new double[4]; double py[] = new double[4]; int edges[] = new int[4]; boolean addCell(Position p, int direction) { debug("addCell: %f %d %d %d %d", level, p.ix, p.iy, p.edge, direction); int c = 0; for (int k=0; k<4; k++) { if (k == p.edge) continue; int x0 = p.ix + off0[k][0]; int y0 = p.iy + off0[k][1]; int x1 = p.ix + off1[k][0]; int y1 = p.iy + off1[k][1]; double l0 = elevation(x0, y0) - level; double l1 = elevation(x1, y1) - level; if (Math.abs(l1) < epsilon || l0*l1 < 0) { edges[c] = k; Brent.Function f = new Edge(data.lat + y0 * DEM.res, data.lon + x0 * DEM.res, data.lat + y1 * DEM.res, data.lon + x1 * DEM.res); double f0 = elevation(x0, y0) - level; double f1 = elevation(x1, y1) - level; double delta; if (Math.abs(1) < epsilon) { delta = 1; } else if (Math.abs(f0) < epsilon) throw new RuntimeException("implementation error!"); else delta = Brent.zero(f, epsilon, 1-epsilon); px[c] = data.lon + (x0+delta*(x1-x0))*DEM.res; py[c] = data.lat + (y0+delta*(y1-y0))*DEM.res; c++; } } if (c == 1) { p.edge = edges[0]; double px0 = p.x; double py0 = p.y; p.x = px[0]; p.y = py[0]; double px1 = p.x; double py1 = p.y; double xMin = data.lon + p.ix * DEM.res; double xMax = xMin + DEM.res; double yMin = data.lat + p.iy * DEM.res; double yMax = yMin + DEM.res; refineAdaptively(xMin, yMin, xMax, yMax, px0, py0, px1, py1, direction, maxDist); addPoint(p.x, p.y, direction); p.moveCell(); return true; } else { debug("addCellByStepping: %d", c); return addCellByStepping(p, direction, c, edges, px, py); } } private void refineAdaptively(double xMin, double yMin, double xMax, double yMax, double x0, double y0, double x1, double y1, int direction, double maxDist) { double dist = orthodromicDistance(x0, y0, x1, y1); if (dist > maxDist) { double dx = x1-x0; double dy = y1-y0; double xm, ym, f0, f1, t0, t1, t; Brent.Function f; xm = x0 + 0.5*dx; ym = y0 + 0.5*dy; double n = Math.sqrt(dx*dx+dy*dy); fn.setParameter(xm, ym, -dy/n, dx/n); f = fn; t0 = -0.05*res; t1 = 0.05*res; f0 = f.eval(t0); f1 = f.eval(t1); int count = 0; while (f0 * f1 > 0 && count++ < 20) { if ((count & 1) > 0) t0 -= 0.05*res; else t1 += 0.05*res; f0 = f.eval(t0); f1 = f.eval(t1); debug("refine: %f %f %f %f", t0, t1, f0, f1); } if (f0 * f1 < 0) { t = Brent.zero(f, t0, t1); xm -= t*dy; ym += t*dx; } else { debug("refine failed: %f %f %f %f", t0, t1, f0, f1); if (false) throw new RuntimeException(String.format("refine failed: %f %f %f %f %f %f %f %f", xMin, yMin, xMax, yMax, x0, y0, x1, y1)); return; } if (xm > xMin && xm < xMax && ym > yMin && ym < yMax) refineAdaptively(xMin, yMin, xMax, yMax, x0, y0, xm, ym, direction, maxDist*1.1); addPoint(xm, ym, direction); if (xm > xMin && xm < xMax && ym > yMin && ym < yMax) refineAdaptively(xMin, yMin, xMax, yMax, xm, ym, x1, y1, direction, maxDist*1.1); } } boolean addCellByStepping(Position p, int direction, int numEdges, int[] edges, double[] px, double[] py) { debug("addCellByStepping: %f %d %d %d %d", level, p.ix, p.iy, p.edge, direction); double xMin = data.lon + p.ix * DEM.res; double xMax = xMin + DEM.res; double yMin = data.lat + p.iy * DEM.res; double yMax = yMin + DEM.res; double dt, t0, t1; double f0, f1; boolean edgeHit = false; double h[] = new double[4]; int dir; double n2 = Math.sqrt(1.0/(grad[0]*grad[0] + grad[1]*grad[1])); double dx; double dy; int count=0; int iMin = -1; double md = 5000; for (int i=0; i<numEdges; i++) { gradient(p.y, p.x, grad); double dist = orthodromicDistance(p.x, p.y, px[i], py[i]); //double dist = (px[i]-p.x)*(px[i]-p.x)*grad[0]*grad[0] + (py[i]-p.y)*(py[i]-p.y)*grad[1]*grad[1]; debug("distance %d: %f", i, dist); if (dist < md && (visited[p.iy*(N+1)+p.ix] & brd[edges[i]]) == 0) { md = dist; iMin = i; } } p.edge = edges[iMin]; double px0 = p.x; double py0 = p.y; p.x = px[iMin]; p.y = py[iMin]; double px1 = p.x; double py1 = p.y; xMin = data.lon + p.ix * DEM.res; xMax = xMin + DEM.res; yMin = data.lat + p.iy * DEM.res; yMax = yMin + DEM.res; refineAdaptively(xMin, yMin, xMax, yMax, px0, py0, px1, py1, direction, maxDist); addPoint(p.x, p.y, direction); p.moveCell(); return true; } private void addPoint(double x, double y, int direction) { double dist = orthodromicDistance(x, y, lastX, lastY); debug("addPoint: %f %f %f", x, y, dist); if (dist > minDist) { if (direction > 0) points.add(0, new Coord(y, x)); else points.add(points.size(), new Coord(y, x)); lastX = x; lastY = y; } } private void close() { points.add(points.size(), points.get(0)); isClosed = true; } private void addMove(int x0, int y0, int x1, int y1, int direction) { double x; double y; if (x0 < minX || x0 >= maxX || y0 < minY || y0 >= maxY) return; double l0 = data.elevation(x0, y0); double l1 = data.elevation(x1, y1); debug("addMove: %d %d %d %d %.2f %.2f %d", x0, y0, x1, y1, l0, l1, direction); if ((l0 < level && l1 < level) || (l0 > level && l1 > level)) throw new RuntimeException("implementation error"); if (l0 == level) { x = data.lon + x0 * DEM.res; y = data.lat + y0 * DEM.res; } else { double delta = (l0-level)/(l0-l1); x = data.lon + (x0 + delta * (x1-x0)) * DEM.res; y = data.lat + (y0 + delta * (y1-y0)) * DEM.res; } double dist = orthodromicDistance(x, y, lastX, lastY); debug("levels: %d %d %f, %d %d %f: %f %f %f", x0, y0, l0, x1, y1, l1, x, y, dist); if (dist > 1) { if (direction > 0) points.add(0, new Coord(y, x)); else points.add(points.size(), new Coord(y, x)); lastX = x; lastY = y; } } } public Isolines(DEM data, int minX, int maxX, int minY, int maxY) { this.data = data; this.minX = minX; this.maxX = maxX; this.minY = minY; this.maxY = maxY; init(); } public Isolines(DEM data, double minLat, double minLon, double maxLat, double maxLon) { System.out.printf("init: %f %f %f %f\n", minLat, minLon, maxLat, maxLon); this.data = data; this.minX = (int) ((minLon-data.lon)/data.res); this.minY = (int) ((minLat-data.lat)/data.res); this.maxX = (int) ((maxLon-data.lon)/data.res); this.maxY = (int) ((maxLat-data.lat)/data.res); init(); } private void init() { System.out.printf("init: %d %d %d %d\n", minX, minY, maxX, maxY); data.read(minX-2, minY-2, maxX+2, maxY+2); // we need some overlap for bicubic interpolation max = -1000; min = 10000; for (int i=minX; i<maxX; i++) for (int j=minY; j<maxY; j++) { if (data.elevation(i, j) < min) min = data.elevation(i, j); if (data.elevation(i, j) > max) max = data.elevation(i, j); } debug("min: %f, max: %f\n", min, max); } double getMinHeight() { return min; } double getMaxHeight() { return max; } public void addLevels(int start, int increment) { for (int level=start; level<max; level+=increment) addLevel(level); } private class Edge implements Brent.Function { double x0, y0, x1, y1, level; Edge(double x0, double y0, double x1, double y1, double level) { this.x0 = x0; this.y0 = y0; this.x1 = x1; this.y1 = y1; this.level = level; } public double eval(double d) { double f = data.elevation(x0 + d * (x1-x0), y0 + d * (y1-y0)) - level; //System.out.printf("evalEdge: %f %f\n", d, f); return f; } } private class Position { int ix, iy; double x, y; int edge; Position(int ix, int iy, double x, double y, int edge) { this.ix = ix; this.iy = iy; this.x = x; this.y = y; this.edge = edge; } Position(Position p) { this.ix = p.ix; this.iy = p.iy; this.x = p.x; this.y = p.y; this.edge = p.edge; } void markEdge() { debug("marking edge: %d %d %d %d", ix, iy, edge, brd[edge]); visited[iy*(N+1)+ix] |= brd[edge]; } void moveCell() { markEdge(); ix += mov[edge][0]; iy += mov[edge][1]; edge = rev[edge]; markEdge(); } } byte visited[] = new byte[(N+1)*(N+1)]; public void addLevel(double level) { if (level < min || level > max) return; System.out.printf("addLevel: %f\n", level); java.util.Arrays.fill(visited, (byte) 0); for (int y=minY; y<maxY; y++) { for (int x=minX; x<maxX; x++) { byte v = 0; // Mark the borders of the cell, represented by the four points (i, j), (i+1, j), (i, j+1), (i+1, j+1), // which are intersected by the contour. The values are: // 1: top // 2: left // 4: bottom // 8: right if (data.elevation(x, y) >= level) { if (data.elevation(x+1, y) < level) { v |= 1; } if (data.elevation(x, y+1) < level) { v |= 2; } } else { if (data.elevation(x+1, y) > level) { v |= 1; } if (data.elevation(x, y+1) > level) { v |= 2; } } int k=-1; if ((v&1) > 0 && (visited[y*(N+1)+x]&1) == 0) { k=0; } else if ((v&2) > 0 && (visited[y*(N+1)+x]&2) == 0) { k=1; } if (k>=0) { int x0 = x + off0[k][0]; int y0 = y + off0[k][1]; int x1 = x + off1[k][0]; int y1 = y + off1[k][1]; try { Brent.Function f = new Edge(data.lat + y0 * DEM.res, data.lon + x0 * DEM.res, data.lat + y1 * DEM.res, data.lon + x1 * DEM.res, level); double f0 = elevation(x0, y0) - level; double f1 = elevation(x1, y1) - level; double delta; if (Math.abs(f0) < epsilon) { delta = 0; } else if (Math.abs(f1) < epsilon) continue; else delta = Brent.zero(f, 0, 1-epsilon); Position p = new Position(x, y, data.lon + (x0+delta*(x1-x0))*DEM.res, data.lat + (y0+delta*(y1-y0))*DEM.res, k); p.markEdge(); isolines.add(traceByStepping(level, p)); } catch (RuntimeException ex) { debug("error: %s", ex.toString()); ex.printStackTrace(); return; } } } } } private Isoline traceByStepping(double level, Position p) { debug("traceByStepping: starting contour %f %d %d %f %f %d", level, p.ix, p.iy, p.x, p.y, p.edge); int direction = 1; int n = 0; Position startP = new Position(p); boolean foundEnd = false; Isoline line = new Isoline(level); while (true) { debug("traceByStepping: %f %d %d %f %f %d", level, p.ix, p.iy, p.x, p.y, p.edge); visited[p.iy*(N+1)+p.ix] |= brd[p.edge]; if (n>0 && p.ix == startP.ix && p.iy == startP.iy && orthodromicDistance(p.x, p.y, startP.x, startP.y) < 5) { debug("closed curve!"); line.close(); break; } else if (p.ix < minX || p.iy < minY || p.ix >= maxX || p.iy >= maxY) { if (foundEnd) // did we already reach one end? { debug("second border reached!"); break; } else { debug("border reached!"); foundEnd = true; n = 0; direction *= -1; p = new Position(startP); p.moveCell(); continue; } } n++; if (!line.addCell(p, direction) || line.points.size() > maxPoints) { debug("ending contour"); isolines.add(line); return line; } } return line; } } /** * Returns the orthodromic distance between two geographic coordinates in WGS84 datum. * The orthodromic distance is the shortest distance between two points * on a sphere's surface. The orthodromic path is always on a great circle. * * @param x1 Longitude of first point (in degrees). * @param y1 Latitude of first point (in degrees). * @param x2 Longitude of second point (in degrees). * @param y2 Latitude of second point (in degrees). * @return The orthodromic distance (in meters). * */ public static double orthodromicDistance(double x1, double y1, double x2, double y2) { x1 = Math.toRadians(x1); y1 = Math.toRadians(y1); x2 = Math.toRadians(x2); y2 = Math.toRadians(y2); final int MAX_ITERATIONS = 100; final double EPS = 5.0E-14; final double F = 1 / inverseFlattening; final double R = 1 - F; double tu1 = R * Math.sin(y1) / Math.cos(y1); double tu2 = R * Math.sin(y2) / Math.cos(y2); double cu1 = 1 / Math.sqrt(tu1 * tu1 + 1); double cu2 = 1 / Math.sqrt(tu2 * tu2 + 1); double su1 = cu1 * tu1; double s = cu1 * cu2; double baz = s * tu2; double faz = baz * tu1; double x = x2 - x1; for (int i = 0; i < MAX_ITERATIONS; i++) { final double sx = Math.sin(x); final double cx = Math.cos(x); tu1 = cu2 * sx; tu2 = baz - su1 * cu2 * cx; final double sy = Math.sqrt(tu1*tu1 + tu2*tu2);; final double cy = s * cx + faz; final double y = Math.atan2(sy, cy); final double SA = s * sx / sy; final double c2a = 1 - SA * SA; double cz = faz + faz; if (c2a > 0) { cz = -cz / c2a + cy; } double e = cz * cz * 2 - 1; double c = ((-3 * c2a + 4) * F + 4) * c2a * F / 16; double d = x; x = ((e * cy * c + cz) * sy * c + y) * SA; x = (1 - c) * x * F + x2 - x1; if (Math.abs(d - x) <= EPS) { x = Math.sqrt((1 / (R * R) - 1) * c2a + 1) + 1; x = (x - 2) / x; c = 1 - x; c = (x * x / 4 + 1) / c; d = (0.375 * x * x - 1) * x; x = e * cy; s = 1 - 2 * e; s = ((((sy * sy * 4 - 3) * s * cz * d / 6 - x) * d / 4 + cz) * sy * d + y) * c * R * semiMajorAxis; return s; } } final double LEPS = 1.0E-10; if (Math.abs(x1 - x2) <= LEPS && Math.abs(y1 - y2) <= LEPS) { return 0; } if (Math.abs(y1) <= LEPS && Math.abs(y2) <= LEPS) { return Math.abs(x1 - x2) * semiMajorAxis; } throw new ArithmeticException("no convergence"); } } /* * Copyright (C) 2009 Christian Gawron * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License version 2 as * published by the Free Software Foundation. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * * Author: Christian Gawron * Create date: 03-Jul-2009 */ package uk.me.parabola.mkgmap.reader.dem; /** * Find zero of a function using Brent's method. */ public class Brent { public interface Function { public double eval(double x); } static double epsilon = 3.0e-10; public static double zero(Function f, double x1, double x2) { return zero(f, x1, x2, 1e-8, 100); } public static double zero(Function f, double x1, double x2, double tol, int maxit) { double a=x1, b=x2, c=x2, d=0, e=0, min1, min2; double fa=f.eval(a), fb=f.eval(b), fc; double p, q, r, s, tolerance, xm; if ((fa > 0.0 && fb > 0.0) || (fa < 0.0 && fb < 0.0)) throw new ArithmeticException("Root must be bracketed"); fc=fb; for (int iterations=0; iterations < maxit; iterations++) { if ((fb > 0.0 && fc > 0.0) || (fb < 0.0 && fc < 0.0)) { c=a; fc=fa; e=d=b-a; } if (Math.abs(fc) < Math.abs(fb)) { a=b; b=c; c=a; fa=fb; fb=fc; fc=fa; } tolerance=2.0*epsilon*Math.abs(b)+0.5*tol; xm=0.5*(c-b); if (Math.abs(xm) <= tolerance || fb == 0.0) return b; if (Math.abs(e) >= tolerance && Math.abs(fa) > Math.abs(fb)) { s=fb/fa; if (a == c) { p=2.0*xm*s; q=1.0-s; } else { q=fa/fc; r=fb/fc; p=s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0)); q=(q-1.0)*(r-1.0)*(s-1.0); } if (p > 0.0) q = -q; p=Math.abs(p); min1=3.0*xm*q-Math.abs(tolerance*q); min2=Math.abs(e*q); if (2.0*p < (min1 < min2 ? min1 : min2)) { e=d; d=p/q; } else { d=xm; e=d; } } else { d=xm; e=d; } a=b; fa=fb; if (Math.abs(d) > tolerance) b += d; else b += xm >= 0 ? tolerance : -tolerance; fb=f.eval(b); } throw new ArithmeticException("Maximum number of iterations exceeded"); } }

8 23

[PATCH] contours v3
by Christian Gawron 05 Jul '09

05 Jul '09

This version allows to store the contours in separate IMG files using the --dem-separate-img option. The name (8-digit number) of the IMG file will be the original mapname + 10000000. I also replaced orthodromicDistance by quickDistance (I only use the distance calculation to find out if I should refine a contour line segment by splitting it). Best wishes Christian Index: src/uk/me/parabola/mkgmap/Option.java =================================================================== --- src/uk/me/parabola/mkgmap/Option.java (Revision 1078) +++ src/uk/me/parabola/mkgmap/Option.java (Arbeitskopie) @@ -34,7 +34,7 @@ } } - protected Option(String option, String value) { + public Option(String option, String value) { this.option = option; this.value = value; } Index: src/uk/me/parabola/mkgmap/CommandArgsReader.java =================================================================== --- src/uk/me/parabola/mkgmap/CommandArgsReader.java (Revision 1078) +++ src/uk/me/parabola/mkgmap/CommandArgsReader.java (Arbeitskopie) @@ -129,7 +129,7 @@ * @param option The option name. * @param value Its value. */ - private void addOption(String option, String value) { + public void addOption(String option, String value) { CommandOption opt = new CommandOption(option, value); addOption(opt); } @@ -181,13 +181,21 @@ return arglist.iterator(); } + public ArgList getArgList() { + return arglist; + } + + public EnhancedProperties getArgs() { + return args; + } + /** * Read a config file that contains more options. When the number of * options becomes large it is more convenient to place them in a file. * * @param filename The filename to obtain options from. */ - private void readConfigFile(String filename) { + public void readConfigFile(String filename) { Options opts = new Options(new OptionProcessor() { public void processOption(Option opt) { log.debug("incoming opt", opt.getOption(), opt.getValue()); @@ -206,14 +214,14 @@ * the argument to be processed in order. Options can be intersperced with * filenames. The options take effect where they appear. */ - interface ArgType { + public interface ArgType { public abstract void processArg(); } /** * A filename. */ - class Filename implements ArgType { + public class Filename implements ArgType { private final String name; private boolean useFilenameAsMapname = true; @@ -270,7 +278,7 @@ /** * An option argument. A key value pair. */ - class CommandOption implements ArgType { + public class CommandOption implements ArgType { private final Option option; private CommandOption(Option option) { @@ -297,7 +305,7 @@ /** * The arguments are held in this list. */ - class ArgList implements Iterable<CommandArgsReader.ArgType> { + public class ArgList implements Iterable<CommandArgsReader.ArgType> { private final List<ArgType> alist; private int filenameCount; Index: src/uk/me/parabola/mkgmap/reader/MapperBasedMapDataSource.java =================================================================== --- src/uk/me/parabola/mkgmap/reader/MapperBasedMapDataSource.java (Revision 1078) +++ src/uk/me/parabola/mkgmap/reader/MapperBasedMapDataSource.java (Arbeitskopie) @@ -27,12 +27,14 @@ import uk.me.parabola.imgfmt.app.Area; import uk.me.parabola.imgfmt.app.Coord; import uk.me.parabola.imgfmt.app.trergn.Overview; +import uk.me.parabola.mkgmap.general.LoadableMapDataSource; import uk.me.parabola.mkgmap.general.MapDataSource; import uk.me.parabola.mkgmap.general.MapDetails; import uk.me.parabola.mkgmap.general.MapLine; import uk.me.parabola.mkgmap.general.MapPoint; import uk.me.parabola.mkgmap.general.MapShape; import uk.me.parabola.mkgmap.general.RoadNetwork; +import uk.me.parabola.mkgmap.reader.dem.DEM; import uk.me.parabola.util.Configurable; import uk.me.parabola.util.EnhancedProperties; @@ -118,6 +120,11 @@ return configProps; } + public MapDetails getMapper() + { + return mapper; + } + /** * We add the background polygons if the map is not transparent. */ @@ -146,5 +153,8 @@ // the overview section. mapper.addShape(background); } + if (getConfig().getProperty("contours", false)) { + DEM.createContours((LoadableMapDataSource) this, getConfig()); + } } } Index: src/uk/me/parabola/mkgmap/reader/osm/Way.java =================================================================== --- src/uk/me/parabola/mkgmap/reader/osm/Way.java (Revision 1078) +++ src/uk/me/parabola/mkgmap/reader/osm/Way.java (Arbeitskopie) @@ -30,12 +30,18 @@ */ public class Way extends Element { - private final List<Coord> points = new ArrayList<Coord>(); + private final List<Coord> points; public Way(long id) { + points = new ArrayList<Coord>(); setId(id); } + public Way(long id, List<Coord> points) { + this.points = points; + setId(id); + } + /** * Get the points that make up the way. We attempt to re-order the segments * and return a list of points that traces the route of the way. /* * Copyright (C) 2009 Christian Gawron * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License version 2 as * published by the Free Software Foundation. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * * Author: Christian Gawron * Create date: 03-Jul-2009 */ package uk.me.parabola.mkgmap.reader.dem; import java.io.*; import java.nio.channels.FileChannel; import java.nio.MappedByteBuffer; import java.util.Vector; import java.util.List; import java.util.ArrayList; import java.util.Locale; import com.sun.media.jai.codec.*; import javax.media.jai.*; import java.awt.image.renderable.ParameterBlock; import java.awt.image.Raster; import java.awt.image.DataBuffer; import uk.me.parabola.imgfmt.Utils; import uk.me.parabola.imgfmt.ExitException; import uk.me.parabola.imgfmt.FormatException; import uk.me.parabola.imgfmt.FileSystemParam; import uk.me.parabola.imgfmt.app.Area; import uk.me.parabola.imgfmt.app.Coord; import uk.me.parabola.imgfmt.app.map.Map; import uk.me.parabola.mkgmap.build.MapBuilder; import uk.me.parabola.mkgmap.general.MapDetails; import uk.me.parabola.mkgmap.general.MapPoint; import uk.me.parabola.mkgmap.general.MapLine; import uk.me.parabola.mkgmap.general.LevelInfo; import uk.me.parabola.mkgmap.general.LoadableMapDataSource; import uk.me.parabola.mkgmap.reader.MapperBasedMapDataSource; import uk.me.parabola.mkgmap.reader.osm.OsmConverter; import uk.me.parabola.mkgmap.reader.osm.Style; import uk.me.parabola.mkgmap.reader.osm.Way; import uk.me.parabola.mkgmap.osmstyle.StyleImpl; import uk.me.parabola.mkgmap.osmstyle.StyledConverter; import uk.me.parabola.mkgmap.osmstyle.eval.SyntaxException; import uk.me.parabola.util.EnhancedProperties; import static java.nio.channels.FileChannel.MapMode.READ_ONLY; /** * Create contour lines using an algorithm similar to that described in * <a href="http://mapcontext.com/autocarto/proceedings/auto-carto-5/pdf/an-adaptive-gr…">An Adaptive Grid Contouring Algorithm</a> by Downing and Zoraster. */ public abstract class DEM { final static double epsilon = 1e-9; final static double delta = 1.5; final static int maxPoints=200000; final static double minDist = 15; final static double maxDist = 21; final static double step = 0.01; final static double semiMajorAxis = 6378137.0; final static double inverseFlattening = 298.257223563; static int M=1200; static int N=M; static double res = 1.0/N; static int id = -1; short values[] = null; int lat; int lon; abstract double ele(int x, int y); abstract void read(int minLon, int minLat, int maxLon, int maxLat); abstract void serializeCopyRight(Writer out) throws IOException; public static void createContours(LoadableMapDataSource mapData, EnhancedProperties config) { Area bounds = mapData.getBounds(); double minLat = Utils.toDegrees(bounds.getMinLat()); double minLon = Utils.toDegrees(bounds.getMinLong()); double maxLat = Utils.toDegrees(bounds.getMaxLat()); double maxLon = Utils.toDegrees(bounds.getMaxLong()); System.out.printf("bounds: %f %f %f %f\n", minLat, minLon, maxLat, maxLon); DEM data; String demType = config.getProperty("dem-type", "CGIAR"); String dataPath; if (demType.equals("ASTER")) { dataPath = config.getProperty("dem-path", "ASTER"); data = new ASTERGeoTiff(dataPath, minLat, minLon, maxLat, maxLon); } else if (demType.equals("CGIAR")) { dataPath = config.getProperty("dem-path", "CGIAR"); data = new CGIARGeoTiff(dataPath, minLat, minLon, maxLat, maxLon); } else { dataPath = config.getProperty("dem-path", "SRTM"); data = new HGTDEM(dataPath, minLat, minLon, maxLat, maxLon); } Isolines lines = data.new Isolines(data, minLat, minLon, maxLat, maxLon); int increment = config.getProperty("dem-increment", 10); double minHeight = lines.getMinHeight(); double maxHeight = lines.getMaxHeight(); int maxLevels = config.getProperty("dem-maxlevels", 100); while ((maxHeight-minHeight)/increment > maxLevels) increment *= 2; String loc = config.getProperty("style-file"); if (loc == null) loc = config.getProperty("map-features"); String name = config.getProperty("style"); if (loc == null && name == null) name = "default"; LoadableMapDataSource dest = mapData; if (config.getProperty("dem-separate-img", false)) { dest = new DEMMapDataSource(mapData, config); } OsmConverter converter; try { Style style = new StyleImpl(loc, name); style.applyOptionOverride(config); converter = new StyledConverter(style, ((MapperBasedMapDataSource) dest).getMapper(), config); } catch (SyntaxException e) { System.err.println("Error in style: " + e.getMessage()); throw new ExitException("Could not open style " + name); } catch (FileNotFoundException e) { String name1 = (name != null)? name: loc; throw new ExitException("Could not open style " + name1); } for (int level=0; level<maxHeight; level+=increment) { if (level < minHeight) continue; // create isolines lines.addLevel(level); for (Isolines.Isoline line : lines.isolines) { Way way = new Way(id--, line.points); way.addTag("contour", "elevation"); way.addTag("ele", String.format("%d", (int) line.level)); converter.convertWay(way); } lines.isolines.clear(); } if (config.getProperty("dem-separate-img", false)) { MapBuilder builder = new MapBuilder(); builder.config(config); FileSystemParam params = new FileSystemParam(); params.setMapDescription("contour lines"); long mapName = Integer.valueOf(config.getProperty("mapname", "63240000")); try { Map map = Map.createMap(String.format("%08d", mapName+10000000), params); builder.makeMap(map, dest); map.close(); } catch (Exception ex) { throw new RuntimeException(ex); } } } public static class HGTDEM extends DEM { MappedByteBuffer buffer = null; public HGTDEM(String dataPath, double minLat, double minLon, double maxLat, double maxLon) { this.lat = (int) minLat; this.lon = (int) minLon; if (maxLat > lat+1 || maxLon > lon+1) throw new RuntimeException("Area too large (must not span more than one SRTM file)"); String northSouth = lat < 0 ? "S" : "N"; String eastWest = lon > 0 ? "E" : "W"; String fileName = String.format("%s/%s%02d%s%03d.hgt", dataPath, northSouth, lat < 0 ? -lat : lat, eastWest, lon < 0 ? -lon : lon); try { FileInputStream is = new FileInputStream(fileName); buffer = is.getChannel().map(READ_ONLY, 0, 2*(M+1)*(M+1)); } catch (Exception e) { throw new RuntimeException(e); } } public void read(int minLon, int minLat, int maxLon, int maxLat) { } public double ele(int x, int y) { return buffer.getShort(2*((M-y)*(M+1)+x))+delta; } public void serializeCopyRight(Writer out) throws IOException { out.write(" <copyright>\n"); out.write(" Contour lines generated from DEM data by NASA\n"); out.write(" </copyright>\n"); } } public static class CGIARGeoTiff extends DEM { Raster raster; String fileName; int minLat, minLon, maxLat, maxLon; PlanarImage image; public CGIARGeoTiff(String dataPath, double minLat, double minLon, double maxLat, double maxLon) { this.lat = ((int) (minLat/5))*5; this.lon = ((int) (minLon/5))*5; if (maxLat > lat+5 || maxLon > lon+5) throw new RuntimeException("Area too large (must not span more than one CGIAR GeoTIFF)"); int tileX, tileY; tileX = (180 + lon)/5 + 1; tileY = (60 - lat)/5; this.fileName = String.format("%s/srtm_%02d_%02d.tif", dataPath, tileX, tileY); init(); } public void serializeCopyRight(Writer out) throws IOException { out.write(" <copyright>\n"); out.write(" Contour lines generated from improved SRTM data by CIAT-CSI (see http://srtm.csi.cgiar.org)\n"); out.write(" </copyright>\n"); } public void read(int minLon, int minLat, int maxLon, int maxLat) { this.minLon = minLon; this.minLat = minLat; this.maxLon = maxLon; this.maxLat = maxLat; raster = image.getData(new java.awt.Rectangle(minLon, 6000-maxLat-1, maxLon-minLon+1, maxLat-minLat+1)); System.out.printf("read: %d %d %d %d\n", minLon, 6000-maxLat-1, maxLon-minLon+1, maxLat-minLat+1); } void init() { System.out.printf("CGIAR GeoTIFF: %s\n", fileName); N = 6000; M = 6000; res = 5.0/M; try { SeekableStream s = new FileSeekableStream(fileName); ParameterBlock pb = new ParameterBlock(); pb.add(s); TIFFDecodeParam param = new TIFFDecodeParam(); pb.add(param); RenderedOp op = JAI.create("tiff", pb); image = op.createInstance(); System.out.printf("Image: %d %d %d %d\n", image.getWidth(), image.getHeight(), image.getNumXTiles(), image.getNumYTiles()); } catch (Exception e) { throw new RuntimeException(e); } } public double ele(int x, int y) { try { int elevation = raster.getPixel(x, 6000-y-1, (int[])null)[0]; return elevation+delta; } catch (ArrayIndexOutOfBoundsException ex) { System.out.printf("ele: (%d, %d) (%d, %d, %d, %d) %s\n", x, 6000-y-1, raster.getMinX(), raster.getMinY(), raster.getWidth(), raster.getHeight(), ex.toString()); throw ex; } } } public static class ASTERGeoTiff extends DEM { Raster raster; String fileName; int minLat, minLon, maxLat, maxLon; PlanarImage image; public ASTERGeoTiff(String dataPath, double minLat, double minLon, double maxLat, double maxLon) { this.lat = (int) minLat; this.lon = (int) minLon; if (maxLat > lat+1 || maxLon > lon+1) throw new RuntimeException("Area too large (must not span more than one ASTER GeoTIFF)"); String northSouth = lat < 0 ? "S" : "N"; String eastWest = lon > 0 ? "E" : "W"; fileName = String.format("%s/ASTGTM_%s%02d%s%03d_dem.tif", dataPath, northSouth, lat < 0 ? -lat : lat, eastWest, lon < 0 ? -lon : lon); init(); } public void serializeCopyRight(Writer out) throws IOException { out.write(" <copyright>\n"); out.write(" Contour lines generated from DGM data by ASTER (see https://wist.echo.nasa.gov/~wist/api/imswelcome)\n"); out.write(" </copyright>\n"); } public void read(int minLon, int minLat, int maxLon, int maxLat) { this.minLon = minLon; this.minLat = minLat; this.maxLon = maxLon; this.maxLat = maxLat; raster = image.getData(new java.awt.Rectangle(minLon, 3601-maxLat-1, maxLon-minLon+1, maxLat-minLat+1)); System.out.printf("read: %d %d %d %d\n", minLon, 3601-maxLat-1, maxLon-minLon+1, maxLat-minLat+1); } void init() { System.out.printf("ASTER GeoTIFF: %s\n", fileName); N = 3600; M = 3600; res = 1.0/M; try { SeekableStream s = new FileSeekableStream(fileName); ParameterBlock pb = new ParameterBlock(); pb.add(s); TIFFDecodeParam param = new TIFFDecodeParam(); pb.add(param); RenderedOp op = JAI.create("tiff", pb); image = op.createInstance(); System.out.printf("Image: %d %d %d %d\n", image.getWidth(), image.getHeight(), image.getNumXTiles(), image.getNumYTiles()); } catch (Exception e) { throw new RuntimeException(e); } } public double ele(int x, int y) { try { int elevation = raster.getPixel(x, 3601-y-1, (int[])null)[0]; return elevation+delta; } catch (ArrayIndexOutOfBoundsException ex) { System.out.printf("ele: (%d, %d) (%d, %d, %d, %d) %s\n", x, 3601-y-1, raster.getMinX(), raster.getMinY(), raster.getWidth(), raster.getHeight(), ex.toString()); throw ex; } } } private static void debug(String format, Object ... args) { //System.out.printf(format + "\n", args); } int lastXi = -1; int lastYi = -1; final static int bcInv[][]= { { 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, { 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0}, { -3, 0, 0, 3, 0, 0, 0, 0, -2, 0, 0, -1, 0, 0, 0, 0}, { 2, 0, 0, -2, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0}, { 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0}, { 0, 0, 0, 0, -3, 0, 0, 3, 0, 0, 0, 0, -2, 0, 0, -1}, { 0, 0, 0, 0, 2, 0, 0, -2, 0, 0, 0, 0, 1, 0, 0, 1}, { -3, 3, 0, 0, -2, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, { 0, 0, 0, 0, 0, 0, 0, 0, -3, 3, 0, 0, -2, -1, 0, 0}, { 9, -9, 9, -9, 6, 3, -3, -6, 6, -6, -3, 3, 4, 2, 1, 2}, { -6, 6, -6, 6, -4, -2, 2, 4, -3, 3, 3, -3, -2, -1, -1, -2}, { 2, -2, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, { 0, 0, 0, 0, 0, 0, 0, 0, 2, -2, 0, 0, 1, 1, 0, 0}, { -6, 6, -6, 6, -3, -3, 3, 3, -4, 4, 2, -2, -2, -2, -1, -1}, { 4, -4, 4, -4, 2, 2, -2, -2, 2, -2, -2, 2, 1, 1, 1, 1} }; static int lastId=1000000000; static double lastX = 0; static double lastY = 0; static int edge[] = new int[2]; static int off0[][] = {{ 0, 0}, { 0, 0}, { 0, 1}, { 1, 1}}; static int off1[][] = {{ 1, 0}, { 0, 1}, { 1, 1}, { 1, 0}}; static int brd[] = { 1, 2, 4, 8 }; static int inv[] = { 4, 8, 1, 2 }; static int rev[] = { 2, 3, 0, 1 }; static int mov[][] = {{ 0, -1}, {-1, 0}, { 0, 1}, { 1, 0}}; double bc[][] = new double[4][4]; double bc_y[] = new double[4]; double bc_y1[] = new double[4]; double bc_y2[] = new double[4]; double bc_y12[] = new double[4]; double bc_Coeff[] = new double[16]; double bc_x[] = new double[16]; private void recalculateCoefficients(int xi, int yi) { double v00, vp0, v0p, vpp; double vm0, v0m, vpm, vmp, vmm, vmP, vPm; double vP0, v0P, vPp, vpP, vPP; v00 = ele(xi, yi); v0p = ele(xi, yi+1); vpp = ele(xi+1, yi+1); vp0 = ele(xi+1, yi); vm0 = ele(xi-1, yi); v0m = ele(xi, yi-1); vmp = ele(xi-1, yi+1); vpm = ele(xi+1, yi-1); vmm = ele(xi-1, yi-1); vmP = ele(xi+2, yi-1); vPm = ele(xi-1, yi+2); vP0 = ele(xi+2, yi); v0P = ele(xi, yi+2); vPp = ele(xi+2, yi+1); vpP = ele(xi+1, yi+2); vPP = ele(xi+2, yi+2); bc_y[0] = v00; bc_y[1] = vp0; bc_y[2] = vpp; bc_y[3] = v0p; bc_y1[0] = (vp0-vm0)/2; bc_y1[1] = (vP0-v00)/2; bc_y1[2] = (vPp-v0p)/2; bc_y1[3] = (vpp-vmp)/2; bc_y2[0] = (v0p-v0m)/2; bc_y2[1] = (vpp-vpm)/2; bc_y2[2] = (vpP-vp0)/2; bc_y2[3] = (v0P-v00)/2; bc_y12[0] = (vpp - vpm - vmp + vmm) / 4; bc_y12[0] = (vPp - vPm - v0p + v0m) / 4; bc_y12[2] = (vPP - vP0 - v0P + v00) / 4; bc_y12[0] = (vpP - vp0 - vmP + vm0) / 4; int j, i; double s; for (i=0; i<4; i++) { bc_x[i]=bc_y[i]; bc_x[i+4]=bc_y1[i]; bc_x[i+8]=bc_y2[i]; bc_x[i+12]=bc_y12[i]; } for (i=0; i<16; i++) { s = 0; for (int k=0; k<16; k++) s += bcInv[i][k]*bc_x[k]; bc_Coeff[i] = s; } int l = 0; for (i=0; i<4; i++) for (j=0; j<4; j++) bc[i][j] = bc_Coeff[l++]; } public double gradient(double lat, double lon, double[] grad) { grad[0] = 0; grad[1] = 0; double x = (lon-this.lon)/res; double y = (lat-this.lat)/res; int xi = (int) x; int yi = (int) y; if (lastXi != xi || lastYi != yi) { debug("new Cell for interpolation: %d %d", xi, yi); recalculateCoefficients(xi, yi); lastXi = xi; lastYi = yi; } double t = x - xi; double u = y - yi; if (xi < 0 || xi > N+1 || yi < 0 || yi > N+1) throw new IndexOutOfBoundsException(String.format("(%f, %f)->(%d, %d)", lat, lon, xi, yi)); double val = 0; for (int i=3; i>=0; i--) { val = t*val + ((bc[i][3]*u + bc[i][2])*u + bc[i][1])*u + bc[i][0]; grad[0] = u*grad[0] + (3*bc[3][i]*t + 2*bc[2][i])*t + bc[1][i]; grad[1] = t*grad[1] + (3*bc[i][3]*t + 2*bc[i][2])*t + bc[i][1]; } return val; } public double elevation(double lat, double lon) { double x = (lon-this.lon)/res; double y = (lat-this.lat)/res; int xi = (int) x; int yi = (int) y; if (lastXi != xi || lastYi != yi) { debug("new Cell for interpolation: %d %d", xi, yi); recalculateCoefficients(xi, yi); lastXi = xi; lastYi = yi; } double t = x - xi; double u = y - yi; if (xi < 0 || xi > N+1 || yi < 0 || yi > N+1) throw new IndexOutOfBoundsException(String.format("(%f, %f)->(%d, %d)", lat, lon, xi, yi)); double val = 0; for (int i=3; i>=0; i--) { val = t*val + ((bc[i][3]*u + bc[i][2])*u + bc[i][1])*u + bc[i][0]; } return val; } public double elevation(int x, int y) { if (x < 0 || x > N || y < 0 || y > N) throw new IndexOutOfBoundsException(String.format("elevation: %d %d", x, y)); return ele(x, y); } class Isolines { DEM data; int minX; int maxX; int minY; int maxY; double min; double max; Vector<Isoline> isolines = new Vector<Isoline>(); class Isoline { int id; Vector<Coord> points; double level; boolean isClosed; private Isoline(double level) { this.level = level; isClosed = false; id = lastId++; points = new Vector<Coord>(); } private class Edge implements Brent.Function { double x0, y0, x1, y1; Edge(double x0, double y0, double x1, double y1) { this.x0 = x0; this.y0 = y0; this.x1 = x1; this.y1 = y1; } public double eval(double d) { double f = data.elevation(x0 + d * (x1-x0), y0 + d * (y1-y0)) - level; //System.out.printf("evalEdge: %f %f\n", d, f); return f; } } private class FN implements Brent.Function { double x0, y0; double dx, dy; public void setParameter(double x0, double y0, double dx, double dy) { this.x0 = x0; this.y0 = y0; this.dx = dx; this.dy = dy; } public double eval(double t) { double f = data.elevation(y0+t*dy, x0+t*dx) - level; return f; } } private FN fn = new FN(); double grad[] = new double[2]; double px[] = new double[4]; double py[] = new double[4]; int edges[] = new int[4]; boolean addCell(Position p, int direction) { debug("addCell: %f %d %d %d %d", level, p.ix, p.iy, p.edge, direction); int c = 0; for (int k=0; k<4; k++) { if (k == p.edge) continue; int x0 = p.ix + off0[k][0]; int y0 = p.iy + off0[k][1]; int x1 = p.ix + off1[k][0]; int y1 = p.iy + off1[k][1]; double l0 = elevation(x0, y0) - level; double l1 = elevation(x1, y1) - level; if (Math.abs(l1) < epsilon || l0*l1 < 0) { edges[c] = k; Brent.Function f = new Edge(data.lat + y0 * DEM.res, data.lon + x0 * DEM.res, data.lat + y1 * DEM.res, data.lon + x1 * DEM.res); double f0 = elevation(x0, y0) - level; double f1 = elevation(x1, y1) - level; double delta; if (Math.abs(1) < epsilon) { delta = 1; } else if (Math.abs(f0) < epsilon) throw new RuntimeException("implementation error!"); else delta = Brent.zero(f, epsilon, 1-epsilon); px[c] = data.lon + (x0+delta*(x1-x0))*DEM.res; py[c] = data.lat + (y0+delta*(y1-y0))*DEM.res; c++; } } if (c == 1) { p.edge = edges[0]; double px0 = p.x; double py0 = p.y; p.x = px[0]; p.y = py[0]; double px1 = p.x; double py1 = p.y; double xMin = data.lon + p.ix * DEM.res; double xMax = xMin + DEM.res; double yMin = data.lat + p.iy * DEM.res; double yMax = yMin + DEM.res; refineAdaptively(xMin, yMin, xMax, yMax, px0, py0, px1, py1, direction, maxDist); addPoint(p.x, p.y, direction); p.moveCell(); return true; } else { debug("addCellByStepping: %d", c); return addCellByStepping(p, direction, c, edges, px, py); } } private void refineAdaptively(double xMin, double yMin, double xMax, double yMax, double x0, double y0, double x1, double y1, int direction, double maxDist) { double dist = quickDistance(x0, y0, x1, y1); if (dist > maxDist) { double dx = x1-x0; double dy = y1-y0; double xm, ym, f0, f1, t0, t1, t; Brent.Function f; xm = x0 + 0.5*dx; ym = y0 + 0.5*dy; double n = Math.sqrt(dx*dx+dy*dy); fn.setParameter(xm, ym, -dy/n, dx/n); f = fn; t0 = -0.05*res; t1 = 0.05*res; f0 = f.eval(t0); f1 = f.eval(t1); int count = 0; while (f0 * f1 > 0 && count++ < 20) { if ((count & 1) > 0) t0 -= 0.05*res; else t1 += 0.05*res; f0 = f.eval(t0); f1 = f.eval(t1); debug("refine: %f %f %f %f", t0, t1, f0, f1); } if (f0 * f1 < 0) { t = Brent.zero(f, t0, t1); xm -= t*dy; ym += t*dx; } else { debug("refine failed: %f %f %f %f", t0, t1, f0, f1); if (false) throw new RuntimeException(String.format("refine failed: %f %f %f %f %f %f %f %f", xMin, yMin, xMax, yMax, x0, y0, x1, y1)); return; } if (xm > xMin && xm < xMax && ym > yMin && ym < yMax) refineAdaptively(xMin, yMin, xMax, yMax, x0, y0, xm, ym, direction, maxDist*1.1); addPoint(xm, ym, direction); if (xm > xMin && xm < xMax && ym > yMin && ym < yMax) refineAdaptively(xMin, yMin, xMax, yMax, xm, ym, x1, y1, direction, maxDist*1.1); } } boolean addCellByStepping(Position p, int direction, int numEdges, int[] edges, double[] px, double[] py) { debug("addCellByStepping: %f %d %d %d %d", level, p.ix, p.iy, p.edge, direction); double xMin = data.lon + p.ix * DEM.res; double xMax = xMin + DEM.res; double yMin = data.lat + p.iy * DEM.res; double yMax = yMin + DEM.res; double dt, t0, t1; double f0, f1; boolean edgeHit = false; double h[] = new double[4]; int dir; double n2 = Math.sqrt(1.0/(grad[0]*grad[0] + grad[1]*grad[1])); double dx; double dy; int count=0; int iMin = -1; double md = 5000; for (int i=0; i<numEdges; i++) { gradient(p.y, p.x, grad); double dist = quickDistance(p.x, p.y, px[i], py[i]); debug("distance %d: %f", i, dist); if (dist < md && (visited[p.iy*(N+1)+p.ix] & brd[edges[i]]) == 0) { md = dist; iMin = i; } } p.edge = edges[iMin]; double px0 = p.x; double py0 = p.y; p.x = px[iMin]; p.y = py[iMin]; double px1 = p.x; double py1 = p.y; xMin = data.lon + p.ix * DEM.res; xMax = xMin + DEM.res; yMin = data.lat + p.iy * DEM.res; yMax = yMin + DEM.res; refineAdaptively(xMin, yMin, xMax, yMax, px0, py0, px1, py1, direction, maxDist); addPoint(p.x, p.y, direction); p.moveCell(); return true; } private void addPoint(double x, double y, int direction) { double dist = quickDistance(x, y, lastX, lastY); debug("addPoint: %f %f %f", x, y, dist); if (dist > minDist) { if (direction > 0) points.add(0, new Coord(y, x)); else points.add(points.size(), new Coord(y, x)); lastX = x; lastY = y; } } private void close() { points.add(points.size(), points.get(0)); isClosed = true; } private void addMove(int x0, int y0, int x1, int y1, int direction) { double x; double y; if (x0 < minX || x0 >= maxX || y0 < minY || y0 >= maxY) return; double l0 = data.elevation(x0, y0); double l1 = data.elevation(x1, y1); debug("addMove: %d %d %d %d %.2f %.2f %d", x0, y0, x1, y1, l0, l1, direction); if ((l0 < level && l1 < level) || (l0 > level && l1 > level)) throw new RuntimeException("implementation error"); if (l0 == level) { x = data.lon + x0 * DEM.res; y = data.lat + y0 * DEM.res; } else { double delta = (l0-level)/(l0-l1); x = data.lon + (x0 + delta * (x1-x0)) * DEM.res; y = data.lat + (y0 + delta * (y1-y0)) * DEM.res; } double dist = quickDistance(x, y, lastX, lastY); debug("levels: %d %d %f, %d %d %f: %f %f %f", x0, y0, l0, x1, y1, l1, x, y, dist); if (dist > 1) { if (direction > 0) points.add(0, new Coord(y, x)); else points.add(points.size(), new Coord(y, x)); lastX = x; lastY = y; } } } public Isolines(DEM data, int minX, int maxX, int minY, int maxY) { this.data = data; this.minX = minX; this.maxX = maxX; this.minY = minY; this.maxY = maxY; init(); } public Isolines(DEM data, double minLat, double minLon, double maxLat, double maxLon) { System.out.printf("init: %f %f %f %f\n", minLat, minLon, maxLat, maxLon); this.data = data; this.minX = (int) ((minLon-data.lon)/data.res); this.minY = (int) ((minLat-data.lat)/data.res); this.maxX = (int) ((maxLon-data.lon)/data.res); this.maxY = (int) ((maxLat-data.lat)/data.res); init(); } private void init() { System.out.printf("init: %d %d %d %d\n", minX, minY, maxX, maxY); data.read(minX-2, minY-2, maxX+2, maxY+2); // we need some overlap for bicubic interpolation max = -1000; min = 10000; for (int i=minX; i<maxX; i++) for (int j=minY; j<maxY; j++) { if (data.elevation(i, j) < min) min = data.elevation(i, j); if (data.elevation(i, j) > max) max = data.elevation(i, j); } debug("min: %f, max: %f\n", min, max); } double getMinHeight() { return min; } double getMaxHeight() { return max; } public void addLevels(int start, int increment) { for (int level=start; level<max; level+=increment) addLevel(level); } private class Edge implements Brent.Function { double x0, y0, x1, y1, level; Edge(double x0, double y0, double x1, double y1, double level) { this.x0 = x0; this.y0 = y0; this.x1 = x1; this.y1 = y1; this.level = level; } public double eval(double d) { double f = data.elevation(x0 + d * (x1-x0), y0 + d * (y1-y0)) - level; //System.out.printf("evalEdge: %f %f\n", d, f); return f; } } private class Position { int ix, iy; double x, y; int edge; Position(int ix, int iy, double x, double y, int edge) { this.ix = ix; this.iy = iy; this.x = x; this.y = y; this.edge = edge; } Position(Position p) { this.ix = p.ix; this.iy = p.iy; this.x = p.x; this.y = p.y; this.edge = p.edge; } void markEdge() { debug("marking edge: %d %d %d %d", ix, iy, edge, brd[edge]); visited[iy*(N+1)+ix] |= brd[edge]; } void moveCell() { markEdge(); ix += mov[edge][0]; iy += mov[edge][1]; edge = rev[edge]; markEdge(); } } byte visited[] = new byte[(N+1)*(N+1)]; public void addLevel(double level) { if (level < min || level > max) return; System.out.printf("addLevel: %f\n", level); java.util.Arrays.fill(visited, (byte) 0); for (int y=minY; y<maxY; y++) { for (int x=minX; x<maxX; x++) { byte v = 0; // Mark the borders of the cell, represented by the four points (i, j), (i+1, j), (i, j+1), (i+1, j+1), // which are intersected by the contour. The values are: // 1: top // 2: left // 4: bottom // 8: right if (data.elevation(x, y) >= level) { if (data.elevation(x+1, y) < level) { v |= 1; } if (data.elevation(x, y+1) < level) { v |= 2; } } else { if (data.elevation(x+1, y) > level) { v |= 1; } if (data.elevation(x, y+1) > level) { v |= 2; } } int k=-1; if ((v&1) > 0 && (visited[y*(N+1)+x]&1) == 0) { k=0; } else if ((v&2) > 0 && (visited[y*(N+1)+x]&2) == 0) { k=1; } if (k>=0) { int x0 = x + off0[k][0]; int y0 = y + off0[k][1]; int x1 = x + off1[k][0]; int y1 = y + off1[k][1]; try { Brent.Function f = new Edge(data.lat + y0 * DEM.res, data.lon + x0 * DEM.res, data.lat + y1 * DEM.res, data.lon + x1 * DEM.res, level); double f0 = elevation(x0, y0) - level; double f1 = elevation(x1, y1) - level; double delta; if (Math.abs(f0) < epsilon) { delta = 0; } else if (Math.abs(f1) < epsilon) continue; else delta = Brent.zero(f, 0, 1-epsilon); Position p = new Position(x, y, data.lon + (x0+delta*(x1-x0))*DEM.res, data.lat + (y0+delta*(y1-y0))*DEM.res, k); p.markEdge(); isolines.add(traceByStepping(level, p)); } catch (RuntimeException ex) { debug("error: %s", ex.toString()); ex.printStackTrace(); return; } } } } } private Isoline traceByStepping(double level, Position p) { debug("traceByStepping: starting contour %f %d %d %f %f %d", level, p.ix, p.iy, p.x, p.y, p.edge); int direction = 1; int n = 0; Position startP = new Position(p); boolean foundEnd = false; Isoline line = new Isoline(level); while (true) { debug("traceByStepping: %f %d %d %f %f %d", level, p.ix, p.iy, p.x, p.y, p.edge); visited[p.iy*(N+1)+p.ix] |= brd[p.edge]; if (n>0 && p.ix == startP.ix && p.iy == startP.iy && quickDistance(p.x, p.y, startP.x, startP.y) < 5) { debug("closed curve!"); line.close(); break; } else if (p.ix < minX || p.iy < minY || p.ix >= maxX || p.iy >= maxY) { if (foundEnd) // did we already reach one end? { debug("second border reached!"); break; } else { debug("border reached!"); foundEnd = true; n = 0; direction *= -1; p = new Position(startP); p.moveCell(); continue; } } n++; if (!line.addCell(p, direction) || line.points.size() > maxPoints) { debug("ending contour"); isolines.add(line); return line; } } return line; } } public static double quickDistance(double long1, double lat1, double long2, double lat2) { double latDiff; if (lat1 < lat2) latDiff = lat2 - lat1; else latDiff = lat1 - lat2; if (latDiff > 90) latDiff -= 180; double longDiff; if (long1 < long2) longDiff = long2 - long1; else longDiff = long1 - long2; if (longDiff > 180) longDiff -= 360; // scale longDiff by cosine of average latitude longDiff *= Math.cos(Math.PI / 180 * Math.abs((lat1 + lat2) / 2)); double distDegSq = (latDiff * latDiff) + (longDiff * longDiff); return 40075000 * Math.sqrt(distDegSq) / 360; } /** * Returns the orthodromic distance between two geographic coordinates in WGS84 datum. * The orthodromic distance is the shortest distance between two points * on a sphere's surface. The orthodromic path is always on a great circle. * * @param x1 Longitude of first point (in degrees). * @param y1 Latitude of first point (in degrees). * @param x2 Longitude of second point (in degrees). * @param y2 Latitude of second point (in degrees). * @return The orthodromic distance (in meters). * */ public static double orthodromicDistance(double x1, double y1, double x2, double y2) { x1 = Math.toRadians(x1); y1 = Math.toRadians(y1); x2 = Math.toRadians(x2); y2 = Math.toRadians(y2); final int MAX_ITERATIONS = 100; final double EPS = 5.0E-14; final double F = 1 / inverseFlattening; final double R = 1 - F; double tu1 = R * Math.sin(y1) / Math.cos(y1); double tu2 = R * Math.sin(y2) / Math.cos(y2); double cu1 = 1 / Math.sqrt(tu1 * tu1 + 1); double cu2 = 1 / Math.sqrt(tu2 * tu2 + 1); double su1 = cu1 * tu1; double s = cu1 * cu2; double baz = s * tu2; double faz = baz * tu1; double x = x2 - x1; for (int i = 0; i < MAX_ITERATIONS; i++) { final double sx = Math.sin(x); final double cx = Math.cos(x); tu1 = cu2 * sx; tu2 = baz - su1 * cu2 * cx; final double sy = Math.sqrt(tu1*tu1 + tu2*tu2);; final double cy = s * cx + faz; final double y = Math.atan2(sy, cy); final double SA = s * sx / sy; final double c2a = 1 - SA * SA; double cz = faz + faz; if (c2a > 0) { cz = -cz / c2a + cy; } double e = cz * cz * 2 - 1; double c = ((-3 * c2a + 4) * F + 4) * c2a * F / 16; double d = x; x = ((e * cy * c + cz) * sy * c + y) * SA; x = (1 - c) * x * F + x2 - x1; if (Math.abs(d - x) <= EPS) { x = Math.sqrt((1 / (R * R) - 1) * c2a + 1) + 1; x = (x - 2) / x; c = 1 - x; c = (x * x / 4 + 1) / c; d = (0.375 * x * x - 1) * x; x = e * cy; s = 1 - 2 * e; s = ((((sy * sy * 4 - 3) * s * cz * d / 6 - x) * d / 4 + cz) * sy * d + y) * c * R * semiMajorAxis; return s; } } final double LEPS = 1.0E-10; if (Math.abs(x1 - x2) <= LEPS && Math.abs(y1 - y2) <= LEPS) { return 0; } if (Math.abs(y1) <= LEPS && Math.abs(y2) <= LEPS) { return Math.abs(x1 - x2) * semiMajorAxis; } throw new ArithmeticException("no convergence"); } private static class DEMMapDataSource extends MapperBasedMapDataSource implements LoadableMapDataSource { LoadableMapDataSource parent; List<String> copyright = new ArrayList<String>(); DEMMapDataSource(LoadableMapDataSource parent, EnhancedProperties props) { this.parent = parent; config(props); } public boolean isFileSupported(String name) { return false; } public void load(String name) throws FileNotFoundException, FormatException { throw new FormatException("load not supported"); } public LevelInfo[] mapLevels() { return parent.mapLevels(); } public String[] copyrightMessages() { return copyright.toArray(new String[1]); } } } /* * Copyright (C) 2009 Christian Gawron * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License version 2 as * published by the Free Software Foundation. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * * Author: Christian Gawron * Create date: 03-Jul-2009 */ package uk.me.parabola.mkgmap.reader.dem; /** * Find zero of a function using Brent's method. */ public class Brent { public interface Function { public double eval(double x); } static double epsilon = 3.0e-10; public static double zero(Function f, double x1, double x2) { return zero(f, x1, x2, 1e-8, 100); } public static double zero(Function f, double x1, double x2, double tol, int maxit) { double a=x1, b=x2, c=x2, d=0, e=0, min1, min2; double fa=f.eval(a), fb=f.eval(b), fc; double p, q, r, s, tolerance, xm; if ((fa > 0.0 && fb > 0.0) || (fa < 0.0 && fb < 0.0)) throw new ArithmeticException("Root must be bracketed"); fc=fb; for (int iterations=0; iterations < maxit; iterations++) { if ((fb > 0.0 && fc > 0.0) || (fb < 0.0 && fc < 0.0)) { c=a; fc=fa; e=d=b-a; } if (Math.abs(fc) < Math.abs(fb)) { a=b; b=c; c=a; fa=fb; fb=fc; fc=fa; } tolerance=2.0*epsilon*Math.abs(b)+0.5*tol; xm=0.5*(c-b); if (Math.abs(xm) <= tolerance || fb == 0.0) return b; if (Math.abs(e) >= tolerance && Math.abs(fa) > Math.abs(fb)) { s=fb/fa; if (a == c) { p=2.0*xm*s; q=1.0-s; } else { q=fa/fc; r=fb/fc; p=s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0)); q=(q-1.0)*(r-1.0)*(s-1.0); } if (p > 0.0) q = -q; p=Math.abs(p); min1=3.0*xm*q-Math.abs(tolerance*q); min2=Math.abs(e*q); if (2.0*p < (min1 < min2 ? min1 : min2)) { e=d; d=p/q; } else { d=xm; e=d; } } else { d=xm; e=d; } a=b; fa=fb; if (Math.abs(d) > tolerance) b += d; else b += xm >= 0 ? tolerance : -tolerance; fb=f.eval(b); } throw new ArithmeticException("Maximum number of iterations exceeded"); } }

1 0

[PATCH] contour lines v2
by Christian Gawron 04 Jul '09

04 Jul '09

Dear all, thank you for your feedback! I changed two points: - Using the option --dem-maxlevels you can specify the maximum number of levels created (the increment is doubled untils this restriction holds). The default is 100. - Now the OSM styling engine is used to set labels, types & resolution of the contour lines. > output osm format data from the DEM (for other use, or later use, or > for intermediate reprocessing) > Actually this was the way my code worked before, but I don't like the size of the created osm files. > make an img that's transparent with the DEM data, not only for > licensing reasons but also because it really seems like something the > user would want to flip off sometimes. I would like to implement that - can one of the experts please give me a hint how I create an IMG file from a MapDataSource? Right now, I plugged my code into the addBackground method of MapperBasedMapDataSource. It would be now problem to create a MapDetails object containing only the contour lines, but how can I write this into an IMG file? I looked at the MapBuilder code, and it seems the makeMap method there needs a LoadableMapDataSource. But I don't want to implement a load(String fileName) method ... Best wishes Christian /* * Copyright (C) 2009 Christian Gawron * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License version 2 as * published by the Free Software Foundation. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * * Author: Christian Gawron * Create date: 03-Jul-2009 */ package uk.me.parabola.mkgmap.reader.dem; import java.io.*; import java.nio.channels.FileChannel; import java.nio.MappedByteBuffer; import java.util.Vector; import java.util.Locale; import com.sun.media.jai.codec.*; import javax.media.jai.*; import java.awt.image.renderable.ParameterBlock; import java.awt.image.Raster; import java.awt.image.DataBuffer; import uk.me.parabola.imgfmt.Utils; import uk.me.parabola.imgfmt.ExitException; import uk.me.parabola.imgfmt.FormatException; import uk.me.parabola.imgfmt.app.Area; import uk.me.parabola.imgfmt.app.Coord; import uk.me.parabola.mkgmap.general.MapDetails; import uk.me.parabola.mkgmap.general.MapPoint; import uk.me.parabola.mkgmap.general.MapLine; import uk.me.parabola.mkgmap.reader.osm.OsmConverter; import uk.me.parabola.mkgmap.reader.osm.Style; import uk.me.parabola.mkgmap.reader.osm.Way; import uk.me.parabola.mkgmap.osmstyle.StyleImpl; import uk.me.parabola.mkgmap.osmstyle.StyledConverter; import uk.me.parabola.mkgmap.osmstyle.eval.SyntaxException; import uk.me.parabola.util.EnhancedProperties; import static java.nio.channels.FileChannel.MapMode.READ_ONLY; /** * Create contour lines using an algorithm similar to that described in * <a href="http://mapcontext.com/autocarto/proceedings/auto-carto-5/pdf/an-adaptive-gr…">An Adaptive Grid Contouring Algorithm</a> by Downing and Zoraster. */ public abstract class DEM { final static double epsilon = 1e-9; final static double delta = 1.5; final static int maxPoints=200000; final static double minDist = 15; final static double maxDist = 21; final static double step = 0.01; final static double semiMajorAxis = 6378137.0; final static double inverseFlattening = 298.257223563; static int M=1200; static int N=M; static double res = 1.0/N; static int id = -1; short values[] = null; int lat; int lon; abstract double ele(int x, int y); abstract void read(int minLon, int minLat, int maxLon, int maxLat); abstract void serializeCopyRight(Writer out) throws IOException; public static void createContours(MapDetails mapper, EnhancedProperties config) { Area bounds = mapper.getBounds(); double minLat = Utils.toDegrees(bounds.getMinLat()); double minLon = Utils.toDegrees(bounds.getMinLong()); double maxLat = Utils.toDegrees(bounds.getMaxLat()); double maxLon = Utils.toDegrees(bounds.getMaxLong()); System.out.printf("bounds: %f %f %f %f\n", minLat, minLon, maxLat, maxLon); DEM data; String demType = config.getProperty("dem-type", "CGIAR"); String dataPath; if (demType.equals("ASTER")) { dataPath = config.getProperty("dem-path", "ASTER"); data = new ASTERGeoTiff(dataPath, minLat, minLon, maxLat, maxLon); } else if (demType.equals("CGIAR")) { dataPath = config.getProperty("dem-path", "CGIAR"); data = new CGIARGeoTiff(dataPath, minLat, minLon, maxLat, maxLon); } else { dataPath = config.getProperty("dem-path", "SRTM"); data = new HGTDEM(dataPath, minLat, minLon, maxLat, maxLon); } Isolines lines = data.new Isolines(data, minLat, minLon, maxLat, maxLon); int increment = config.getProperty("dem-increment", 10); double minHeight = lines.getMinHeight(); double maxHeight = lines.getMaxHeight(); int maxLevels = config.getProperty("dem-maxlevels", 100); while ((maxHeight-minHeight)/increment > maxLevels) increment *= 2; String loc = config.getProperty("style-file"); if (loc == null) loc = config.getProperty("map-features"); String name = config.getProperty("style"); if (loc == null && name == null) name = "default"; OsmConverter converter; try { Style style = new StyleImpl(loc, name); style.applyOptionOverride(config); converter = new StyledConverter(style, mapper, config); } catch (SyntaxException e) { System.err.println("Error in style: " + e.getMessage()); throw new ExitException("Could not open style " + name); } catch (FileNotFoundException e) { String name1 = (name != null)? name: loc; throw new ExitException("Could not open style " + name1); } for (int level=0; level<maxHeight; level+=increment) { if (level < minHeight) continue; // create isolines lines.addLevel(level); for (Isolines.Isoline line : lines.isolines) { Way way = new Way(id--, line.points); way.addTag("contour", "elevation"); way.addTag("ele", String.format("%d", (int) line.level)); converter.convertWay(way); } lines.isolines.clear(); } } public static class HGTDEM extends DEM { MappedByteBuffer buffer = null; public HGTDEM(String dataPath, double minLat, double minLon, double maxLat, double maxLon) { this.lat = (int) minLat; this.lon = (int) minLon; if (maxLat > lat+1 || maxLon > lon+1) throw new RuntimeException("Area too large (must not span more than one SRTM file)"); String northSouth = lat < 0 ? "S" : "N"; String eastWest = lon > 0 ? "E" : "W"; String fileName = String.format("%s/%s%02d%s%03d.hgt", dataPath, northSouth, lat < 0 ? -lat : lat, eastWest, lon < 0 ? -lon : lon); try { FileInputStream is = new FileInputStream(fileName); buffer = is.getChannel().map(READ_ONLY, 0, 2*(M+1)*(M+1)); } catch (Exception e) { throw new RuntimeException(e); } } public void read(int minLon, int minLat, int maxLon, int maxLat) { } public double ele(int x, int y) { return buffer.getShort(2*((M-y)*(M+1)+x))+delta; } public void serializeCopyRight(Writer out) throws IOException { out.write(" <copyright>\n"); out.write(" Contour lines generated from DEM data by NASA\n"); out.write(" </copyright>\n"); } } public static class CGIARGeoTiff extends DEM { Raster raster; String fileName; int minLat, minLon, maxLat, maxLon; PlanarImage image; public CGIARGeoTiff(String dataPath, double minLat, double minLon, double maxLat, double maxLon) { this.lat = ((int) (minLat/5))*5; this.lon = ((int) (minLon/5))*5; if (maxLat > lat+5 || maxLon > lon+5) throw new RuntimeException("Area too large (must not span more than one CGIAR GeoTIFF)"); int tileX, tileY; tileX = (180 + lon)/5 + 1; tileY = (60 - lat)/5; this.fileName = String.format("%s/srtm_%02d_%02d.tif", dataPath, tileX, tileY); init(); } public void serializeCopyRight(Writer out) throws IOException { out.write(" <copyright>\n"); out.write(" Contour lines generated from improved SRTM data by CIAT-CSI (see http://srtm.csi.cgiar.org)\n"); out.write(" </copyright>\n"); } public void read(int minLon, int minLat, int maxLon, int maxLat) { this.minLon = minLon; this.minLat = minLat; this.maxLon = maxLon; this.maxLat = maxLat; raster = image.getData(new java.awt.Rectangle(minLon, 6000-maxLat-1, maxLon-minLon+1, maxLat-minLat+1)); System.out.printf("read: %d %d %d %d\n", minLon, 6000-maxLat-1, maxLon-minLon+1, maxLat-minLat+1); } void init() { System.out.printf("CGIAR GeoTIFF: %s\n", fileName); N = 6000; M = 6000; res = 5.0/M; try { SeekableStream s = new FileSeekableStream(fileName); ParameterBlock pb = new ParameterBlock(); pb.add(s); TIFFDecodeParam param = new TIFFDecodeParam(); pb.add(param); RenderedOp op = JAI.create("tiff", pb); image = op.createInstance(); System.out.printf("Image: %d %d %d %d\n", image.getWidth(), image.getHeight(), image.getNumXTiles(), image.getNumYTiles()); } catch (Exception e) { throw new RuntimeException(e); } } public double ele(int x, int y) { try { int elevation = raster.getPixel(x, 6000-y-1, (int[])null)[0]; return elevation+delta; } catch (ArrayIndexOutOfBoundsException ex) { System.out.printf("ele: (%d, %d) (%d, %d, %d, %d) %s\n", x, 6000-y-1, raster.getMinX(), raster.getMinY(), raster.getWidth(), raster.getHeight(), ex.toString()); throw ex; } } } public static class ASTERGeoTiff extends DEM { Raster raster; String fileName; int minLat, minLon, maxLat, maxLon; PlanarImage image; public ASTERGeoTiff(String dataPath, double minLat, double minLon, double maxLat, double maxLon) { this.lat = (int) minLat; this.lon = (int) minLon; if (maxLat > lat+1 || maxLon > lon+1) throw new RuntimeException("Area too large (must not span more than one ASTER GeoTIFF)"); String northSouth = lat < 0 ? "S" : "N"; String eastWest = lon > 0 ? "E" : "W"; fileName = String.format("%s/ASTGTM_%s%02d%s%03d_dem.tif", dataPath, northSouth, lat < 0 ? -lat : lat, eastWest, lon < 0 ? -lon : lon); init(); } public void serializeCopyRight(Writer out) throws IOException { out.write(" <copyright>\n"); out.write(" Contour lines generated from DGM data by ASTER (see https://wist.echo.nasa.gov/~wist/api/imswelcome)\n"); out.write(" </copyright>\n"); } public void read(int minLon, int minLat, int maxLon, int maxLat) { this.minLon = minLon; this.minLat = minLat; this.maxLon = maxLon; this.maxLat = maxLat; raster = image.getData(new java.awt.Rectangle(minLon, 3601-maxLat-1, maxLon-minLon+1, maxLat-minLat+1)); System.out.printf("read: %d %d %d %d\n", minLon, 3601-maxLat-1, maxLon-minLon+1, maxLat-minLat+1); } void init() { System.out.printf("ASTER GeoTIFF: %s\n", fileName); N = 3600; M = 3600; res = 1.0/M; try { SeekableStream s = new FileSeekableStream(fileName); ParameterBlock pb = new ParameterBlock(); pb.add(s); TIFFDecodeParam param = new TIFFDecodeParam(); pb.add(param); RenderedOp op = JAI.create("tiff", pb); image = op.createInstance(); System.out.printf("Image: %d %d %d %d\n", image.getWidth(), image.getHeight(), image.getNumXTiles(), image.getNumYTiles()); } catch (Exception e) { throw new RuntimeException(e); } } public double ele(int x, int y) { try { int elevation = raster.getPixel(x, 3601-y-1, (int[])null)[0]; return elevation+delta; } catch (ArrayIndexOutOfBoundsException ex) { System.out.printf("ele: (%d, %d) (%d, %d, %d, %d) %s\n", x, 3601-y-1, raster.getMinX(), raster.getMinY(), raster.getWidth(), raster.getHeight(), ex.toString()); throw ex; } } } private static void debug(String format, Object ... args) { //System.out.printf(format + "\n", args); } int lastXi = -1; int lastYi = -1; final static int bcInv[][]= { { 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, { 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0}, { -3, 0, 0, 3, 0, 0, 0, 0, -2, 0, 0, -1, 0, 0, 0, 0}, { 2, 0, 0, -2, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0}, { 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0}, { 0, 0, 0, 0, -3, 0, 0, 3, 0, 0, 0, 0, -2, 0, 0, -1}, { 0, 0, 0, 0, 2, 0, 0, -2, 0, 0, 0, 0, 1, 0, 0, 1}, { -3, 3, 0, 0, -2, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, { 0, 0, 0, 0, 0, 0, 0, 0, -3, 3, 0, 0, -2, -1, 0, 0}, { 9, -9, 9, -9, 6, 3, -3, -6, 6, -6, -3, 3, 4, 2, 1, 2}, { -6, 6, -6, 6, -4, -2, 2, 4, -3, 3, 3, -3, -2, -1, -1, -2}, { 2, -2, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, { 0, 0, 0, 0, 0, 0, 0, 0, 2, -2, 0, 0, 1, 1, 0, 0}, { -6, 6, -6, 6, -3, -3, 3, 3, -4, 4, 2, -2, -2, -2, -1, -1}, { 4, -4, 4, -4, 2, 2, -2, -2, 2, -2, -2, 2, 1, 1, 1, 1} }; static int lastId=1000000000; static double lastX = 0; static double lastY = 0; static int edge[] = new int[2]; static int off0[][] = {{ 0, 0}, { 0, 0}, { 0, 1}, { 1, 1}}; static int off1[][] = {{ 1, 0}, { 0, 1}, { 1, 1}, { 1, 0}}; static int brd[] = { 1, 2, 4, 8 }; static int inv[] = { 4, 8, 1, 2 }; static int rev[] = { 2, 3, 0, 1 }; static int mov[][] = {{ 0, -1}, {-1, 0}, { 0, 1}, { 1, 0}}; double bc[][] = new double[4][4]; double bc_y[] = new double[4]; double bc_y1[] = new double[4]; double bc_y2[] = new double[4]; double bc_y12[] = new double[4]; double bc_Coeff[] = new double[16]; double bc_x[] = new double[16]; private void recalculateCoefficients(int xi, int yi) { double v00, vp0, v0p, vpp; double vm0, v0m, vpm, vmp, vmm, vmP, vPm; double vP0, v0P, vPp, vpP, vPP; v00 = ele(xi, yi); v0p = ele(xi, yi+1); vpp = ele(xi+1, yi+1); vp0 = ele(xi+1, yi); vm0 = ele(xi-1, yi); v0m = ele(xi, yi-1); vmp = ele(xi-1, yi+1); vpm = ele(xi+1, yi-1); vmm = ele(xi-1, yi-1); vmP = ele(xi+2, yi-1); vPm = ele(xi-1, yi+2); vP0 = ele(xi+2, yi); v0P = ele(xi, yi+2); vPp = ele(xi+2, yi+1); vpP = ele(xi+1, yi+2); vPP = ele(xi+2, yi+2); bc_y[0] = v00; bc_y[1] = vp0; bc_y[2] = vpp; bc_y[3] = v0p; bc_y1[0] = (vp0-vm0)/2; bc_y1[1] = (vP0-v00)/2; bc_y1[2] = (vPp-v0p)/2; bc_y1[3] = (vpp-vmp)/2; bc_y2[0] = (v0p-v0m)/2; bc_y2[1] = (vpp-vpm)/2; bc_y2[2] = (vpP-vp0)/2; bc_y2[3] = (v0P-v00)/2; bc_y12[0] = (vpp - vpm - vmp + vmm) / 4; bc_y12[0] = (vPp - vPm - v0p + v0m) / 4; bc_y12[2] = (vPP - vP0 - v0P + v00) / 4; bc_y12[0] = (vpP - vp0 - vmP + vm0) / 4; int j, i; double s; for (i=0; i<4; i++) { bc_x[i]=bc_y[i]; bc_x[i+4]=bc_y1[i]; bc_x[i+8]=bc_y2[i]; bc_x[i+12]=bc_y12[i]; } for (i=0; i<16; i++) { s = 0; for (int k=0; k<16; k++) s += bcInv[i][k]*bc_x[k]; bc_Coeff[i] = s; } int l = 0; for (i=0; i<4; i++) for (j=0; j<4; j++) bc[i][j] = bc_Coeff[l++]; } public double gradient(double lat, double lon, double[] grad) { grad[0] = 0; grad[1] = 0; double x = (lon-this.lon)/res; double y = (lat-this.lat)/res; int xi = (int) x; int yi = (int) y; if (lastXi != xi || lastYi != yi) { debug("new Cell for interpolation: %d %d", xi, yi); recalculateCoefficients(xi, yi); lastXi = xi; lastYi = yi; } double t = x - xi; double u = y - yi; if (xi < 0 || xi > N+1 || yi < 0 || yi > N+1) throw new IndexOutOfBoundsException(String.format("(%f, %f)->(%d, %d)", lat, lon, xi, yi)); double val = 0; for (int i=3; i>=0; i--) { val = t*val + ((bc[i][3]*u + bc[i][2])*u + bc[i][1])*u + bc[i][0]; grad[0] = u*grad[0] + (3*bc[3][i]*t + 2*bc[2][i])*t + bc[1][i]; grad[1] = t*grad[1] + (3*bc[i][3]*t + 2*bc[i][2])*t + bc[i][1]; } return val; } public double elevation(double lat, double lon) { double x = (lon-this.lon)/res; double y = (lat-this.lat)/res; int xi = (int) x; int yi = (int) y; if (lastXi != xi || lastYi != yi) { debug("new Cell for interpolation: %d %d", xi, yi); recalculateCoefficients(xi, yi); lastXi = xi; lastYi = yi; } double t = x - xi; double u = y - yi; if (xi < 0 || xi > N+1 || yi < 0 || yi > N+1) throw new IndexOutOfBoundsException(String.format("(%f, %f)->(%d, %d)", lat, lon, xi, yi)); double val = 0; for (int i=3; i>=0; i--) { val = t*val + ((bc[i][3]*u + bc[i][2])*u + bc[i][1])*u + bc[i][0]; } return val; } public double elevation(int x, int y) { if (x < 0 || x > N || y < 0 || y > N) throw new IndexOutOfBoundsException(String.format("elevation: %d %d", x, y)); return ele(x, y); } class Isolines { DEM data; int minX; int maxX; int minY; int maxY; double min; double max; Vector<Isoline> isolines = new Vector<Isoline>(); class Isoline { int id; Vector<Coord> points; double level; boolean isClosed; private Isoline(double level) { this.level = level; isClosed = false; id = lastId++; points = new Vector<Coord>(); } private class Edge implements Brent.Function { double x0, y0, x1, y1; Edge(double x0, double y0, double x1, double y1) { this.x0 = x0; this.y0 = y0; this.x1 = x1; this.y1 = y1; } public double eval(double d) { double f = data.elevation(x0 + d * (x1-x0), y0 + d * (y1-y0)) - level; //System.out.printf("evalEdge: %f %f\n", d, f); return f; } } private class FN implements Brent.Function { double x0, y0; double dx, dy; public void setParameter(double x0, double y0, double dx, double dy) { this.x0 = x0; this.y0 = y0; this.dx = dx; this.dy = dy; } public double eval(double t) { double f = data.elevation(y0+t*dy, x0+t*dx) - level; return f; } } private FN fn = new FN(); double grad[] = new double[2]; double px[] = new double[4]; double py[] = new double[4]; int edges[] = new int[4]; boolean addCell(Position p, int direction) { debug("addCell: %f %d %d %d %d", level, p.ix, p.iy, p.edge, direction); int c = 0; for (int k=0; k<4; k++) { if (k == p.edge) continue; int x0 = p.ix + off0[k][0]; int y0 = p.iy + off0[k][1]; int x1 = p.ix + off1[k][0]; int y1 = p.iy + off1[k][1]; double l0 = elevation(x0, y0) - level; double l1 = elevation(x1, y1) - level; if (Math.abs(l1) < epsilon || l0*l1 < 0) { edges[c] = k; Brent.Function f = new Edge(data.lat + y0 * DEM.res, data.lon + x0 * DEM.res, data.lat + y1 * DEM.res, data.lon + x1 * DEM.res); double f0 = elevation(x0, y0) - level; double f1 = elevation(x1, y1) - level; double delta; if (Math.abs(1) < epsilon) { delta = 1; } else if (Math.abs(f0) < epsilon) throw new RuntimeException("implementation error!"); else delta = Brent.zero(f, epsilon, 1-epsilon); px[c] = data.lon + (x0+delta*(x1-x0))*DEM.res; py[c] = data.lat + (y0+delta*(y1-y0))*DEM.res; c++; } } if (c == 1) { p.edge = edges[0]; double px0 = p.x; double py0 = p.y; p.x = px[0]; p.y = py[0]; double px1 = p.x; double py1 = p.y; double xMin = data.lon + p.ix * DEM.res; double xMax = xMin + DEM.res; double yMin = data.lat + p.iy * DEM.res; double yMax = yMin + DEM.res; refineAdaptively(xMin, yMin, xMax, yMax, px0, py0, px1, py1, direction, maxDist); addPoint(p.x, p.y, direction); p.moveCell(); return true; } else { debug("addCellByStepping: %d", c); return addCellByStepping(p, direction, c, edges, px, py); } } private void refineAdaptively(double xMin, double yMin, double xMax, double yMax, double x0, double y0, double x1, double y1, int direction, double maxDist) { double dist = orthodromicDistance(x0, y0, x1, y1); if (dist > maxDist) { double dx = x1-x0; double dy = y1-y0; double xm, ym, f0, f1, t0, t1, t; Brent.Function f; xm = x0 + 0.5*dx; ym = y0 + 0.5*dy; double n = Math.sqrt(dx*dx+dy*dy); fn.setParameter(xm, ym, -dy/n, dx/n); f = fn; t0 = -0.05*res; t1 = 0.05*res; f0 = f.eval(t0); f1 = f.eval(t1); int count = 0; while (f0 * f1 > 0 && count++ < 20) { if ((count & 1) > 0) t0 -= 0.05*res; else t1 += 0.05*res; f0 = f.eval(t0); f1 = f.eval(t1); debug("refine: %f %f %f %f", t0, t1, f0, f1); } if (f0 * f1 < 0) { t = Brent.zero(f, t0, t1); xm -= t*dy; ym += t*dx; } else { debug("refine failed: %f %f %f %f", t0, t1, f0, f1); if (false) throw new RuntimeException(String.format("refine failed: %f %f %f %f %f %f %f %f", xMin, yMin, xMax, yMax, x0, y0, x1, y1)); return; } if (xm > xMin && xm < xMax && ym > yMin && ym < yMax) refineAdaptively(xMin, yMin, xMax, yMax, x0, y0, xm, ym, direction, maxDist*1.1); addPoint(xm, ym, direction); if (xm > xMin && xm < xMax && ym > yMin && ym < yMax) refineAdaptively(xMin, yMin, xMax, yMax, xm, ym, x1, y1, direction, maxDist*1.1); } } boolean addCellByStepping(Position p, int direction, int numEdges, int[] edges, double[] px, double[] py) { debug("addCellByStepping: %f %d %d %d %d", level, p.ix, p.iy, p.edge, direction); double xMin = data.lon + p.ix * DEM.res; double xMax = xMin + DEM.res; double yMin = data.lat + p.iy * DEM.res; double yMax = yMin + DEM.res; double dt, t0, t1; double f0, f1; boolean edgeHit = false; double h[] = new double[4]; int dir; double n2 = Math.sqrt(1.0/(grad[0]*grad[0] + grad[1]*grad[1])); double dx; double dy; int count=0; int iMin = -1; double md = 5000; for (int i=0; i<numEdges; i++) { gradient(p.y, p.x, grad); double dist = orthodromicDistance(p.x, p.y, px[i], py[i]); //double dist = (px[i]-p.x)*(px[i]-p.x)*grad[0]*grad[0] + (py[i]-p.y)*(py[i]-p.y)*grad[1]*grad[1]; debug("distance %d: %f", i, dist); if (dist < md && (visited[p.iy*(N+1)+p.ix] & brd[edges[i]]) == 0) { md = dist; iMin = i; } } p.edge = edges[iMin]; double px0 = p.x; double py0 = p.y; p.x = px[iMin]; p.y = py[iMin]; double px1 = p.x; double py1 = p.y; xMin = data.lon + p.ix * DEM.res; xMax = xMin + DEM.res; yMin = data.lat + p.iy * DEM.res; yMax = yMin + DEM.res; refineAdaptively(xMin, yMin, xMax, yMax, px0, py0, px1, py1, direction, maxDist); addPoint(p.x, p.y, direction); p.moveCell(); return true; } private void addPoint(double x, double y, int direction) { double dist = orthodromicDistance(x, y, lastX, lastY); debug("addPoint: %f %f %f", x, y, dist); if (dist > minDist) { if (direction > 0) points.add(0, new Coord(y, x)); else points.add(points.size(), new Coord(y, x)); lastX = x; lastY = y; } } private void close() { points.add(points.size(), points.get(0)); isClosed = true; } private void addMove(int x0, int y0, int x1, int y1, int direction) { double x; double y; if (x0 < minX || x0 >= maxX || y0 < minY || y0 >= maxY) return; double l0 = data.elevation(x0, y0); double l1 = data.elevation(x1, y1); debug("addMove: %d %d %d %d %.2f %.2f %d", x0, y0, x1, y1, l0, l1, direction); if ((l0 < level && l1 < level) || (l0 > level && l1 > level)) throw new RuntimeException("implementation error"); if (l0 == level) { x = data.lon + x0 * DEM.res; y = data.lat + y0 * DEM.res; } else { double delta = (l0-level)/(l0-l1); x = data.lon + (x0 + delta * (x1-x0)) * DEM.res; y = data.lat + (y0 + delta * (y1-y0)) * DEM.res; } double dist = orthodromicDistance(x, y, lastX, lastY); debug("levels: %d %d %f, %d %d %f: %f %f %f", x0, y0, l0, x1, y1, l1, x, y, dist); if (dist > 1) { if (direction > 0) points.add(0, new Coord(y, x)); else points.add(points.size(), new Coord(y, x)); lastX = x; lastY = y; } } } public Isolines(DEM data, int minX, int maxX, int minY, int maxY) { this.data = data; this.minX = minX; this.maxX = maxX; this.minY = minY; this.maxY = maxY; init(); } public Isolines(DEM data, double minLat, double minLon, double maxLat, double maxLon) { System.out.printf("init: %f %f %f %f\n", minLat, minLon, maxLat, maxLon); this.data = data; this.minX = (int) ((minLon-data.lon)/data.res); this.minY = (int) ((minLat-data.lat)/data.res); this.maxX = (int) ((maxLon-data.lon)/data.res); this.maxY = (int) ((maxLat-data.lat)/data.res); init(); } private void init() { System.out.printf("init: %d %d %d %d\n", minX, minY, maxX, maxY); data.read(minX-2, minY-2, maxX+2, maxY+2); // we need some overlap for bicubic interpolation max = -1000; min = 10000; for (int i=minX; i<maxX; i++) for (int j=minY; j<maxY; j++) { if (data.elevation(i, j) < min) min = data.elevation(i, j); if (data.elevation(i, j) > max) max = data.elevation(i, j); } debug("min: %f, max: %f\n", min, max); } double getMinHeight() { return min; } double getMaxHeight() { return max; } public void addLevels(int start, int increment) { for (int level=start; level<max; level+=increment) addLevel(level); } private class Edge implements Brent.Function { double x0, y0, x1, y1, level; Edge(double x0, double y0, double x1, double y1, double level) { this.x0 = x0; this.y0 = y0; this.x1 = x1; this.y1 = y1; this.level = level; } public double eval(double d) { double f = data.elevation(x0 + d * (x1-x0), y0 + d * (y1-y0)) - level; //System.out.printf("evalEdge: %f %f\n", d, f); return f; } } private class Position { int ix, iy; double x, y; int edge; Position(int ix, int iy, double x, double y, int edge) { this.ix = ix; this.iy = iy; this.x = x; this.y = y; this.edge = edge; } Position(Position p) { this.ix = p.ix; this.iy = p.iy; this.x = p.x; this.y = p.y; this.edge = p.edge; } void markEdge() { debug("marking edge: %d %d %d %d", ix, iy, edge, brd[edge]); visited[iy*(N+1)+ix] |= brd[edge]; } void moveCell() { markEdge(); ix += mov[edge][0]; iy += mov[edge][1]; edge = rev[edge]; markEdge(); } } byte visited[] = new byte[(N+1)*(N+1)]; public void addLevel(double level) { if (level < min || level > max) return; System.out.printf("addLevel: %f\n", level); java.util.Arrays.fill(visited, (byte) 0); for (int y=minY; y<maxY; y++) { for (int x=minX; x<maxX; x++) { byte v = 0; // Mark the borders of the cell, represented by the four points (i, j), (i+1, j), (i, j+1), (i+1, j+1), // which are intersected by the contour. The values are: // 1: top // 2: left // 4: bottom // 8: right if (data.elevation(x, y) >= level) { if (data.elevation(x+1, y) < level) { v |= 1; } if (data.elevation(x, y+1) < level) { v |= 2; } } else { if (data.elevation(x+1, y) > level) { v |= 1; } if (data.elevation(x, y+1) > level) { v |= 2; } } int k=-1; if ((v&1) > 0 && (visited[y*(N+1)+x]&1) == 0) { k=0; } else if ((v&2) > 0 && (visited[y*(N+1)+x]&2) == 0) { k=1; } if (k>=0) { int x0 = x + off0[k][0]; int y0 = y + off0[k][1]; int x1 = x + off1[k][0]; int y1 = y + off1[k][1]; try { Brent.Function f = new Edge(data.lat + y0 * DEM.res, data.lon + x0 * DEM.res, data.lat + y1 * DEM.res, data.lon + x1 * DEM.res, level); double f0 = elevation(x0, y0) - level; double f1 = elevation(x1, y1) - level; double delta; if (Math.abs(f0) < epsilon) { delta = 0; } else if (Math.abs(f1) < epsilon) continue; else delta = Brent.zero(f, 0, 1-epsilon); Position p = new Position(x, y, data.lon + (x0+delta*(x1-x0))*DEM.res, data.lat + (y0+delta*(y1-y0))*DEM.res, k); p.markEdge(); isolines.add(traceByStepping(level, p)); } catch (RuntimeException ex) { debug("error: %s", ex.toString()); ex.printStackTrace(); return; } } } } } private Isoline traceByStepping(double level, Position p) { debug("traceByStepping: starting contour %f %d %d %f %f %d", level, p.ix, p.iy, p.x, p.y, p.edge); int direction = 1; int n = 0; Position startP = new Position(p); boolean foundEnd = false; Isoline line = new Isoline(level); while (true) { debug("traceByStepping: %f %d %d %f %f %d", level, p.ix, p.iy, p.x, p.y, p.edge); visited[p.iy*(N+1)+p.ix] |= brd[p.edge]; if (n>0 && p.ix == startP.ix && p.iy == startP.iy && orthodromicDistance(p.x, p.y, startP.x, startP.y) < 5) { debug("closed curve!"); line.close(); break; } else if (p.ix < minX || p.iy < minY || p.ix >= maxX || p.iy >= maxY) { if (foundEnd) // did we already reach one end? { debug("second border reached!"); break; } else { debug("border reached!"); foundEnd = true; n = 0; direction *= -1; p = new Position(startP); p.moveCell(); continue; } } n++; if (!line.addCell(p, direction) || line.points.size() > maxPoints) { debug("ending contour"); isolines.add(line); return line; } } return line; } } /** * Returns the orthodromic distance between two geographic coordinates in WGS84 datum. * The orthodromic distance is the shortest distance between two points * on a sphere's surface. The orthodromic path is always on a great circle. * * @param x1 Longitude of first point (in degrees). * @param y1 Latitude of first point (in degrees). * @param x2 Longitude of second point (in degrees). * @param y2 Latitude of second point (in degrees). * @return The orthodromic distance (in meters). * */ public static double orthodromicDistance(double x1, double y1, double x2, double y2) { x1 = Math.toRadians(x1); y1 = Math.toRadians(y1); x2 = Math.toRadians(x2); y2 = Math.toRadians(y2); final int MAX_ITERATIONS = 100; final double EPS = 5.0E-14; final double F = 1 / inverseFlattening; final double R = 1 - F; double tu1 = R * Math.sin(y1) / Math.cos(y1); double tu2 = R * Math.sin(y2) / Math.cos(y2); double cu1 = 1 / Math.sqrt(tu1 * tu1 + 1); double cu2 = 1 / Math.sqrt(tu2 * tu2 + 1); double su1 = cu1 * tu1; double s = cu1 * cu2; double baz = s * tu2; double faz = baz * tu1; double x = x2 - x1; for (int i = 0; i < MAX_ITERATIONS; i++) { final double sx = Math.sin(x); final double cx = Math.cos(x); tu1 = cu2 * sx; tu2 = baz - su1 * cu2 * cx; final double sy = Math.sqrt(tu1*tu1 + tu2*tu2);; final double cy = s * cx + faz; final double y = Math.atan2(sy, cy); final double SA = s * sx / sy; final double c2a = 1 - SA * SA; double cz = faz + faz; if (c2a > 0) { cz = -cz / c2a + cy; } double e = cz * cz * 2 - 1; double c = ((-3 * c2a + 4) * F + 4) * c2a * F / 16; double d = x; x = ((e * cy * c + cz) * sy * c + y) * SA; x = (1 - c) * x * F + x2 - x1; if (Math.abs(d - x) <= EPS) { x = Math.sqrt((1 / (R * R) - 1) * c2a + 1) + 1; x = (x - 2) / x; c = 1 - x; c = (x * x / 4 + 1) / c; d = (0.375 * x * x - 1) * x; x = e * cy; s = 1 - 2 * e; s = ((((sy * sy * 4 - 3) * s * cz * d / 6 - x) * d / 4 + cz) * sy * d + y) * c * R * semiMajorAxis; return s; } } final double LEPS = 1.0E-10; if (Math.abs(x1 - x2) <= LEPS && Math.abs(y1 - y2) <= LEPS) { return 0; } if (Math.abs(y1) <= LEPS && Math.abs(y2) <= LEPS) { return Math.abs(x1 - x2) * semiMajorAxis; } throw new ArithmeticException("no convergence"); } } /* * Copyright (C) 2009 Christian Gawron * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License version 2 as * published by the Free Software Foundation. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * * Author: Christian Gawron * Create date: 03-Jul-2009 */ package uk.me.parabola.mkgmap.reader.dem; /** * Find zero of a function using Brent's method. */ public class Brent { public interface Function { public double eval(double x); } static double epsilon = 3.0e-10; public static double zero(Function f, double x1, double x2) { return zero(f, x1, x2, 1e-8, 100); } public static double zero(Function f, double x1, double x2, double tol, int maxit) { double a=x1, b=x2, c=x2, d=0, e=0, min1, min2; double fa=f.eval(a), fb=f.eval(b), fc; double p, q, r, s, tolerance, xm; if ((fa > 0.0 && fb > 0.0) || (fa < 0.0 && fb < 0.0)) throw new ArithmeticException("Root must be bracketed"); fc=fb; for (int iterations=0; iterations < maxit; iterations++) { if ((fb > 0.0 && fc > 0.0) || (fb < 0.0 && fc < 0.0)) { c=a; fc=fa; e=d=b-a; } if (Math.abs(fc) < Math.abs(fb)) { a=b; b=c; c=a; fa=fb; fb=fc; fc=fa; } tolerance=2.0*epsilon*Math.abs(b)+0.5*tol; xm=0.5*(c-b); if (Math.abs(xm) <= tolerance || fb == 0.0) return b; if (Math.abs(e) >= tolerance && Math.abs(fa) > Math.abs(fb)) { s=fb/fa; if (a == c) { p=2.0*xm*s; q=1.0-s; } else { q=fa/fc; r=fb/fc; p=s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0)); q=(q-1.0)*(r-1.0)*(s-1.0); } if (p > 0.0) q = -q; p=Math.abs(p); min1=3.0*xm*q-Math.abs(tolerance*q); min2=Math.abs(e*q); if (2.0*p < (min1 < min2 ? min1 : min2)) { e=d; d=p/q; } else { d=xm; e=d; } } else { d=xm; e=d; } a=b; fa=fb; if (Math.abs(d) > tolerance) b += d; else b += xm >= 0 ? tolerance : -tolerance; fb=f.eval(b); } throw new ArithmeticException("Maximum number of iterations exceeded"); } } Index: src/uk/me/parabola/mkgmap/Option.java =================================================================== --- src/uk/me/parabola/mkgmap/Option.java (Revision 1078) +++ src/uk/me/parabola/mkgmap/Option.java (Arbeitskopie) @@ -34,7 +34,7 @@ } } - protected Option(String option, String value) { + public Option(String option, String value) { this.option = option; this.value = value; } Index: src/uk/me/parabola/mkgmap/CommandArgsReader.java =================================================================== --- src/uk/me/parabola/mkgmap/CommandArgsReader.java (Revision 1078) +++ src/uk/me/parabola/mkgmap/CommandArgsReader.java (Arbeitskopie) @@ -129,7 +129,7 @@ * @param option The option name. * @param value Its value. */ - private void addOption(String option, String value) { + public void addOption(String option, String value) { CommandOption opt = new CommandOption(option, value); addOption(opt); } @@ -181,13 +181,21 @@ return arglist.iterator(); } + public ArgList getArgList() { + return arglist; + } + + public EnhancedProperties getArgs() { + return args; + } + /** * Read a config file that contains more options. When the number of * options becomes large it is more convenient to place them in a file. * * @param filename The filename to obtain options from. */ - private void readConfigFile(String filename) { + public void readConfigFile(String filename) { Options opts = new Options(new OptionProcessor() { public void processOption(Option opt) { log.debug("incoming opt", opt.getOption(), opt.getValue()); @@ -206,14 +214,14 @@ * the argument to be processed in order. Options can be intersperced with * filenames. The options take effect where they appear. */ - interface ArgType { + public interface ArgType { public abstract void processArg(); } /** * A filename. */ - class Filename implements ArgType { + public class Filename implements ArgType { private final String name; private boolean useFilenameAsMapname = true; @@ -270,7 +278,7 @@ /** * An option argument. A key value pair. */ - class CommandOption implements ArgType { + public class CommandOption implements ArgType { private final Option option; private CommandOption(Option option) { @@ -297,7 +305,7 @@ /** * The arguments are held in this list. */ - class ArgList implements Iterable<CommandArgsReader.ArgType> { + public class ArgList implements Iterable<CommandArgsReader.ArgType> { private final List<ArgType> alist; private int filenameCount; Index: src/uk/me/parabola/mkgmap/reader/MapperBasedMapDataSource.java =================================================================== --- src/uk/me/parabola/mkgmap/reader/MapperBasedMapDataSource.java (Revision 1078) +++ src/uk/me/parabola/mkgmap/reader/MapperBasedMapDataSource.java (Arbeitskopie) @@ -33,6 +33,7 @@ import uk.me.parabola.mkgmap.general.MapPoint; import uk.me.parabola.mkgmap.general.MapShape; import uk.me.parabola.mkgmap.general.RoadNetwork; +import uk.me.parabola.mkgmap.reader.dem.DEM; import uk.me.parabola.util.Configurable; import uk.me.parabola.util.EnhancedProperties; @@ -146,5 +147,8 @@ // the overview section. mapper.addShape(background); } + if (getConfig().getProperty("contours", false)) { + DEM.createContours(mapper, getConfig()); + } } } Index: src/uk/me/parabola/mkgmap/reader/osm/Way.java =================================================================== --- src/uk/me/parabola/mkgmap/reader/osm/Way.java (Revision 1078) +++ src/uk/me/parabola/mkgmap/reader/osm/Way.java (Arbeitskopie) @@ -30,12 +30,18 @@ */ public class Way extends Element { - private final List<Coord> points = new ArrayList<Coord>(); + private final List<Coord> points; public Way(long id) { + points = new ArrayList<Coord>(); setId(id); } + public Way(long id, List<Coord> points) { + this.points = points; + setId(id); + } + /** * Get the points that make up the way. We attempt to re-order the segments * and return a list of points that traces the route of the way.

2 1

Route recalculation; error circle
by Marko Mäkelä 04 Jul '09

04 Jul '09

I made some observations on my Edge 705 yesterday. The route recalculation sensitivity varies in bicycle mode. Sometimes, I can ride a couple hundred metres until getting a notice; other times, I get the notice within a couple dozen metres. This might have something to do with junctions and parallel roads. Once, I rode on a cycleway that is almost parallel to the calculated route. In the next junction, it told me to turn the right way, even though the current position was indicated correctly. Another observation: I can confirm what someone said earlier about the error circle around the current position marker. Even though the GPS accuracy is something like 5 metres most of the time, the radius of the error circle on the map is something like 20 metres. It will get bigger when there are obstructions, such as tunnels or tall buildings. Marko

2 1

Commit: r1078: Fix bugs in date/time encoding in headers.
by svn commit 04 Jul '09

04 Jul '09

Version 1078 was commited by markb on 2009-07-04 01:43:14 +0100 (Sat, 04 Jul 2009) Fix bugs in date/time encoding in headers. Hour is now modulo 24 and not modulo 12. In one place, year had 1900 erroneously subtracted from it.

1 0

Patch so that the --add-pois-to-areas option applies actions defined in the points style file and takes the POI name from there
by Thilo Hannemann 03 Jul '09

03 Jul '09

This one is a little complicated: When using the --add-pois-to-areas option, for every shape whose tags will match an entry in the points style file, a POI entry in the map will be created in the center of the shape (if there is not already one POI with the same name inside the shape). The POI will have the name of the *shape* and no actions from the points style file will be applied. I think that is not how it should be. One of the most important cases for this option is when you have an POI that is in the OSM data as a building (close way with "building=yes") carrying the tags for the POI. You might now have lots of name-tweaking rules in you points style file to generate nice looking names for hotels, etc. But in the polygons style file there is a single entry that handles all tags with "building=yes". What you now get is that all the name-tweaking rules *are not applied* if the POI comes in the form of a shape, only if it comes as a single node. To get around this you had to duplicate your name-tweaking for each POI in the polygons style file. This will need lots of rules that from the perspective of the polygons style file all will do the same: generate a shape for a building. So this patch will modify the behaviour so that a POI that is generated from a shape will have all the actions defined in the points style file applied to it and that it will also get the name defined by the points style file. Regards Thilo

1 0

Commit: r1077: Polish input reader now recognises codes like ~[0x05].
by svn commit 03 Jul '09

03 Jul '09

Version 1077 was commited by steve on 2009-07-03 14:42:32 +0100 (Fri, 03 Jul 2009) Polish input reader now recognises codes like ~[0x05].

1 0

Recommended "stable" version?
by Nop 03 Jul '09

03 Jul '09

Hi! According to the Wiki page, the last version recommended as reasonably stable is 858. As I needed some particular improvements, I am using 1001, which works reasonably well. However, there are some weird glitches in my map that may be caused by problems in my preprocessing or in mkgmap with equal propability. So an update to a later version might in order but quite a number of the posts here seem to indicate that many things are very experimental so this might be a change for the worse at the moment. My map is not intended for routing, so will not profit from any features to that end, I am rather hoping for general improvements of processing and reliability. So the question is: Are there any plans to work towards another recommended reasonably stable version? Or do any of you have a personal recommendation for a build that in your experience should not cause trouble? thanks Nop _______________________________________________

2 2

Commit: r1076: Make highway shield work with 6-bit labels.
by svn commit 03 Jul '09

03 Jul '09

Version 1076 was commited by steve on 2009-07-03 12:26:02 +0100 (Fri, 03 Jul 2009) Make highway shield work with 6-bit labels. Add a test for the shield substitution.

1 0

Commit: r1073: Add highway shields and other ways of modifying values in the style files.
by svn commit 03 Jul '09

03 Jul '09

Version 1073 was commited by steve on 2009-07-02 12:55:28 +0100 (Thu, 02 Jul 2009) Add highway shields and other ways of modifying values in the style files. Implemented as a series of modifiers that can be used in the style files. HighwaySymbolFilter: inserts the highway symbol markers after cleaning up the name. HeightFilter: Prepends the value that makes the garmin display as an elevation. SubstitutionFilter: allows arbitary substitutions within a value. PrependFilter: Prepend one the garmin magic values to a string. Patch by Toby Speight as the end result of a discussion involving a number of people see http://www.mkgmap.org.uk/pipermail/mkgmap-dev/2009q2/001573.html for the full details. This only appears to work with 8bit outputs.

5 7