July 2009 - mkgmap-dev - The mkgmap lists

R: [mkgmap-dev] [PATCH v2] - beware of the bollards!
by marco_certelli＠yahoo.it 07 Jul '09

07 Jul '09

Mark, may I suggest a totally different approach? why not to use a couple of restriction relations? a no_straight_on between the segment before (from) and after (to) the bollard node (via) and a second no_straight_on relation the way round. Of course with the right "except": bicycle, foot, ... You could also invent an "internal" relation restriction like "no_passage" The issue here seems that the "except" key does not foresees "foot". Would it be feasible? Ciao. --- Mar 7/7/09, Mark Burton <markb(a)ordern.com> ha scritto: > Da: Mark Burton <markb(a)ordern.com> > Oggetto: [mkgmap-dev] [PATCH v2] - beware of the bollards! > A: mkgmap-dev(a)lists.mkgmap.org.uk > Data: Martedì 7 luglio 2009, 16:58 > > v2 - quick update based on instantaneous feedback from ML! > > Now works for any POI that sets "access=no" (could use the > more general > test of any of the access tags being set but let's see how > this works > for now). > > Default access rights now set in style file. > > Any suggestions for a better code for cycle_barrier? > > ------------ > > Fed up of being routed in your car down city streets only > to find the > way is blocked by a bollard? Well, if so, this is the patch > for you. > > If a way has a bollard on a point, the segments of the way > that > connect to the bollard have access restrictions placed on > them. By > default, a bollard implies: access=no, foot=yes, > bicycle=yes. > > Testing using mapsource shows that it generally works as > expected > although if the destination cannot be reached by any other > route, it > routes straight through the bollard rather than failing! If > the > destination can be reached by some other route, even if the > route is > really long, it will avoid the bollard. > > I have chosen Garmin code 0x660f (pillar) for the bollard. > On my etrex > it appears as a small dot in the way. > > As usual, all feedback is welcome. > > Cheers, > > Mark > > -----Segue allegato----- > > _______________________________________________ > mkgmap-dev mailing list > mkgmap-dev(a)lists.mkgmap.org.uk > http://www.mkgmap.org.uk/mailman/listinfo/mkgmap-dev

2 1

[PATCH v2] - beware of the bollards!
by Mark Burton 07 Jul '09

07 Jul '09

v2 - quick update based on instantaneous feedback from ML! Now works for any POI that sets "access=no" (could use the more general test of any of the access tags being set but let's see how this works for now). Default access rights now set in style file. Any suggestions for a better code for cycle_barrier? ------------ Fed up of being routed in your car down city streets only to find the way is blocked by a bollard? Well, if so, this is the patch for you. If a way has a bollard on a point, the segments of the way that connect to the bollard have access restrictions placed on them. By default, a bollard implies: access=no, foot=yes, bicycle=yes. Testing using mapsource shows that it generally works as expected although if the destination cannot be reached by any other route, it routes straight through the bollard rather than failing! If the destination can be reached by some other route, even if the route is really long, it will avoid the bollard. I have chosen Garmin code 0x660f (pillar) for the bollard. On my etrex it appears as a small dot in the way. As usual, all feedback is welcome. Cheers, Mark

1 0

R: [mkgmap-dev] [PATCH v1] - beware of the bollards!
by Marco Certelli 07 Jul '09

07 Jul '09

--- Mar 7/7/09, Mark Burton <markb(a)ordern.com> ha scritto: > Da: Mark Burton <markb(a)ordern.com> > Oggetto: [mkgmap-dev] [PATCH v1] - beware of the bollards! > A: mkgmap-dev(a)lists.mkgmap.org.uk > Data: Martedì 7 luglio 2009, 15:30 > > Fed up of being routed in your car down city streets only > to find the > way is blocked by a bollard? Well, if so, this is the patch > for you. > > If a way has a bollard on a point, the segments of the way > that > connect to the bollard have access restrictions placed on > them. By > default, a bollard implies: access=no, foot=yes, > bicycle=yes. > > Testing using mapsource shows that it generally works as > expected > although if the destination cannot be reached by any other > route, it > routes straight through the bollard rather than failing! If > the > destination can be reached by some other route, even if the > route is > really long, it will avoid the bollard. > > I have chosen Garmin code 0x660f (pillar) for the bollard. > On my etrex > it appears as a small dot in the way. > > As usual, all feedback is welcome. Mark, I read on wiki that cycle_barrier has no "default" access rule. I think you should be more conservative and, in absence of explicit tags, let bicycles passing through (wiki says that cycle_barrier should only imply motor_vehicle=no) Ciao. > > Cheers, > > Mark > > -----Segue allegato----- > > _______________________________________________ > mkgmap-dev mailing list > mkgmap-dev(a)lists.mkgmap.org.uk > http://www.mkgmap.org.uk/mailman/listinfo/mkgmap-dev

2 1

resource usage?
by Greg Troxel 06 Jul '09

06 Jul '09

I am trying to make an img/tdb/overview from massachusetts.osm (downloaded from cloudmade). My computer is a mac with 2G of ram. The .osm file is 2.4G. While running with 1024m heap size java ran out of heap. With 2048, it is moving along, but paging a lot wiht a RSIZE of 1750M and getting only about 2% on cpu, and now up to 4m35s of runtime. Should I expect this to finish? I would like a routable map - the cloudmade ones aren't routable, and I'd like to be able to tweak styles to emphasize hiking trails. I think it would be nice to give some rules of thumb about input sizes somewhere, prehaps in a README in the sources, in terms of what the max sizes that work are, how long they take, and what the needed RAM/RSS is. It would also be cool to have an output file or stdout that counts the resources used somehow.

1 0

Re: mkgmap Contours
by Christian Gawron 06 Jul '09

06 Jul '09

The patch was disabled because it introduces a dependency on JAI. I'm just creating a new version which just includes a reader for HGT files and makes the CGIAR/ASTER code (which uses JAI to read the GeoTIFFs) optional. Best wishes Christian Stuart Poulton schrieb: > Hi Christian.... > > I'm trying the mkgmap-r1080 binary release, trying to use the contour > tools > > Here's the command line: > > java -Xmx2000m -jar mkgmap.jar --contours --dem-type=SRTM > --dem-path=/home/q-maps/data/SrtmCache/ --dem-separate-img *.osm > > SRTM directory: > > gateway:/home/g-dev/mkgmap-r1080 # ls -l /home/q-maps/data/SrtmCache/ > total 11296 > -rw-r--r-- 1 root root 2884802 May 30 08:19 N52W003.hgt > -rw-r--r-- 1 root root 2884802 May 30 08:18 N52W004.hgt > -rw-r--r-- 1 root root 2884802 May 30 08:20 N53W003.hgt > -rw-r--r-- 1 root root 2884802 May 30 08:19 N53W004.hgt > > I don't seem to get any contours, clearly I'm missing something :( Any > pointers ? > > Cheers > > Stuart

3 3

[PATCH v1] - clone ways that contain some tag pairs
by Mark Burton 06 Jul '09

06 Jul '09

Further to my earlier email about the problem of ways that have both highway and boundary tags defined, here's a patch that possibly provides a workaround. It lets you specify pairs of tags. If a way has both of the tags defined, it is cloned and one of the clones has one of the tags and the other clone has the other tag. So for the example in question, you could specify the option: --clone-ways-with-these-tags=highway/boundary Another possible tag pair could be waterway/boundary All feedback is welcome. Cheers, Mark

3 4

Splitter on MacOS?
by Dermot McNally 06 Jul '09

06 Jul '09

Hi folks, Up until now, I've never bothered using the splitter (you can get all of Ireland onto one .img file). But now I want to build a map of NL in anticipation of SOTM. I'm on MacOS Leopard and can't seem to make the splitter work. When I try to use "ant dist" to build from latest svn I get the following: compile: [javac] Compiling 16 source files to /Users/dermot/Mapping/osm/Garmin Maps/routable/splitter/trunk/build/classes [javac] /Users/dermot/Mapping/osm/Garmin Maps/routable/splitter/trunk/src/uk/me/parabola/splitter/AreaList.java:73: cannot find symbol [javac] symbol : variable ROOT [javac] location: class java.util.Locale [javac] pw.format(Locale.ROOT, "%d: %d,%d to %d,%d\n", [javac] ^ [javac] /Users/dermot/Mapping/osm/Garmin Maps/routable/splitter/trunk/src/uk/me/parabola/splitter/AreaList.java:77: cannot find symbol [javac] symbol : variable ROOT [javac] location: class java.util.Locale [javac] pw.format(Locale.ROOT, "# : %f,%f to %f,%f\n", [javac] ^ [javac] /Users/dermot/Mapping/osm/Garmin Maps/routable/splitter/trunk/src/uk/me/parabola/splitter/AreaList.java:110: cannot find symbol [javac] symbol : method isEmpty() [javac] location: class java.lang.String [javac] if (line.isEmpty() || line.charAt(0) == '#') [javac] ^ [javac] /Users/dermot/Mapping/osm/Garmin Maps/routable/splitter/trunk/src/uk/me/parabola/splitter/SubArea.java:93: cannot find symbol [javac] symbol : variable ROOT [javac] location: class java.util.Locale [javac] String filename = new Formatter().format(Locale.ROOT, "%08d.osm.gz", mapid).toString(); [javac] ^ [javac] 4 errors Trying the prebuilt jar file doesn't work either: Exception in thread "main" java.lang.UnsupportedClassVersionError: Bad version number in .class file at java.lang.ClassLoader.defineClass1(Native Method) at java.lang.ClassLoader.defineClass(ClassLoader.java:675) at java.security.SecureClassLoader.defineClass(SecureClassLoader.java:124) at java.net.URLClassLoader.defineClass(URLClassLoader.java:260) at java.net.URLClassLoader.access$100(URLClassLoader.java:56) at java.net.URLClassLoader$1.run(URLClassLoader.java:195) at java.security.AccessController.doPrivileged(Native Method) at java.net.URLClassLoader.findClass(URLClassLoader.java:188) at java.lang.ClassLoader.loadClass(ClassLoader.java:316) at sun.misc.Launcher$AppClassLoader.loadClass(Launcher.java:280) at java.lang.ClassLoader.loadClass(ClassLoader.java:251) at java.lang.ClassLoader.loadClassInternal(ClassLoader.java:374) I daresay the reasons for these failures will be obvious to somebody with more idea of java than me, so apologies if I'm missing something obvious. Cheers, Dermot -- -------------------------------------- Iren sind menschlich

3 4

Commit: r1080: Exclude the dem files as they do not compile without extra packages.
by svn commit 06 Jul '09

06 Jul '09

Version 1080 was commited by steve on 2009-07-06 00:03:12 +0100 (Mon, 06 Jul 2009) Exclude the dem files as they do not compile without extra packages.

2 1

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

05 Jul '09

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

8 23

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

05 Jul '09

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

1 0