Am 22.02.2010 23:07, schrieb Ronny Klier:
This patch adds support for maps crossing borders of SRTM (hgt) files.
So its possible to create contour lines for larger areas.
Additional:
- description of contour relevant options
- "dem-increment" option has been modified so it is possible to create
different contour levels (major, medium, minor)
Open issues:
- I am still struggling with using multiple GeoTiff files (CGIAR + ASTER).
- data files have to be loaded manually.
Ronny
Finally got it to work with CGIAR data.
The patch is yet untested with ASTER data. If anyone has such data
please try this patch.
Index: resources/help/en/options
===================================================================
--- resources/help/en/options (revision 1580)
+++ resources/help/en/options (working copy)
@@ -370,11 +370,49 @@
--tdbfile
Write a .tdb file.
---show-profiles=1
+--show-profiles=0|1
Sets a flag in tdb file which marks set mapset as having contour
lines and allows showing profile in MapSource. Default is 0
which means disabled.
-
+ DON'T use this option when combining separate tiles with contour
+ lines and real data tiles for the same area. MapSource will crash
+ on them.
+
+--contour
+ Create contour lines. Behaviour is specified by the dem-* options.
+
+--dem-type=[SRTM|ASTER|CGIAR]
+ Defines kind of processed data. Default is SRTM.
+ SRTM data covers earth surface from 60° north to 56° south in 1°x1°
+ plain data files (1201x1201 short int values).
+ See http://www2.jpl.nasa.gov/srtm/
+ CGIAR contains optimized SRTM to fill missing data areas.
+ Data is provided as 5°x5° GeoTiff files.
+ See http://srtm.csi.cgiar.org/
+ ASTER data covers earth surface from 83° north to 83° south and is
+ provided as 1°x1° GeoTiff files.
+ See http://www.gdem.aster.ersdac.or.jp/
+
+--dem-path
+ Location of DEM source data. Currently mkgmap does not load required
+ files from the web.
+
+--dem-increment=10[,0[,0]]
+ Create contour lines for steps with specified increment in meter.
+ The really used increment depends on option --dem-maxlevels.
+ The optional parameters define a factor which contour lines are
+ marked as major (first parameter) and as medium (second). If both
+ are omitted all lines are just marked as "contour=elevation".
+ Example: 20,10,5 results in minor contour every 20 meter, medium
+ contour every 100 meter and major contour every 200 meter.
+
+--dem-maxlevels=100
+ Defines maximum number of different elevation levels. Increment
+ will be doubled while maxlevel is exceeded. Using the default
+ values of 10 for dem-increment and 100 for dem-maxlevels, the
+ active increment will be set to 20 if lowest point of tile is 0
+ and highest point is 1200: 10*100 < 1200 < 20*100
+
--draw-priority=25
When two maps cover the same area, this option controls what
order they are drawn in and therefore which map is on top of
Index: src/uk/me/parabola/mkgmap/reader/dem/DEM.java
===================================================================
--- src/uk/me/parabola/mkgmap/reader/dem/DEM.java (revision 1580)
+++ src/uk/me/parabola/mkgmap/reader/dem/DEM.java (working copy)
@@ -28,6 +28,7 @@
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.log.Logger;
import uk.me.parabola.mkgmap.build.MapBuilder;
import uk.me.parabola.mkgmap.general.LevelInfo;
import uk.me.parabola.mkgmap.general.LoadableMapDataSource;
@@ -38,6 +39,7 @@
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.reader.osm.xml.Osm5XmlHandler;
import uk.me.parabola.util.EnhancedProperties;
@@ -47,6 +49,8 @@
* Adaptive Grid Contouring Algorithm</a> by Downing and Zoraster.
*/
public abstract class DEM {
+ private static final Logger log = Logger.getLogger(DEM.class);
+
final static double epsilon = 1e-9;
final protected static double delta = 1.5;
final static int maxPoints = 200000;
@@ -104,11 +108,43 @@
}
Isolines lines = data.new Isolines(data, minLat, minLon,
maxLat, maxLon);
- int increment = config.getProperty("dem-increment", 10);
+ // get base increment between isolines (meter)
+ int increment = 10;
+ // which levels are marked as major elevation lines (default is
all)
+ int maj_increment = 0;
+ // which levels are marked as medium elevation lines (default
is all)
+ int med_increment = 0;
+ String strInc = config.getProperty("dem-increment", null);
+ if (strInc != null) {
+ String steps[] = strInc.split(",");
+ try {
+ increment = Integer.parseInt(steps[0]);
+ } catch (NumberFormatException ex) {
+ log.error("dem-increment expects list of
numerical data");
+ }
+ if (steps.length > 1) {
+ try {
+ maj_increment =
Integer.parseInt(steps[1]);
+ } catch (NumberFormatException ex) {
+ log.error("dem-increment expects list
of numerical data");
+ }
+ }
+ if (steps.length > 2) {
+ try {
+ med_increment =
Integer.parseInt(steps[1]);
+ } catch (NumberFormatException ex) {
+ log.error("dem-increment expects list
of numerical data");
+ }
+ }
+ }
double minHeight = lines.getMinHeight();
double maxHeight = lines.getMaxHeight();
+ // get maximum number of isolines
int maxLevels = config.getProperty("dem-maxlevels", 100);
+ // adjust increment to ensure no more than maxLevels isolines
are created
+ // e.g. if we have a difference between lowest and highest
point of 1980 increment is increased to 20.
+ // So the only way to ensure fixed increments is to define
dem-maxlevels to a value which will be reached never.
while ((maxHeight - minHeight) / increment > maxLevels)
increment *= 2;
@@ -149,6 +185,14 @@
for (Isolines.Isoline line : lines.isolines) {
Way way = new Way(id--, line.points);
way.addTag("contour", "elevation");
+ if (maj_increment != 0) {
+ if ((level / increment) % maj_increment
== 0)
+ way.addTag("contour_ext",
"elevation_major");
+ else if (med_increment != 0 && ((level
/ increment) % med_increment == 0))
+ way.addTag("contour_ext",
"elevation_medium");
+ else
+ way.addTag("contour_ext",
"elevation_minor");
+ }
way.addTag("ele", String.format("%d", (int)
line.level));
converter.convertWay(way);
}
@@ -682,15 +726,26 @@
}
private void init() {
+ // check for valid bounds/prevent negativ bounds
+ if (minX < 2)
+ minX = 2;
+ if (minY < 2)
+ minY = 2;
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
+ data.read(minX - 2, minY - 2, maxX + 2, maxY + 2);
max = -1000;
min = 10000;
+ double e;
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);
+ e = data.elevation(i, j);
+ if (e < - 1000) // SRTM marks invalid
values with Short.MIN_VALUE, so don't get confused by this
+ continue;
+ if (e < min)
+ min = e;
+ if (e > max)
+ max = e;
}
debug("min: %f, max: %f\n", min, max);
Index: src/uk/me/parabola/mkgmap/reader/dem/HGTDEM.java
===================================================================
--- src/uk/me/parabola/mkgmap/reader/dem/HGTDEM.java (revision 1580)
+++ src/uk/me/parabola/mkgmap/reader/dem/HGTDEM.java (working copy)
@@ -22,44 +22,69 @@
import static java.nio.channels.FileChannel.MapMode.READ_ONLY;
+import java.util.HashMap;
+import java.util.Map;
+
public class HGTDEM extends DEM
{
- MappedByteBuffer buffer = null;
-
+ Map<Integer, MappedByteBuffer> buffers = new HashMap<Integer,
MappedByteBuffer>();
+
+ int nLon; // number of touched files in longitude direction
+ int nLat; // number of touched files in latitude direction
+
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);
- }
+ res = 1.0 / M;
+ this.lat = (int) minLat;
+ this.lon = (int) minLon;
+ nLat = (int)maxLat - this.lat + 1;
+ nLon = (int)maxLon - this.lon + 1;
+ if (nLat > nLon)
+ N = M * nLat;
+ else
+ N = M * nLon;
+ for (int x = 0; x < nLon; ++x) {
+ String eastWest = "E";
+ int actLon = lon + x;
+ if ((lon + x) < 0) {
+ actLon = -(lon + x);
+ eastWest = "W";
+ }
+ for (int y = 0; y < nLat; ++y) {
+ String northSouth = (lat + y) < 0 ? "S" : "N";
+ String fileName = String.format("%s/%s%02d%s%03d.hgt",
dataPath,
+ northSouth, (lat + y) <
0 ? -(lat + y) : (lat + y),
+ eastWest, actLon);
+ try {
+ FileInputStream is = new
FileInputStream(fileName);
+ buffers.put(x * nLat + y,
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)
+ @Override
+ protected 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;
- }
+ @Override
+ protected double ele(int x, int y) {
+ int actLon = x / M;
+ int actLat = y / M;
+ if (actLon > nLon || actLat > nLat)
+ throw new IndexOutOfBoundsException("point not inside
area");
+ MappedByteBuffer buffer = buffers.get(actLon*nLat+actLat);
+ return buffer.getShort(2*((M-y%M)*(M+1)+x%M))+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");
+ out.write(" <copyright>\n");
+ out.write(" Contour lines generated from DEM data by NASA\n");
+ out.write(" </copyright>\n");
}
}
Index: src/uk/me/parabola/mkgmap/reader/dem/optional/GeoTiffDEM.java
===================================================================
--- src/uk/me/parabola/mkgmap/reader/dem/optional/GeoTiffDEM.java
(revision 1580)
+++ src/uk/me/parabola/mkgmap/reader/dem/optional/GeoTiffDEM.java
(working copy)
@@ -20,60 +20,157 @@
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 java.util.HashMap;
+import java.util.Map;
import com.sun.media.jai.codec.*;
import javax.media.jai.*;
public abstract class GeoTiffDEM extends DEM
{
- Raster raster;
- String fileName;
+ /*
+ * GeoTiff represents a single GeoTiff file
+ */
+ class GeoTiff
+ {
+ Raster raster;
+ String fileName;
+ PlanarImage image;
+
+ GeoTiff(String name) {
+ fileName = name;
+ }
+
+ 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);
+ }
+ }
+
+ /*
+ * read required part of image file.
+ *
+ * (minLon, minLat) start point of data relative to west-north
corner of tile
+ */
+ void read(int minLon, int minLat, int sizeLon, int sizeLat) {
+ raster = image.getData(new java.awt.Rectangle(minLon, minLat,
sizeLon, sizeLat));
+ System.out.printf("read: %d %d %d %d\n", minLon,
minLat, sizeLon, sizeLat);
+ }
+
+ double ele(int x, int y) {
+ try {
+ int elevation = raster.getPixel(x, M-y-1,
(int[])null)[0];
+ return elevation+delta;
+ }
+ catch (ArrayIndexOutOfBoundsException ex) {
+ System.out.printf("ele: (%d, %d) (%d, %d, %d, %d)
%s\n",
+ x, M-y-1,
+ raster.getMinX(), raster.getMinY(),
+ raster.getWidth(),
raster.getHeight(), ex.toString());
+ throw ex;
+ }
+ }
+ }
+
+ Map<Integer, GeoTiff> files = new HashMap<Integer, GeoTiff>();
int minLat, minLon, maxLat, maxLon;
- PlanarImage image;
+ int nLon; // number of touched files in longitude direction
+ int nLat; // number of touched files in latitude direction
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());
+ for (GeoTiff tiff : files.values()) {
+ tiff.initData();
+ }
+ }
+
+ protected void read(int minLon, int minLat, int maxLon, int maxLat)
+ {
+ this.minLon = minLon;
+ this.minLat = minLat;
+ this.maxLon = maxLon;
+ this.maxLat = maxLat;
+ // dem data starts in west-north corner, so latitude has to be
reversed
+ for (int x = 0; x < nLon; ++x) {
+ // start at 0 except first map where to start at minLon
+ int actLon = (x == 0) ? minLon : 0;
+ // read all values except for last file. For this read only to
maxLon
+ int sizeLon = M - actLon + 1;
+ if (x == (nLon - 1))
+ sizeLon = maxLon - x * M - actLon + 1;
+ for (int y = 0; y < nLat; ++y) {
+ int actLat = 0; // start at northern edge of dem tile
+ if (y == (nLat - 1)) {
+ // for last dem tile start at northern edge of
requested data
+ actLat = nLat * M - maxLat - 1;
+ }
+ int sizeLat = M - actLat + 1; // read to southern edge
of tile
+ if (y == 0) {
+ // for first tile read only to southern edge of
requested data
+ sizeLat -= minLat;
+ }
+ files.get(x*nLat+y).read(actLon, actLat, sizeLon,
sizeLat);
+ }
+ }
}
- catch (Exception e) {
- throw new RuntimeException(e);
+
+ protected double ele(int x, int y)
+ {
+ int iLon;
+ int iLat = nLat;
+ for (iLon = 0; iLon < nLon; ++iLon, x-=M) {
+ if (x >= M)
+ continue;
+ for (iLat = 0; iLat < nLat; ++iLat, y-=M) {
+ if (y >= M)
+ continue;
+ break;
+ }
+ break;
+ }
+ if (iLon >= nLon || iLat >= nLat)
+ throw new IndexOutOfBoundsException("point not inside
area");
+ return files.get(iLon*nLat+iLat).ele(x, y);
}
- }
- public static class CGIAR extends GeoTiffDEM
+ 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;
+ this.nLat = (((int) (maxLat/5))*5 - this.lat) / 5 + 1;
+ this.nLon = (((int) (maxLon/5))*5 - this.lon) / 5 + 1;
M = 6000;
res = 5.0/M;
-
+ if (nLat > nLon)
+ N = M * nLat;
+ else
+ N = M * nLon;
+ for (int x = 0; x < nLon; ++x) {
+ int tileX = (180 + lon)/5 + x + 1;
+ for (int y = 0; y < nLat; ++y) {
+ int tileY = (60 - lat)/5 - y; // we have to count
reverse
+ String fileName =
String.format("%s/srtm_%02d_%02d.tif", dataPath, tileX, tileY);
+ files.put(x * nLat + y, new GeoTiff(fileName));
+ System.out.printf("CGIAR GeoTIFF: %s\n", fileName);
+ }
+ }
initData();
}
@@ -84,30 +181,6 @@
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
@@ -117,19 +190,30 @@
{
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;
+ nLat = (int)maxLat - this.lat + 1;
+ nLon = (int)maxLon - this.lon + 1;
+ M = 3601;
res = 1.0/M;
+ if (nLat > nLon)
+ N = M * nLat;
+ else
+ N = M * nLon;
+ for (int x = 0; x < nLon; ++x) {
+ String eastWest = "E";
+ int actLon = lon + x;
+ if ((lon + x) < 0) {
+ actLon = -(lon + x);
+ eastWest = "W";
+ }
+ for (int y = 0; y < nLat; ++y) {
+ String northSouth = (lat + y) < 0 ? "S" : "N";
+ String fileName =
String.format("%s/ASTGTM_%s%02d%s%03d_dem.tif", dataPath,
+ northSouth, (lat + y) < 0
? -(lat + y) : (lat + y),
+ eastWest, actLon);
+ files.put(x * nLat + y, new GeoTiff(fileName));
+ System.out.printf("ASTER GeoTIFF: %s\n", fileName);
+ }
+ }
initData();
}
@@ -140,30 +224,6 @@
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;
- }
- }
}
}
_______________________________________________
mkgmap-dev mailing list
[email protected]
http://www.mkgmap.org.uk/mailman/listinfo/mkgmap-dev
