[PATCH] Creating contour lines from DEM data
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-gri...) 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-grid-contouring-algorithm.pdf">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"); } }
Hi Christian, congratulations, that is very good work. As there are some people on this list who will publish their created maps I'll have to add a caveat though: You should make sure that the data source you use for the generation of the contour lines is compatible with the OSM license. See http://wiki.openstreetmap.org/wiki/Legal_FAQ and http://wiki.openstreetmap.org/wiki/Common_licence_interpretations. The reason why I bring that up is that for example the improved SRTM data by CIAT-CSI is in my opinion *not* compatible with the OSM license! Regards Thilo
I have some cycle networks that share the same ways. This means that a road belongs to more than one cycle network. In this case I want to see the ref of the cycle networks in the name of the road, e.g. "<road_name> (RRM, VBT, MWW)". My file 'relations' contains this line: type=route & route=bicycle { apply { set route_name='${ref}' }} In this case the tag route_name of the way is set for every route and the last one wins. So I tried: type=route & route=bicycle & route_name!=* { apply { set route_name='${ref}' }} type=route & route=bicycle & route_name=* { apply { set route_name='${route_name}, ${ref}' }} This doesn't work because we have no access to the tags of the way. Who has an idea how I can join the different names? Rudi
I want to show road signs for restrictions. So I tried the following: relations: type=restriction { apply { set relation_restriction='${restriction}' }} points: relation_restriction=only_straight_on [0x300f level 0] This creates a sign at the crossing. But I'd prefer to see the sign at the 'from' part of the relation. So I added the following line to the file polygons: relation_restriction=only_straight_on [0x21 level 0] In this case I get 3 signs for every restriction (from, via, to). (This is a 'hack' to get a POI for a way with --add-pois-to-areas. The polygon is not displayed because the way consists of 2 nodes). How can I restrict the POI to the 'from' member of the relation? The following line doesn't work because 'role' doesn't belong to the relation, it belongs to the members of the relation: type=restriction & role=from { apply { set relation_restriction='${restriction}' }} Rudi
Hi! Thilo Hannemann schrieb:
As there are some people on this list who will publish their created maps I'll have to add a caveat though: You should make sure that the data source you use for the generation of the contour lines is compatible with the OSM license. See http://wiki.openstreetmap.org/wiki/Legal_FAQ and http://wiki.openstreetmap.org/wiki/Common_licence_interpretations.
The reason why I bring that up is that for example the improved SRTM data by CIAT-CSI is in my opinion *not* compatible with the OSM license!
A very good hint. The understanding of all discussions on the matter is that the current OSM licence does not allow mixing with CIAT data nor with ASTER data. So you cannot publish those maps. bye Nop
A very good hint. The understanding of all discussions on the matter is that the current OSM licence does not allow mixing with CIAT data nor with ASTER data. So you cannot publish those maps. Could you elaborate on this? I would understand that e.g. the license of the ASTER data forbids publishing of derived work (although I haven't checked this). But is there really a point in the OSM license which forbids creating derived works (maps) which as an "add-on" contain features created from proprietary data sources?
Best wishes Christian
Am 04.07.2009 um 08:04 schrieb Christian Gawron:
A very good hint. The understanding of all discussions on the matter is that the current OSM licence does not allow mixing with CIAT data nor with ASTER data. So you cannot publish those maps. Could you elaborate on this? I would understand that e.g. the license of the ASTER data forbids publishing of derived work (although I haven't checked this). But is there really a point in the OSM license which forbids creating derived works (maps) which as an "add-on" contain features created from proprietary data sources?
Yes it is. See the section "I created a layer on top of an OSM map. What do I have to put under your license?" from http://wiki.openstreetmap.org/wiki/Common_licence_interpretations : -- You have to determine whether what you have created is a Collective Work or a Derivative work, under the terms of the OSM licence. • If what you create is based on OSM data (for example if you create a new layer by looking at the OSM data and refering to locations on it) then it is likely you have created a derivative work. • If you generate a merged work with OSM data and other data (such as a printed map or pdf map) where the non-OSM data can no longer be considered to be separate and independent from the OSM data, is is likely you have created a derivative work. • If you overlay OSM data with your own data created from other sources (for example you going out there with a GPS receiver) and the layers are kept separate and independent, and the OSM layer is unchanged, then you may have created a collective work. If you have created a derivative work, the work as a whole must be subject to the OSM licence. If you have created a collective work, then only the OSM component of the work must be subject to the OSM licence. -- So especially with your solution to integrate the contour lines into the map this is IMHO a derivative work, so that all components of it must be subject to the OSM license. Regards Thilo
0> In article <68B4C1B8-4426-4453-8B09-00B90F870F01@gmx.de>, 0> Thilo Hannemann <URL:mailto:thannema@gmx.de> ("Thilo") wrote: Thilo> So especially with your solution to integrate the contour lines Thilo> into the map this is IMHO a derivative work, so that all components Thilo> of it must be subject to the OSM license. Agreed. I believe (but I'm not any kind of lawyer!). I'd say that separate invocations of mkgmap (one with OSM data and one with contour data) would be fine - it's perfectly possible for a user to then combine them into a single gmapsupp for their device. In fact, separate maps is probably also the best plan for those of us generating maps for our own use (and not distributing them), as the contour maps should be a one-shot task. Things would get complicated if the contour generation algorithm uses water body data to discover flat regions, though. P.S. it would be interesting to see contour regions (as opposed to lines) suitable for elevation shading layers; another idea is to compute slope aspect for a hill shading layer - maybe I'll look into that when I get time...
Am 04.07.2009 um 10:04 schrieb Toby Speight:
P.S. it would be interesting to see contour regions (as opposed to lines) suitable for elevation shading layers; another idea is to compute slope aspect for a hill shading layer - maybe I'll look into that when I get time...
Actually I've done that already using octave (matrix calculation application, among other uses). I can share the scripts if you like. The problem there is not to generate the data, but to display it in a usable way on the GPS unit. The DEM maps with elevation and hill shading look "Wow" before you try to convert them into an .IMG file. Then they look like crap, because you have nowhere near enough different layers available. There is an article about how to create a bigger number of shadings by dithering: http://www.cgpsmapper.com/download/Creating%20custom%20types%20to%20represen... But even with this method I think the resulting maps look crap, especially if you look at the maps that the GPS units shows if a DEM model for the map is available. So, anybody interested in figuring out the Garmin DEM file format? I would assume that it might not be too complicated if the DEM data is stored in a sampled format like the SRMT data (hopefully not compressed). This could have the added benefit that you might also be able to view the elevation data and use the 3D display mode on the GPS unit. Regards Thilo
Thilo Hannemann schrieb:
-- You have to determine whether what you have created is a Collective Work or a Derivative work, under the terms of the OSM licence. • If what you create is based on OSM data (for example if you create a new layer by looking at the OSM data and refering to locations on it) then it is likely you have created a derivative work. • If you generate a merged work with OSM data and other data (such as a printed map or pdf map) where the non-OSM data can no longer be considered to be separate and independent from the OSM data, is is likely you have created a derivative work. • If you overlay OSM data with your own data created from other sources (for example you going out there with a GPS receiver) and the layers are kept separate and independent, and the OSM layer is unchanged, then you may have created a collective work. If you have created a derivative work, the work as a whole must be subject to the OSM licence. If you have created a collective work, then only the OSM component of the work must be subject to the OSM licence. --
So especially with your solution to integrate the contour lines into the map this is IMHO a derivative work, so that all components of it must be subject to the OSM license. I don't see any difference between merging externally generated contour lines into one IMG file or generating the contour lines with mkgmap - the resulting IMG file looks the same. And I don't agree that this would be a derived work. OSM does not contain any DEM data or contour lines, so by adding them the OSM layer remains unchanged and clearly discernible from the contour lines, so I would say the result is a collective work. But of course anyone planning to distribute such maps should check this more thoroughly (maybe the guys who are hosting the cycle layer or the "Wanderkarte" did this).
Best wishes Christian
Hi! Christian Gawron schrieb:
I don't see any difference between merging externally generated contour lines into one IMG file or generating the contour lines with mkgmap - the resulting IMG file looks the same.
The difference is what is being published: You may not create a map from both sources and redistribute it. But you may distribute each layer seperately and let the end user combine them into a map, which he may not redistribute.
And I don't agree that this would be a derived work. OSM does not contain any DEM data or contour lines, so by adding them the OSM layer remains unchanged and clearly discernible from the contour lines, so I would say the result is a collective work.
Discernible is not enough. You need to read through the full legalese. A collective work is defined as a collection of unmodified entities. This is definitely not true for a gmapsupp.img. Took me a while to understand and believe how bad the current OSM licence is, too. bye Nop
Nop schrieb:
And I don't agree that this would be a derived work. OSM does not contain any DEM data or contour lines, so by adding them the OSM layer remains unchanged and clearly discernible from the contour lines, so I would say the result is a collective work.
Discernible is not enough. You need to read through the full legalese. A collective work is defined as a collection of unmodified entities. This is definitely not true for a gmapsupp.img.
Took me a while to understand and believe how bad the current OSM licence is, too. Thank you for the clarification. So maybe I should add an option to put the contours into a separate IMG file. Right now you could achieve this by using an input file which contains no data but the bounding box.
I think using one IMG has the advantage that using a TYP file you can specify the drawing order (first the "background" polygons, than the contours, then the streets). Best wishes Christian
Am 04.07.2009 um 13:06 schrieb Christian Gawron:
Nop schrieb:
And I don't agree that this would be a derived work. OSM does not contain any DEM data or contour lines, so by adding them the OSM layer remains unchanged and clearly discernible from the contour lines, so I would say the result is a collective work.
Discernible is not enough. You need to read through the full legalese. A collective work is defined as a collection of unmodified entities. This is definitely not true for a gmapsupp.img.
Took me a while to understand and believe how bad the current OSM licence is, too. Thank you for the clarification. So maybe I should add an option to put the contours into a separate IMG file. Right now you could achieve this by using an input file which contains no data but the bounding box.
I think using one IMG has the advantage that using a TYP file you can specify the drawing order (first the "background" polygons, than the contours, then the streets).
Actually, using different IMG files will allow you to do that, but with one IMG file for everything it will be difficult. You cannot specify the drawing order of ways in the TYP file, only of polygons. And as I understand it, it is not easy to control the drawing order of ways in any other way. Nop seems to be able to do it by controlling the order the ways are written into the IMG file, but I was largely unsucessfull using that method. Regards Thilo
Thilo Hannemann schrieb:
I think using one IMG has the advantage that using a TYP file you can specify the drawing order (first the "background" polygons, than the contours, then the streets).
Actually, using different IMG files will allow you to do that, but with one IMG file for everything it will be difficult. You cannot specify the drawing order of ways in the TYP file, only of polygons. And as I understand it, it is not easy to control the drawing order of ways in any other way. Nop seems to be able to do it by controlling the order the ways are written into the IMG file, but I was largely unsucessfull using that method. But if you put the IMG with the contours on top, wouldn't the contour lines be on top of the ways? And if you do it the other way around, wouldn't the background polygons (eg. forrest) hide the contour lines? How did you solve this?
I actually tried using separate IMG files first, but I didn't get "nice" results, so I used XInclude to include the contour lines from the OSM file (before I wrote the code to create the contour lines in mkgmap). Best wishes Christian
What I actually do is to include the contour lines into the IMG file, because that works best for me. If you put them into a separate IMG file and then the polygons into a separate IMG file and also the ways into their separate IMG file you can create layers and control the draw order - in a way. We just lately had a discussion about this on the mailing list, see the thread 'Option "preserve-element-order"'. Regards Thilo Am 04.07.2009 um 22:42 schrieb Christian Gawron:
Thilo Hannemann schrieb:
I think using one IMG has the advantage that using a TYP file you can specify the drawing order (first the "background" polygons, than the contours, then the streets).
Actually, using different IMG files will allow you to do that, but with one IMG file for everything it will be difficult. You cannot specify the drawing order of ways in the TYP file, only of polygons. And as I understand it, it is not easy to control the drawing order of ways in any other way. Nop seems to be able to do it by controlling the order the ways are written into the IMG file, but I was largely unsucessfull using that method. But if you put the IMG with the contours on top, wouldn't the contour lines be on top of the ways? And if you do it the other way around, wouldn't the background polygons (eg. forrest) hide the contour lines? How did you solve this?
I actually tried using separate IMG files first, but I didn't get "nice" results, so I used XInclude to include the contour lines from the OSM file (before I wrote the code to create the contour lines in mkgmap).
Best wishes Christian _______________________________________________ mkgmap-dev mailing list mkgmap-dev@lists.mkgmap.org.uk http://www.mkgmap.org.uk/mailman/listinfo/mkgmap-dev
Hi! Christian Gawron schrieb:
A very good hint. The understanding of all discussions on the matter is that the current OSM licence does not allow mixing with CIAT data nor with ASTER data. So you cannot publish those maps. Could you elaborate on this? I would understand that e.g. the license of the ASTER data forbids publishing of derived work (although I haven't checked this). But is there really a point in the OSM license which forbids creating derived works (maps) which as an "add-on" contain features created from proprietary data sources?
Yes. To put it very simple: The share-alike cause of the OSM licence requires you to make everything available under the OSM licence again. If you combine OSM with other data which has _any_ form of restrictions or any different interpretation of freeness, you are breaking one of the licences. A common workaround is to add the combined data in an extra layer and claim this makes a collective work. This is not legally sound, but has never been challenged. I am putting my hope in the new licence which would allow you to create such a combined work. bye Nop
Hi! Christian Gawron schrieb:
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).
Sounds interesting. Does it also distinguish minor and major elevation lines? What determines the area in which the contour lines are calculated? Is there a way to avoid an overload with contours e.g. when processing the whole of Germany inclduing very flat country as well as the the alps? bye Nop
Nop schrieb:
Sounds interesting. Does it also distinguish minor and major elevation lines? Currently the types are still hard-coded as follows: multiples of 200m: 0x22 multiples of 50m: 0x21 0x20 for all others This should of course be configurable (as well as a "feet" mode). What determines the area in which the contour lines are calculated? The bounding box of the tile. Is there a way to avoid an overload with contours e.g. when processing the whole of Germany inclduing very flat country as well as the the alps? You may split the map and use diffeent settings for each tile.
Best wishes Christian
Hi! I am very much interested in your work as I used srtm2osm, which is no longer working, and now I am looking for alternate techniques. Christian Gawron schrieb:
Nop schrieb:
Sounds interesting. Does it also distinguish minor and major elevation lines? Currently the types are still hard-coded as follows: multiples of 200m: 0x22 multiples of 50m: 0x21 0x20 for all others This should of course be configurable (as well as a "feet" mode).
Should not be too difficult to make this configurable and also turn off the elevation texts for each level. No need to label each 10m-line.
What determines the area in which the contour lines are calculated? The bounding box of the tile.
So each map tile must have a bounding box? What happens if there is none? Currently I am not generating them, but that's no problem.
Is there a way to avoid an overload with contours e.g. when processing the whole of Germany inclduing very flat country as well as the the alps? You may split the map and use diffeent settings for each tile.
This is still a problem. When creating a map of Germany, I have 440 tiles, but I cannot tell in advance which of them are too mountaineous too show with all altitude lines. Even if nothing crashes, if you show 10m lines for the whole of the alps, the map size will explode and the map consists mostly of altitude lines. In the past, I have handled it by first creating the contour line .osm with srtm2osm and 10m lines. When the size of the output .osm was too large, I doubled the base distance to 20m lines and tried again, until the output size was within limits. Can you think of a way to integrate an upper limit for elevation lines and a try-and-error fitting against this limit to your calculations? Or even better, is there a way to predict the outcome without doing all calculations first? bye Nop
Nop schrieb:
Hi!
I am very much interested in your work as I used srtm2osm, which is no longer working, and now I am looking for alternate techniques.
Christian Gawron schrieb:
Nop schrieb:
Sounds interesting. Does it also distinguish minor and major elevation lines? Currently the types are still hard-coded as follows: multiples of 200m: 0x22 multiples of 50m: 0x21 0x20 for all others This should of course be configurable (as well as a "feet" mode).
Should not be too difficult to make this configurable and also turn off the elevation texts for each level. No need to label each 10m-line.
Maybe I should use the styling engine here - I will have a look at it.
What determines the area in which the contour lines are calculated? The bounding box of the tile.
So each map tile must have a bounding box? What happens if there is none? Currently I am not generating them, but that's no problem. AFAIK the mapper generates a bounding box, but this may be too large.
Is there a way to avoid an overload with contours e.g. when processing the whole of Germany inclduing very flat country as well as the the alps? You may split the map and use diffeent settings for each tile.
This is still a problem. When creating a map of Germany, I have 440 tiles, but I cannot tell in advance which of them are too mountaineous too show with all altitude lines. Even if nothing crashes, if you show 10m lines for the whole of the alps, the map size will explode and the map consists mostly of altitude lines.
In the past, I have handled it by first creating the contour line .osm with srtm2osm and 10m lines. When the size of the output .osm was too large, I doubled the base distance to 20m lines and tried again, until the output size was within limits. Can you think of a way to integrate an upper limit for elevation lines and a try-and-error fitting against this limit to your calculations? Or even better, is there a way to predict the outcome without doing all calculations first?
I could add a feature which looks at the difference between minimum and maximum height in a tile and chose the interval based on this. Best wishes Christian
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. That seems very cool, and I look forward to trying it. It would be nice to have a few intermediate abilities (I've read most the rest of the thread by now): output osm format data from the DEM (for other use, or later use, or for intermediate reprocessing) 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. For the first point, my immediate reaction was that it was too bad this was in mkgmap instead of a standalone utility, but I see the intermediate storage issue. Having the code in mkgmap is not a big deal, and it enables the nice behavior of an internal pipeline. This isn't really related to your work, but it seems to me that mkgmap has grown tons of options, and modern recommended usage is to have a lot of them on. It would be nice if every option had a negative form and there was a --best-practice option that enabled all the ones that are recommended for making the best garmin maps (complete with tdb etc. for turning into gmapi, and everything else one would need). Right now it seems a little hard to use, where everyone has to chase down those options and write a script. I can certainly see where at one time each of these options might have caused trouble and thus was not default, and the notion that behavior shouldn't change for compatibility. But I also wonder how many users really want that rather than just to make a img file for their receiver.
Great, will try out you patch soon. Is it compatible with srtm1" data? - I would like to create maps from 1" SRTM from viewfinderpanoramas.org How big is the output size relative to other tools, have you compared it to other tools like groundtruth (to mp and then mkgmap, or using cgpsmapper)? SRTM2OSM mainly produced huge output files, much bigger than needed. groundtruth is much better in that respect. What happens when running non void-filled data? 2009/7/3 Christian Gawron <christian.gawron@gmx.de>
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-gri...) 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-grid-contouring-algorithm.pdf">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<https://wist.echo.nasa.gov/%7Ewist/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"); } } _______________________________________________ mkgmap-dev mailing list mkgmap-dev@lists.mkgmap.org.uk http://www.mkgmap.org.uk/mailman/listinfo/mkgmap-dev
Hi! Felix Hartmann schrieb:
SRTM2OSM mainly produced huge output files, much bigger than needed. groundtruth is much better in that respect.
Unfortunately, groundtruth is not yet capable of simply creating contour lines within a given bounding box, it always forcibly breaks it up in segments and changes the bounding box size. It looks like Christians work is a good step towards an srtm2osm alternative. bye Nop
Felix Hartmann wrote:
SRTM2OSM mainly produced huge output files, much bigger than needed. groundtruth is much better in that respect.
Groundtruth needs tons of RAM. I didn't manage to run it for whole Germany data, but SRTM2OSM did that just fine on the same 4 GB RAM machine.
participants (8)
-
Christian Gawron -
Felix Hartmann -
Greg Troxel -
Nop -
Ralf Kleineisel -
Rudi -
Thilo Hannemann -
Toby Speight