This version of the patch (to be applied against rev. 1080) has a HGTDEM which reads SRTM .hgt files and an optional GeoTiffDEM which contains inner classes to read CGIAR and ASTER. To compile GeoTiffDEM.java you'll need Java Advanced Imaging (JAI) installed, so it is disabled in the build.xml. BTW: Does anybody know how to switch the exclude on and off depending on a property? Best wishes Christian Index: src/uk/me/parabola/mkgmap/reader/dem/DEM.java =================================================================== --- src/uk/me/parabola/mkgmap/reader/dem/DEM.java (revision 1080) +++ src/uk/me/parabola/mkgmap/reader/dem/DEM.java (working copy) @@ -17,20 +17,11 @@ 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; @@ -54,7 +45,6 @@ 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 @@ -63,7 +53,7 @@ public abstract class DEM { final static double epsilon = 1e-9; - final static double delta = 1.5; + final protected static double delta = 1.5; final static int maxPoints=200000; final static double minDist = 15; final static double maxDist = 21; @@ -71,18 +61,17 @@ 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; + protected static int M=1200; + protected static int N=M; + protected static double res = 1.0/N; static int id = -1; short values[] = null; - int lat; - int lon; + protected int lat; + protected 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; + protected abstract double ele(int x, int y); + protected abstract void read(int minLon, int minLat, int maxLon, int maxLat); public static void createContours(LoadableMapDataSource mapData, EnhancedProperties config) { @@ -95,22 +84,33 @@ System.out.printf("bounds: %f %f %f %f\n", minLat, minLon, maxLat, maxLon); DEM data; - String demType = config.getProperty("dem-type", "CGIAR"); + String demType = config.getProperty("dem-type", "SRTM"); String dataPath; - if (demType.equals("ASTER")) { - dataPath = config.getProperty("dem-path", "ASTER"); - data = new ASTERGeoTiff(dataPath, minLat, minLon, maxLat, maxLon); + Class demClass; + try { + if (demType.equals("ASTER")) { + dataPath = config.getProperty("dem-path", "ASTER"); + demClass = Class.forName("uk.me.parabola.mkgmap.reader.dem.optional.GeoTiffDEM$ASTER"); + } + else if (demType.equals("CGIAR")) { + dataPath = config.getProperty("dem-path", "CGIAR"); + demClass = Class.forName("uk.me.parabola.mkgmap.reader.dem.optional.GeoTiffDEM$CGIAR"); + } + else { + dataPath = config.getProperty("dem-path", "SRTM"); + demClass = Class.forName("uk.me.parabola.mkgmap.reader.dem.HGTDEM"); + } + java.lang.reflect.Constructor<DEM> constructor = demClass.getConstructor(String.class, + Double.TYPE, Double.TYPE, + Double.TYPE, Double.TYPE); + data = constructor.newInstance(dataPath, minLat, minLon, maxLat, maxLon); } - else if (demType.equals("CGIAR")) { - dataPath = config.getProperty("dem-path", "CGIAR"); - data = new CGIARGeoTiff(dataPath, minLat, minLon, maxLat, maxLon); + catch (Exception ex) + { + throw new RuntimeException("failed to create DEM", ex); } - 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); @@ -181,213 +181,11 @@ } } - 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; @@ -1184,78 +982,6 @@ 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 { Index: src/uk/me/parabola/mkgmap/reader/MapperBasedMapDataSource.java =================================================================== --- src/uk/me/parabola/mkgmap/reader/MapperBasedMapDataSource.java (revision 1080) +++ src/uk/me/parabola/mkgmap/reader/MapperBasedMapDataSource.java (working copy) @@ -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 1080) +++ src/uk/me/parabola/mkgmap/reader/osm/Way.java (working copy) @@ -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. Index: build.xml =================================================================== --- build.xml (revision 1080) +++ build.xml (working copy) @@ -72,7 +72,7 @@ <include name="**/*.java" /> <classpath refid="main"/> <!-- <compilerarg value="-Xlint:unchecked"/> --> - <exclude name="**/reader/dem/*.java"/> + <exclude name="**/reader/dem/optional/*.java"/> </javac> </target> /* * 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 static java.nio.channels.FileChannel.MapMode.READ_ONLY; public 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"); } } /* * 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.optional; import uk.me.parabola.mkgmap.reader.dem.DEM; import java.awt.image.renderable.ParameterBlock; import java.awt.image.Raster; import java.awt.image.DataBuffer; import java.io.*; import java.nio.channels.FileChannel; import java.nio.MappedByteBuffer; import com.sun.media.jai.codec.*; import javax.media.jai.*; public abstract class GeoTiffDEM extends DEM { Raster raster; String fileName; int minLat, minLon, maxLat, maxLon; PlanarImage image; void initData() { 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 static class CGIAR extends GeoTiffDEM { public CGIAR(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); System.out.printf("CGIAR GeoTIFF: %s\n", fileName); N = 6000; M = 6000; res = 5.0/M; initData(); } 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"); } protected 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); } 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 ASTER extends GeoTiffDEM { public ASTER(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); System.out.printf("ASTER GeoTIFF: %s\n", fileName); N = 3600; M = 3600; res = 1.0/M; initData(); } 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"); } protected 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); } 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; } } } }