断片的な海岸線の線分を繋ぎ合わせて、閉ループを完成させるまでをサポートしている。
/*
javac SAX.java
java -Xmx5000m SAX c:/osm/kanto.osm
java -Xmx5000m SAX c:/osm/japan.osm
*/
import java.io.IOException;
import java.io.InputStream;
import java.io.BufferedInputStream;
import java.nio.file.Files;
import java.nio.file.Paths;
import javax.xml.parsers.ParserConfigurationException;
import javax.xml.parsers.SAXParser;
import javax.xml.parsers.SAXParserFactory;
import org.xml.sax.SAXException;
import org.xml.sax.Attributes;
import org.xml.sax.helpers.DefaultHandler;
import java.util.*;
import java.util.stream.Collectors;
public class SAX {
static int sumNodes = 0;
static int cntWays = 0;
static List<long[]> listPolygons = new ArrayList<>();
static HashMap<Long, long[]> mapCoastlines = new HashMap<>(5000);
// 海岸線のノードIDから経度・緯度を得る
static HashMap<Long, Long> mapId2LonLat = new HashMap<>(5*1024*1024);
// パス１
static class SAX1 extends DefaultHandler {
boolean fCoastline = false;
String section = "node";
List<Long> listNodes = new ArrayList<>();
long id;
@Override
public void startElement (String namespaceURI,
String localName, String qName, Attributes atts ) {
if (qName.equals("node")) {
return; // パス１では nodeセクションをスキップ
} else if (qName.equals("way")) {
section = "way";
fCoastline = false;
} else if (section.equals("node") && qName.equals("tag")) {
return; // パス１のnodeセクションでは tag も無視する
} else if (qName.equals("osm")) {
return;
}
if (qName.equals("nd")) {
long ndRef = Long.parseLong(atts.getValue("ref"));
listNodes.add(ndRef);
}
if (qName.equals("way")) {
//System.out.println( "Start Element: " + qName );
listNodes.clear(); // 再利用
}
if (section.equals("way") && qName.equals("tag")) {
String key = atts.getValue("k");
String val = atts.getValue("v");
if (key.equals("natural") && val.equals("coastline")) {
fCoastline = true;
}
}
}
@Override
public void endElement(String namespaceURI, String localName, String qName) {
if (!qName.equals("way") || !fCoastline) {
return;
}
sumNodes += listNodes.size();
if (++cntWays % 1000 == 0) System.out.printf("sumNodes=%d\n", sumNodes);
long[] nodes = new long[listNodes.size()];
for (int i = 0; i < nodes.length; i++) {
long ndRef = listNodes.get(i);
mapId2LonLat.put(ndRef, null); // 海岸線ノードであることを登録する
nodes[i] = ndRef;
}
long head = nodes[0];
long tail = nodes[nodes.length-1];
if (head == tail) { // 単独で陸地ポリゴン
listPolygons.add(nodes);
return;
}
if (!mapCoastlines.containsKey(head) && !mapCoastlines.containsKey(tail)) {
mapCoastlines.put(head, nodes); // 先頭ノードをキー
mapCoastlines.put(tail, nodes); // 末尾ノードをキー
return;
}
while ((mapCoastlines.containsKey(head) || mapCoastlines.containsKey(tail))) {
long[] src = mapCoastlines.get(head);
if (src == null) src = mapCoastlines.get(tail);
mapCoastlines.remove(src[0]);
mapCoastlines.remove(src[src.length-1]); // 上と同じ配列
if (src[0] == head || src[src.length-1] == tail) {
reverse(nodes);
head = nodes[0];
tail = nodes[nodes.length-1];
} // 反転する
if (src[src.length-1] == head && src[0] == tail) {
long[] dst = new long[src.length + nodes.length - 2];
System.arraycopy(src, 0, dst, 0, src.length);
System.arraycopy(nodes, 1, dst, src.length, nodes.length-2);
listPolygons.add(dst);
System.out.printf( "閉ループ完成 nodes=%d mapCoastlines数=%d listPolygons数=%d\n",
dst.length, mapCoastlines.size(), listPolygons.size());
return; // 閉ループになる
}
long[] dst = new long[src.length + nodes.length - 1];
if (src[src.length-1] == head) {
System.arraycopy(src, 0, dst, 0, src.length);
System.arraycopy(nodes, 1, dst, src.length, nodes.length-1);
// head は src配列の末尾と同じため、コピーしない
} else if (src[0] == tail) {
System.arraycopy(nodes, 0, dst, 0, nodes.length-1);
System.arraycopy(src, 0, dst, nodes.length-1, src.length);
// tail は src配列の先頭と同じため、コピーしない
}
if (dst[0] == dst[dst.length-1]) { // 閉ループになった
listPolygons.add(dst);
System.out.printf( "閉ループ完成 nodes=%d mapCoastlines数=%d listPolygons数=%d\n",
dst.length, mapCoastlines.size(), listPolygons.size());
return;
}
nodes = Arrays.copyOf(dst, dst.length);
head = nodes[0];
tail = nodes[nodes.length-1];
}
mapCoastlines.put(head, nodes); // HashMapに再登録
mapCoastlines.put(tail, nodes);
System.out.printf( "再登録 dst=%d 登録数=%d\n", nodes.length, mapCoastlines.size() );
}
} // SAX1
//================ パス2 =================//
static class SAX2 extends DefaultHandler {
@Override
public void startElement(String namespaceURI, String localName, String qName, Attributes atts) {
if (atts != null && qName.equals("node")) {
long nid = Long.parseLong(atts.getValue("id"));
if (mapId2LonLat.containsKey(nid)) {
int lon = (int)(Double.parseDouble(atts.getValue("lon"))*10000000);
int lat = (int)(Double.parseDouble(atts.getValue("lat"))*10000000);
long lonlat = (((long)lon)<<32) + lat;
mapId2LonLat.put(nid, lonlat); // キー: ノードID、値: 経度、緯度
} //System.out.println( "lat=" + lat + "lon=" + lon );
}
}
} // SAX2
// 元の配列を並び替える
static void reverse(long[] array) {
for (int left=0, right=array.length-1; left < right; ) {
long temp = array[left];
array[left++] = array[right];
array[right--] = temp;
}
}
static class PointF {
static PointF[] japans = { // 日本地図の領域
new PointF(1.538901E+02f, 2.638211E+01f),
new PointF(1.321529E+02f, 2.646809E+01f),
new PointF(1.316915E+02f, 2.120992E+01f),
new PointF(1.225954E+02f, 2.351966E+01f),
new PointF(1.225607E+02f, 2.584146E+01f),
new PointF(1.288145E+02f, 3.474835E+01f),
new PointF(1.293966E+02f, 3.509403E+01f),
new PointF(1.353079E+02f, 3.754740E+01f),
new PointF(1.405769E+02f, 4.570648E+01f),
new PointF(1.491891E+02f, 4.580245E+01f),
new PointF(1.538901E+02f, 2.638211E+01f)
};
float X;
float Y;
public PointF(float x, float y) {
X = x;
Y = y;
}
static boolean PointInJapan(PointF p) {
return PointInPolygon(p, japans);
}
static boolean PointInPolygon(PointF p, PointF[] poly) {
PointF p1, p2;
boolean inside = false;
PointF oldPoint = poly[poly.length-1];
for (int i = 0; i < poly.length; i++) {
PointF newPoint = poly[i];
if (newPoint.X > oldPoint.X) { p1 = oldPoint; p2 = newPoint; }
else { p1 = newPoint; p2 = oldPoint; }
if ((p1.X < p.X) == (p.X <= p2.X) &&
(p.Y-p1.Y)*(p2.X-p1.X) < (p2.Y-p1.Y)*(p.X-p1.X)) {
inside = !inside;
}
oldPoint = newPoint;
}
return inside;
}
}
public static void main( String[] args ) {
long start = System.currentTimeMillis();
String src = args.length >= 1 ? args[0] : "c:/osm/miyakejima.osm";
System.out.println("==== パス１ ====");
try ( BufferedInputStream bis =
new BufferedInputStream(Files.newInputStream(Paths.get(src))) ) {
SAXParserFactory factory = SAXParserFactory.newInstance();
SAXParser parser = factory.newSAXParser();
SAX1 handler = new SAX1();
parser.parse(bis, handler);
parser = factory.newSAXParser();
} catch ( ParserConfigurationException | SAXException | IOException e ) {
System.out.println( e.getMessage() );
}
System.out.printf("sumNodes=%d\n", sumNodes);
System.out.println("==== パス２ ====");
try ( BufferedInputStream bis =
new BufferedInputStream(Files.newInputStream(Paths.get(src))) ) {
SAXParserFactory factory = SAXParserFactory.newInstance();
SAXParser parser = factory.newSAXParser();
SAX2 handler = new SAX2();
parser.parse(bis, handler);
parser = factory.newSAXParser();
} catch ( ParserConfigurationException | SAXException | IOException e ) {
System.out.println( e.getMessage() );
}
for (long[] val : mapCoastlines.values()) {
float lon0 = (val[0]>>32) / 10000000.0f;
float lat0 = (val[0]&0xffffffff) / 10000000.0f;
PointF p0 = new PointF(lon0, lat0);
float lon1 = (val[val.length-1]>>32) / 10000000.0f;
float lat1 = (val[val.length-1]&0xffffffff) / 10000000.0f;
PointF p1 = new PointF(lon1, lat1);
if (PointF.PointInJapan(p0) && PointF.PointInJapan(p1)) {
System.out.printf("head: %d, tail: %d\n", val[0], val[val.length-1]);
}
}
System.out.printf("実行時間＝%.1f秒", (System.currentTimeMillis() - start)/1000.0);
}
}
最終的には、カスタマイズしたOSM地図向けのバイナリレコード形式でファイル出力することになろうが、 まずはパフォーマンスよりも分かりやすさを重視して、テキストファイル出力とする。
現状では、単純なポリゴンしか存在しないが、PostGISの出力に合わせて、MULTIPOLYGON としておこう。
List<long[]> listPolygons から long[] nids を取り出して、ノードID を座標値に変換して 出力する。
海岸線のノードIDから経度・緯度を得るには HashMap<Long, Long> mapId2LonLat を使う。キーはノードID で、 値は整数化した経度と緯度を合わせたものである。
プログラムとしては、末尾に以下を追加した。
for (long[] nids : listPolygons) {
System.out.printf("MULTIPOLYGON(((");
for (int i = 0; i < nids.length; i++) {
long nid = nids[i];
long lonlat = mapId2LonLat.get(nid);
double lon = (lonlat>>32)/10000000.0;
double lat = (lonlat&0xffffffffL)/10000000.0;
System.out.printf("%s%.7f %.7f", (i>0 ? "," : ""), lon, lat);
}
System.out.printf(")))\n");
}
三宅島のデータに対して実行した結果を以下に示す。 今回のプログラムで生成するデータは全レコードとも単独の閉ループである。従って、 途中に、「(」や「)」は現れない。しかし、一般には、一番外側の（ ）が全体を括るもので、 次の ( ) が複数の outer polygon を括るもので、最後の （ ）が個々の outer または inner polygon を括るものである。
MULTIPOLYGON(((139.4773730 34.0619460,139.4774620 34.0619880,139.4774770 34.0620360,139.4774750 34.0621040,139.4774520 34.0621390,139.4773570 34.0621620,139.4772850 34.0620810,139.4772050 34.0620369,139.4772150 34.0620010,139.4773030 34.0619510,139.4773730 34.0619460))) MULTIPOLYGON(((139.4759690 34.0800760,139.4760780 34.0800740,139.4761840 34.0802520,139.4761650 34.0802970,139.4760330 34.0803920,139.4759860 34.0804040,139.4759320 34.0804090,139.4759010 34.0803810,139.4758500 34.0802810,139.4758440 34.0801440,139.4759690 34.0800760))) [以下略]
wayオブジェクト natural=coastline については以下のようにして、Overpass API[4] で入手できた[2023.1.22]。 ２GB近いが、９分弱で得られた。
日本全体の osmファイルは 40GB近いので、それに比べれば、5％程度に過ぎないので、 パソコンで解析できる。
世界地図は低ズームだけでよいので、簡略化が必須となる。面積の小さいポリゴンは無視する。 大きいポリゴンについてはノードを間引く。 間引きでは、自己交差が生じると、ポリゴンの塗りつぶしにエラーが生じる可能性があるため、要注意である。 描画エラーが生じていないか、または、間引きで自己交差を回避する必要がある。
c:\gisa>curl --globoff -o c:/osm/world_coastline.xml -d @c:/osm/overpass_coastline.xml http://overpass-api.de/api/interpreter
% Total % Received % Xferd Average Speed Time Time Time Current
Dload Upload Total Spent Left Speed
100 1953M 0 1953M 0 176 3764k 0 --:--:-- 0:08:51 --:--:-- 4379k
<osm-script timeout="900" element-limit="2150000000"><!-- 60 minutes、約 2GB -->
<query type="way">
<has-kv k="natural" v="coastline"/>
</query>
<print/>
</osm-script>
element-limit はファイルサイズの制限ではないので、もっと小さい値でよい。
このコマンドでは node の座標が得られていない。 関連ノード情報を一緒に得ようとした場合、
日本地図では問題はないが、 世界全体を対象とすると、メモリ不足となる。
開発したプログラムで調べた結果、海岸線データを構成するノードは 71,131,960 となった。 ノードIDを指定してノード情報を得る場合、クエリのサイズが大きくなるため、分割がいるであろうが、 プログラム的には楽であろう。
curl --globoff -o output.xml -d @c:/osm/query.txt http://overpass-api.de/api/interpreter
% Total % Received % Xferd Average Speed Time Time Time Current
Dload Upload Total Spent Left Speed
100 2456 0 2097 100 359 1369 234 0:00:01 0:00:01 --:--:-- 1604
[query.txt]
node(id:4198135984,4198139995,4198140001,4198139997, [中略] ,4198135984);out;
得られた結果は、query.txt の順番ではなく、id の昇順となっている。 リクエストで id が重複していた場合、出力はひとつである。
<?xml version="1.0" encoding="UTF-8"?> <osm version="0.6" generator="Overpass API 0.7.59 e21c39fe"> <note>The data included in this document is from www.openstreetmap.org. The data is made available under ODbL.</note> <meta osm_base="2023-01-23T02:10:39Z"/> <node id="4198135984" lat="35.5699406" lon="139.5265572"/> <node id="4198139995" lat="35.5700412" lon="139.5267127"/> <node id="4198139997" lat="35.5700618" lon="139.5266279"/> <node id="4198140001" lat="35.5700837" lon="139.5266609"/> <node id="4198140010" lat="35.5701634" lon="139.5269836"/> [中略] <node id="4198140075" lat="35.5705261" lon="139.5268032"/> </osm>