diff --git a/.depend b/.depend index 9483c71..5f7ddd0 100644 --- a/.depend +++ b/.depend @@ -1,168 +1,190 @@ -./database/column.cmo: generic/util.cmi generic/mapped_file.cmi \ +./database/column.cmo : generic/util.cmi generic/mapped_file.cmi \ ./database/column.cmi -./database/column.cmx: generic/util.cmx generic/mapped_file.cmx \ +./database/column.cmx : generic/util.cmx generic/mapped_file.cmx \ ./database/column.cmi -./database/column.cmi: -./database/column_ops.cmo: database/column.cmi ./database/column_ops.cmi -./database/column_ops.cmx: database/column.cmx ./database/column_ops.cmi -./database/column_ops.cmi: database/column.cmi -./database/dictionary.cmo: generic/util.cmi database/sorting.cmi \ +./database/column.cmi : +./database/column_ops.cmo : database/column.cmi ./database/column_ops.cmi +./database/column_ops.cmx : database/column.cmx ./database/column_ops.cmi +./database/column_ops.cmi : database/column.cmi +./database/dictionary.cmo : generic/util.cmi database/sorting.cmi \ generic/mapped_file.cmi database/column.cmi generic/bytearray.cmi \ ./database/dictionary.cmi -./database/dictionary.cmx: generic/util.cmx database/sorting.cmx \ +./database/dictionary.cmx : generic/util.cmx database/sorting.cmx \ generic/mapped_file.cmx database/column.cmx generic/bytearray.cmx \ ./database/dictionary.cmi -./database/dictionary.cmi: database/column.cmi -./database/join.cmo: database/column.cmi ./database/join.cmi -./database/join.cmx: database/column.cmx ./database/join.cmi -./database/join.cmi: database/column.cmi -./database/projection.cmo: database/column.cmi ./database/projection.cmi -./database/projection.cmx: database/column.cmx ./database/projection.cmi -./database/projection.cmi: database/column.cmi -./database/query.cmo: database/dictionary.cmi generic/debug.cmi \ - database/column.cmi -./database/query.cmx: database/dictionary.cmx generic/debug.cmx \ - database/column.cmx -./database/rtree.cmo: generic/util.cmi generic/pqueue.cmi \ +./database/dictionary.cmi : database/column.cmi +./database/join.cmo : database/column.cmi ./database/join.cmi +./database/join.cmx : database/column.cmx ./database/join.cmi +./database/join.cmi : database/column.cmi +./database/projection.cmo : database/column.cmi ./database/projection.cmi +./database/projection.cmx : database/column.cmx ./database/projection.cmi +./database/projection.cmi : database/column.cmi +./database/query.cmo : generic/util.cmi database/dictionary.cmi \ + generic/debug.cmi database/column.cmi +./database/query.cmx : generic/util.cmx database/dictionary.cmx \ + generic/debug.cmx database/column.cmx +./database/rtree.cmo : generic/util.cmi generic/pqueue.cmi \ generic/mapped_file.cmi database/column.cmi ./database/rtree.cmi -./database/rtree.cmx: generic/util.cmx generic/pqueue.cmx \ +./database/rtree.cmx : generic/util.cmx generic/pqueue.cmx \ generic/mapped_file.cmx database/column.cmx ./database/rtree.cmi -./database/rtree.cmi: -./database/sorting.cmo: database/column.cmi ./database/sorting.cmi -./database/sorting.cmx: database/column.cmx ./database/sorting.cmi -./database/sorting.cmi: database/column.cmi -./database/table.cmo: ./database/table.cmi -./database/table.cmx: ./database/table.cmi -./database/table.cmi: -./generic/binary_heap.cmo: ./generic/binary_heap.cmi -./generic/binary_heap.cmx: ./generic/binary_heap.cmi -./generic/binary_heap.cmi: -./generic/bitvect.cmo: ./generic/bitvect.cmi -./generic/bitvect.cmx: ./generic/bitvect.cmi -./generic/bitvect.cmi: -./generic/bytearray.cmo: ./generic/bytearray.cmi -./generic/bytearray.cmx: ./generic/bytearray.cmi -./generic/bytearray.cmi: -./generic/data_stream.cmo: ./generic/data_stream.cmi -./generic/data_stream.cmx: ./generic/data_stream.cmi -./generic/data_stream.cmi: -./generic/debug.cmo: ./generic/debug.cmi -./generic/debug.cmx: ./generic/debug.cmi -./generic/debug.cmi: -./generic/lru_cache.cmo: ./generic/lru_cache.cmi -./generic/lru_cache.cmx: ./generic/lru_cache.cmi -./generic/lru_cache.cmi: -./generic/mapped_file.cmo: generic/debug.cmi ./generic/mapped_file.cmi -./generic/mapped_file.cmx: generic/debug.cmx ./generic/mapped_file.cmi -./generic/mapped_file.cmi: -./generic/pqueue.cmo: ./generic/pqueue.cmi -./generic/pqueue.cmx: ./generic/pqueue.cmi -./generic/pqueue.cmi: -./generic/protobuf.cmo: ./generic/protobuf.cmi -./generic/protobuf.cmx: ./generic/protobuf.cmi -./generic/protobuf.cmi: -./generic/task.cmo: generic/util.cmi generic/debug.cmi generic/bytearray.cmi \ - ./generic/task.cmi -./generic/task.cmx: generic/util.cmx generic/debug.cmx generic/bytearray.cmx \ - ./generic/task.cmi -./generic/task.cmi: -./generic/util.cmo: ./generic/util.cmi -./generic/util.cmx: ./generic/util.cmi -./generic/util.cmi: -./not_used/bloom.cmo: generic/bitvect.cmi ./not_used/bloom.cmi -./not_used/bloom.cmx: generic/bitvect.cmx ./not_used/bloom.cmi -./not_used/bloom.cmi: -./not_used/hierholzer.cmo: database/rtree.cmi database/column.cmi \ +./database/rtree.cmi : +./database/sorting.cmo : database/column.cmi ./database/sorting.cmi +./database/sorting.cmx : database/column.cmx ./database/sorting.cmi +./database/sorting.cmi : database/column.cmi +./database/table.cmo : ./database/table.cmi +./database/table.cmx : ./database/table.cmi +./database/table.cmi : +./generic/binary_heap.cmo : ./generic/binary_heap.cmi +./generic/binary_heap.cmx : ./generic/binary_heap.cmi +./generic/binary_heap.cmi : +./generic/bitvect.cmo : ./generic/bitvect.cmi +./generic/bitvect.cmx : ./generic/bitvect.cmi +./generic/bitvect.cmi : +./generic/bytearray.cmo : ./generic/bytearray.cmi +./generic/bytearray.cmx : ./generic/bytearray.cmi +./generic/bytearray.cmi : +./generic/data_stream.cmo : ./generic/data_stream.cmi +./generic/data_stream.cmx : ./generic/data_stream.cmi +./generic/data_stream.cmi : +./generic/debug.cmo : ./generic/debug.cmi +./generic/debug.cmx : ./generic/debug.cmi +./generic/debug.cmi : +./generic/lru_cache.cmo : ./generic/lru_cache.cmi +./generic/lru_cache.cmx : ./generic/lru_cache.cmi +./generic/lru_cache.cmi : +./generic/mapped_file.cmo : generic/debug.cmi ./generic/mapped_file.cmi +./generic/mapped_file.cmx : generic/debug.cmx ./generic/mapped_file.cmi +./generic/mapped_file.cmi : +./generic/pqueue.cmo : ./generic/pqueue.cmi +./generic/pqueue.cmx : ./generic/pqueue.cmi +./generic/pqueue.cmi : +./generic/protobuf.cmo : ./generic/protobuf.cmi +./generic/protobuf.cmx : ./generic/protobuf.cmi +./generic/protobuf.cmi : +./generic/task.cmo : generic/util.cmi generic/debug.cmi \ + generic/bytearray.cmi ./generic/task.cmi +./generic/task.cmx : generic/util.cmx generic/debug.cmx \ + generic/bytearray.cmx ./generic/task.cmi +./generic/task.cmi : +./generic/util.cmo : ./generic/util.cmi +./generic/util.cmx : ./generic/util.cmi +./generic/util.cmi : +./not_used/bloom.cmo : generic/bitvect.cmi ./not_used/bloom.cmi +./not_used/bloom.cmx : generic/bitvect.cmx ./not_used/bloom.cmi +./not_used/bloom.cmi : +./not_used/hierholzer.cmo : database/rtree.cmi database/column.cmi \ ./not_used/hierholzer.cmi -./not_used/hierholzer.cmx: database/rtree.cmx database/column.cmx \ +./not_used/hierholzer.cmx : database/rtree.cmx database/column.cmx \ ./not_used/hierholzer.cmi -./not_used/hierholzer.cmi: -./not_used/huffman.cmo: generic/pqueue.cmi -./not_used/huffman.cmx: generic/pqueue.cmx -./not_used/osm.cmi: -./not_used/system.cmo: ./not_used/system.cmi -./not_used/system.cmx: ./not_used/system.cmi -./not_used/system.cmi: -./osm/category.cmo: database/dictionary.cmi ./osm/category.cmi -./osm/category.cmx: database/dictionary.cmx ./osm/category.cmi -./osm/category.cmi: database/dictionary.cmi -./osm/clipping.cmo: ./osm/clipping.cmi -./osm/clipping.cmx: ./osm/clipping.cmi -./osm/clipping.cmi: -./osm/coastline.cmo: generic/util.cmi database/rtree.cmi osm/geometry.cmi \ - osm/douglas_peucker.cmi database/column.cmi osm/clipping.cmi -./osm/coastline.cmx: generic/util.cmx database/rtree.cmx osm/geometry.cmx \ - osm/douglas_peucker.cmx database/column.cmx osm/clipping.cmx -./osm/contraction.cmo: generic/util.cmi database/sorting.cmi \ - generic/mapped_file.cmi database/column_ops.cmi database/column.cmi \ - generic/bitvect.cmi generic/binary_heap.cmi -./osm/contraction.cmx: generic/util.cmx database/sorting.cmx \ - generic/mapped_file.cmx database/column_ops.cmx database/column.cmx \ - generic/bitvect.cmx generic/binary_heap.cmx -./osm/display.cmo: database/rtree.cmi osm/routing.cmi generic/lru_cache.cmi \ - osm/line_smoothing.cmi osm/geometry.cmi osm/douglas_peucker.cmi \ - database/column.cmi osm/category.cmi -./osm/display.cmx: database/rtree.cmx osm/routing.cmx generic/lru_cache.cmx \ - osm/line_smoothing.cmx osm/geometry.cmx osm/douglas_peucker.cmx \ - database/column.cmx osm/category.cmx -./osm/douglas_peucker.cmo: ./osm/douglas_peucker.cmi -./osm/douglas_peucker.cmx: ./osm/douglas_peucker.cmi -./osm/douglas_peucker.cmi: -./osm/geometry.cmo: ./osm/geometry.cmi -./osm/geometry.cmx: ./osm/geometry.cmi -./osm/geometry.cmi: -./osm/highway.cmo: generic/util.cmi database/sorting.cmi \ - osm/routing_profile.cmi database/projection.cmi generic/mapped_file.cmi \ - database/join.cmi osm/geometry.cmi database/dictionary.cmi \ +./not_used/hierholzer.cmi : +./not_used/hilbert.cmo : +./not_used/hilbert.cmx : +./not_used/huffman.cmo : generic/pqueue.cmi +./not_used/huffman.cmx : generic/pqueue.cmx +./not_used/osm.cmi : +./not_used/system.cmo : ./not_used/system.cmi +./not_used/system.cmx : ./not_used/system.cmi +./not_used/system.cmi : +./osm/coastline.cmo : generic/util.cmi osm/lib/osm_coastline.cmo \ + database/column.cmi +./osm/coastline.cmx : generic/util.cmx osm/lib/osm_coastline.cmx \ + database/column.cmx +./osm/contraction.cmo : database/sorting.cmi osm/lib/osm_contraction.cmo \ + database/column_ops.cmi database/column.cmi generic/bitvect.cmi +./osm/contraction.cmx : database/sorting.cmx osm/lib/osm_contraction.cmx \ + database/column_ops.cmx database/column.cmx generic/bitvect.cmx +./osm/display.cmo : database/rtree.cmi osm/lib/osm_geometry.cmi \ + osm/lib/osm_display.cmo database/column.cmi +./osm/display.cmx : database/rtree.cmx osm/lib/osm_geometry.cmx \ + osm/lib/osm_display.cmx database/column.cmx +./osm/highway.cmo : generic/util.cmi database/sorting.cmi \ + osm/routing_profile.cmi database/projection.cmi osm/lib/osm_geometry.cmi \ + generic/mapped_file.cmi database/join.cmi database/dictionary.cmi \ database/column_ops.cmi database/column.cmi -./osm/highway.cmx: generic/util.cmx database/sorting.cmx \ - osm/routing_profile.cmx database/projection.cmx generic/mapped_file.cmx \ - database/join.cmx osm/geometry.cmx database/dictionary.cmx \ +./osm/highway.cmx : generic/util.cmx database/sorting.cmx \ + osm/routing_profile.cmx database/projection.cmx osm/lib/osm_geometry.cmx \ + generic/mapped_file.cmx database/join.cmx database/dictionary.cmx \ database/column_ops.cmx database/column.cmx -./osm/line_smoothing.cmo: ./osm/line_smoothing.cmi -./osm/line_smoothing.cmx: ./osm/line_smoothing.cmi -./osm/line_smoothing.cmi: -./osm/linear.cmo: database/sorting.cmi database/rtree.cmi \ - database/projection.cmi database/join.cmi osm/geometry.cmi \ - database/dictionary.cmi database/column_ops.cmi database/column.cmi \ - osm/category.cmi -./osm/linear.cmx: database/sorting.cmx database/rtree.cmx \ - database/projection.cmx database/join.cmx osm/geometry.cmx \ - database/dictionary.cmx database/column_ops.cmx database/column.cmx \ - osm/category.cmx -./osm/multipolygons.cmo: generic/util.cmi database/sorting.cmi \ - database/projection.cmi database/join.cmi osm/geometry.cmi \ +./osm/lib/osm_category.cmo : database/dictionary.cmi \ + ./osm/lib/osm_category.cmi +./osm/lib/osm_category.cmx : database/dictionary.cmx \ + ./osm/lib/osm_category.cmi +./osm/lib/osm_category.cmi : database/dictionary.cmi +./osm/lib/osm_clipping.cmo : ./osm/lib/osm_clipping.cmi +./osm/lib/osm_clipping.cmx : ./osm/lib/osm_clipping.cmi +./osm/lib/osm_clipping.cmi : +./osm/lib/osm_coastline.cmo : generic/util.cmi database/rtree.cmi \ + osm/lib/osm_geometry.cmi osm/lib/osm_douglas_peucker.cmi \ + osm/lib/osm_clipping.cmi database/column.cmi +./osm/lib/osm_coastline.cmx : generic/util.cmx database/rtree.cmx \ + osm/lib/osm_geometry.cmx osm/lib/osm_douglas_peucker.cmx \ + osm/lib/osm_clipping.cmx database/column.cmx +./osm/lib/osm_contraction.cmo : database/column.cmi generic/bitvect.cmi \ + generic/binary_heap.cmi +./osm/lib/osm_contraction.cmx : database/column.cmx generic/bitvect.cmx \ + generic/binary_heap.cmx +./osm/lib/osm_display.cmo : database/rtree.cmi osm/routing.cmi \ + osm/lib/osm_geometry.cmi osm/lib/osm_douglas_peucker.cmi \ + osm/lib/osm_category.cmi generic/lru_cache.cmi osm/line_smoothing.cmi \ + database/column.cmi +./osm/lib/osm_display.cmx : database/rtree.cmx osm/routing.cmx \ + osm/lib/osm_geometry.cmx osm/lib/osm_douglas_peucker.cmx \ + osm/lib/osm_category.cmx generic/lru_cache.cmx osm/line_smoothing.cmx \ + database/column.cmx +./osm/lib/osm_douglas_peucker.cmo : ./osm/lib/osm_douglas_peucker.cmi +./osm/lib/osm_douglas_peucker.cmx : ./osm/lib/osm_douglas_peucker.cmi +./osm/lib/osm_douglas_peucker.cmi : +./osm/lib/osm_geometry.cmo : ./osm/lib/osm_geometry.cmi +./osm/lib/osm_geometry.cmx : ./osm/lib/osm_geometry.cmi +./osm/lib/osm_geometry.cmi : +./osm/line_smoothing.cmo : ./osm/line_smoothing.cmi +./osm/line_smoothing.cmx : ./osm/line_smoothing.cmi +./osm/line_smoothing.cmi : +./osm/linear.cmo : generic/util.cmi database/sorting.cmi database/rtree.cmi \ + database/projection.cmi osm/lib/osm_geometry.cmi osm/lib/osm_category.cmi \ + database/join.cmi database/dictionary.cmi database/column_ops.cmi \ + database/column.cmi +./osm/linear.cmx : generic/util.cmx database/sorting.cmx database/rtree.cmx \ + database/projection.cmx osm/lib/osm_geometry.cmx osm/lib/osm_category.cmx \ + database/join.cmx database/dictionary.cmx database/column_ops.cmx \ + database/column.cmx +./osm/multipolygons.cmo : generic/util.cmi database/sorting.cmi \ + database/projection.cmi osm/lib/osm_geometry.cmi database/join.cmi \ database/dictionary.cmi generic/debug.cmi generic/data_stream.cmi \ database/column_ops.cmi database/column.cmi -./osm/multipolygons.cmx: generic/util.cmx database/sorting.cmx \ - database/projection.cmx database/join.cmx osm/geometry.cmx \ +./osm/multipolygons.cmx : generic/util.cmx database/sorting.cmx \ + database/projection.cmx osm/lib/osm_geometry.cmx database/join.cmx \ database/dictionary.cmx generic/debug.cmx generic/data_stream.cmx \ database/column_ops.cmx database/column.cmx -./osm/parser.cmo: generic/util.cmi generic/task.cmi generic/protobuf.cmi \ +./osm/parser.cmo : generic/util.cmi generic/task.cmi generic/protobuf.cmi \ generic/debug.cmi database/column.cmi -./osm/parser.cmx: generic/util.cmx generic/task.cmx generic/protobuf.cmx \ +./osm/parser.cmx : generic/util.cmx generic/task.cmx generic/protobuf.cmx \ generic/debug.cmx database/column.cmx -./osm/prepare.cmo: generic/util.cmi database/sorting.cmi \ +./osm/prepare.cmo : generic/util.cmi database/sorting.cmi \ database/projection.cmi database/join.cmi database/dictionary.cmi \ database/column_ops.cmi database/column.cmi -./osm/prepare.cmx: generic/util.cmx database/sorting.cmx \ +./osm/prepare.cmx : generic/util.cmx database/sorting.cmx \ database/projection.cmx database/join.cmx database/dictionary.cmx \ database/column_ops.cmx database/column.cmx -./osm/profile_car.cmo: database/table.cmi osm/routing_profile.cmi -./osm/profile_car.cmx: database/table.cmx osm/routing_profile.cmx -./osm/profile_pedestrian.cmo: database/table.cmi osm/routing_profile.cmi -./osm/profile_pedestrian.cmx: database/table.cmx osm/routing_profile.cmx -./osm/routing.cmo: generic/pqueue.cmi database/column.cmi ./osm/routing.cmi -./osm/routing.cmx: generic/pqueue.cmx database/column.cmx ./osm/routing.cmi -./osm/routing.cmi: -./osm/routing_profile.cmo: database/dictionary.cmi ./osm/routing_profile.cmi -./osm/routing_profile.cmx: database/dictionary.cmx ./osm/routing_profile.cmi -./osm/routing_profile.cmi: database/dictionary.cmi -./osm/surfaces.cmo: generic/util.cmi database/sorting.cmi database/rtree.cmi \ - database/projection.cmi database/join.cmi osm/geometry.cmi \ - osm/douglas_peucker.cmi database/dictionary.cmi generic/data_stream.cmi \ - database/column_ops.cmi database/column.cmi osm/category.cmi -./osm/surfaces.cmx: generic/util.cmx database/sorting.cmx database/rtree.cmx \ - database/projection.cmx database/join.cmx osm/geometry.cmx \ - osm/douglas_peucker.cmx database/dictionary.cmx generic/data_stream.cmx \ - database/column_ops.cmx database/column.cmx osm/category.cmx +./osm/profile_car.cmo : database/table.cmi osm/routing_profile.cmi +./osm/profile_car.cmx : database/table.cmx osm/routing_profile.cmx +./osm/profile_pedestrian.cmo : database/table.cmi osm/routing_profile.cmi +./osm/profile_pedestrian.cmx : database/table.cmx osm/routing_profile.cmx +./osm/routing.cmo : generic/pqueue.cmi database/column.cmi ./osm/routing.cmi +./osm/routing.cmx : generic/pqueue.cmx database/column.cmx ./osm/routing.cmi +./osm/routing.cmi : +./osm/routing_profile.cmo : database/dictionary.cmi \ + ./osm/routing_profile.cmi +./osm/routing_profile.cmx : database/dictionary.cmx \ + ./osm/routing_profile.cmi +./osm/routing_profile.cmi : database/dictionary.cmi +./osm/surfaces.cmo : generic/util.cmi database/sorting.cmi \ + database/rtree.cmi database/projection.cmi osm/lib/osm_geometry.cmi \ + osm/lib/osm_douglas_peucker.cmi osm/lib/osm_category.cmi \ + database/join.cmi database/dictionary.cmi generic/data_stream.cmi \ + database/column_ops.cmi database/column.cmi +./osm/surfaces.cmx : generic/util.cmx database/sorting.cmx \ + database/rtree.cmx database/projection.cmx osm/lib/osm_geometry.cmx \ + osm/lib/osm_douglas_peucker.cmx osm/lib/osm_category.cmx \ + database/join.cmx database/dictionary.cmx generic/data_stream.cmx \ + database/column_ops.cmx database/column.cmx diff --git a/Makefile b/Makefile index 01b740e..6c53ab0 100644 --- a/Makefile +++ b/Makefile @@ -23,12 +23,12 @@ DATABASE=\ projection.cmx table.cmx dictionary.cmx column_ops.cmx rtree.cmx OSM=\ - geometry.cmx + lib/osm_geometry.cmx ROUTING=\ routing_profile.cmx profile_car.cmx profile_pedestrian.cmx -DIRS=-I generic -I database -I osm +DIRS=-I generic -I database -I osm -I osm/lib OBJS= $(addprefix generic/,$(GENERIC)) $(addprefix database/,$(DATABASE)) \ $(addprefix osm/,$(OSM)) @@ -46,22 +46,22 @@ load: $(OBJS) osm/parser.cmx osm/prepare.cmx multipolygons: $(OBJS) osm/multipolygons.cmx $(OCAMLOPT) $(OPTLINKFLAGS) -o $@ $^ -linear: $(OBJS) osm/category.cmx osm/linear.cmx +linear: $(OBJS) osm/lib/osm_category.cmx osm/linear.cmx $(OCAMLOPT) $(OPTLINKFLAGS) -o $@ $^ -surfaces: $(OBJS) osm/douglas_peucker.cmx osm/category.cmx osm/surfaces.cmx +surfaces: $(OBJS) osm/lib/osm_douglas_peucker.cmx osm/lib/osm_category.cmx osm/surfaces.cmx $(OCAMLOPT) $(OPTLINKFLAGS) -o $@ $^ highway: $(OBJS) $(addprefix osm/, $(ROUTING)) osm/highway.cmx $(OCAMLOPT) $(OPTLINKFLAGS) -o $@ $^ -contraction: $(OBJS) osm/contraction.cmx +contraction: $(OBJS) osm/lib/osm_contraction.cmx osm/contraction.cmx $(OCAMLOPT) $(OPTLINKFLAGS) -o $@ $^ -display: $(OBJS) osm/routing.cmx osm/line_smoothing.cmx osm/douglas_peucker.cmx osm/category.cmx osm/display.cmx +display: $(OBJS) osm/routing.cmx osm/line_smoothing.cmx osm/lib/osm_douglas_peucker.cmx osm/lib/osm_category.cmx osm/lib/osm_display.cmx osm/display.cmx $(OCAMLOPT) $(OPTLINKFLAGS) -o $@ $^ -coastline: $(OBJS) osm/category.cmx osm/douglas_peucker.cmx osm/clipping.cmx osm/coastline.cmx +coastline: $(OBJS) osm/lib/osm_category.cmx osm/lib/osm_douglas_peucker.cmx osm/lib/osm_clipping.cmx osm/lib/osm_coastline.cmx osm/coastline.cmx $(OCAMLOPT) $(OPTLINKFLAGS) -o $@ $^ clean:: diff --git a/osm/coastline.ml b/osm/coastline.ml index 6cad462..94e3ec5 100644 --- a/osm/coastline.ml +++ b/osm/coastline.ml @@ -16,296 +16,11 @@ * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *) -(* -TODO -==== -- simplification should depend on latitude - (should be done after projection) -*) - -let leaf_size = 2048 -let max_polygon_size = 200 - -(****) - let _ = Printexc.record_backtrace true let _ = Column.set_database "/tmp/osm" (****) -let read_big_int ch = - let i3 = input_byte ch in - let i2 = input_byte ch in - let i1 = input_byte ch in - let i0 = input_byte ch in - i0 lor (i1 lsl 8) lor (i2 lsl 16) lor (i3 lsl 24) - -let read_lit_int ch = - let i0 = input_byte ch in - let i1 = input_byte ch in - let i2 = input_byte ch in - let i3 = input_byte ch in - i0 lor (i1 lsl 8) lor (i2 lsl 16) lor (i3 lsl 24) - -let read_float ch = - let i0 = Int64.of_int (read_lit_int ch) in - let i1 = Int64.of_int (read_lit_int ch) in - Int64.float_of_bits (Int64.logor i0 (Int64.shift_left i1 32)) - -let build_rings l = - let (lst, rem) = - List.partition - (fun (x, y) -> - let l = Array.length x in x.(0) = x.(l - 1) && y.(0) = y.(l - 1)) - l - in - let succ = Hashtbl.create 17 in - let pred = Hashtbl.create 17 in - let rem = Array.of_list rem in - Array.iteri - (fun i (x, y) -> - let l = Array.length x in - Hashtbl.add succ (x.(0), y.(0)) i; - Hashtbl.add pred (x.(l - 1), y.(l - 1)) i) - rem; - let lst = ref lst in - let rec follow i l = - let (x, y) = rem.(i) in - let len = Array.length x in - try - let k = Hashtbl.find pred (x.(0), y.(0)) in - let l = (Array.sub x 1 (len - 1), Array.sub y 1 (len - 1)) :: l in - follow k l - with Not_found -> - (x, y) :: l - in - for i = 0 to Array.length rem - 1 do - let (x, y) = rem.(i) in - let len = Array.length x in - if not (Hashtbl.mem succ (x.(len - 1), y.(len - 1))) then begin - let l = follow i [] in - let x = Array.concat (List.map fst l) in - let y = Array.concat (List.map snd l) in - let len = Array.length x in - let (x, y) = - if x.(0) < -1799999000. && x.(len - 1) < -1799999000. then begin - x.(0) <- -1800000000.; x.(len - 1) <- -1800000000.; - (Array.append x [|x.(0)|], Array.append y [|y.(0)|]) - end else if - x.(0) > 1799999000. && x.(len - 1) > 1799999000. - then begin - x.(0) <- 1800000000.; x.(len - 1) <- 1800000000.; - (Array.append x [|x.(0)|], Array.append y [|y.(0)|]) - end else if - x.(0) = -1800000000. && x.(len -1) = 1800000000. - then begin - (* Antartica *) - (Array.append x [|1800000000.; -1800000000.; x.(0)|], - Array.append y [| 850000000.; 850000000.; y.(0)|]) - end else begin -Format.eprintf "%f %f - %f %f@." x.(0) y.(0) x.(len - 1) y.(len - 1); - assert false - end - in - lst := (x, y) :: !lst - end - done; - !lst - -let read_record ch = - let num = read_big_int ch in - let len = read_big_int ch in - ignore (num, len); - let typ = read_lit_int ch in - assert (typ = 3); (* polyline *) - - for i = 0 to 3 do ignore (read_float ch) done; (* bbox *) - - let num_parts = read_lit_int ch in - assert (num_parts = 1); - - let num_points = read_lit_int ch in - - ignore (read_lit_int ch); (* first and unique part *) - - let lon = Array.make (num_points - 1) 0. in - let lat = Array.make (num_points - 1) 0. in - - for i = 0 to num_points - 2 do - lon.(i) <- 10_000_000. *. read_float ch; - lat.(i) <- 10_000_000. *. read_float ch - done; - let lon0 = 10_000_000. *. read_float ch in - let lat0 = 10_000_000. *. read_float ch in - ((lon, lat), (lon0, lat0)) - -let open_file () = - let ch = open_in Sys.argv.(1) in - let magic = read_big_int ch in - if magic = 9994 then begin - seek_in ch 0; - ch - end else if magic = 0x504b0304 then begin - close_in ch; - let z = Zip.open_in Sys.argv.(1) in - let e = - try - Zip.find_entry z "coastlines-split-4326/lines.shp" - with Not_found -> - Util.fail "Zip entry 'coastlines-split-4326/lines.shp' not found" - in - let (r, w) = Unix.pipe () in - match Unix.fork () with - 0 -> - Unix.close r; - let ch = Unix.out_channel_of_descr w in - Zip.copy_entry_to_channel z e ch; - flush ch; - exit 0 - | _ -> - Unix.close w; - Unix.in_channel_of_descr r - end else - Util.fail (Format.sprintf "bad file format") - -let load () = - Format.eprintf "==== Loading.@."; - - let ch = open_file () in - - for i = 0 to 24 do - ignore (read_big_int ch) - done; - -(* - let a = read_float ch in - Format.eprintf "%f@." a; - let a = read_float ch in - Format.eprintf "%f@." a; - let a = read_float ch in - Format.eprintf "%f@." a; - let a = read_float ch in - Format.eprintf "%f@." a; -*) - - let lst = ref [] in - let cur = ref [] in - let lon0' = ref 0. in - let lat0' = ref 0. in - begin try - while true do - let ((lon, lat), (lon0, lat0)) = read_record ch in - if !cur <> [] && (lat.(0) <> !lat0' || lon.(0) <> !lon0') then begin - let path = List.rev (([|!lon0'|], [|!lat0'|]) :: !cur) in - let lon'' = Array.concat (List.map fst path) in - let lat'' = Array.concat (List.map snd path) in - lst := (lon'', lat'') :: !lst; - cur := [] - end; - cur := (lon, lat) :: !cur; - lon0' := lon0; - lat0' := lat0 - done - with End_of_file -> () end; - let polys = build_rings !lst in -(* -Format.eprintf "%d@." (List.length polys); -*) - polys - -(**** Quickselect ****) - -let swap (tbl : float array) i j = - let v = tbl.(i) in tbl.(i) <- tbl.(j); tbl.(j) <- v - -let partition (tbl : float array) left right pivot = - let p = tbl.(pivot) in - swap tbl pivot right; - let j = ref left in - for i = left to right - 1 do - if tbl.(i) < p then begin - swap tbl i !j; - incr j - end - done; - swap tbl right !j; - !j - -let rec select tbl left right k = - if left = right then - tbl.(left) - else begin - let pivot = Random.int (right - left + 1) + left in - let pivot = partition tbl left right pivot in - if pivot = k then - tbl.(pivot) - else if k < pivot then - select tbl left (pivot - 1) k - else - select tbl (pivot + 1) right k - end - -let median tbl = - let len = Array.length tbl in - select tbl 0 (len - 1) (len / 2) - -let poly_median polys = - let l = List.rev_map fst polys in (* Tail recursive! *) - median (Array.concat l) - -let stats polys = - let a = Array.concat (List.rev_map fst polys) in (* Tail recursive! *) - let min = ref max_float in - let max = ref (-. max_float) in - for i = 0 to Array.length a - 1 do - if a.(i) < !min then min := a.(i); - if a.(i) > !max then max := a.(i) - done; - let min' = ref max_float in - let max' = ref (-. max_float) in - for i = 0 to Array.length a - 1 do - if a.(i) < !min' && a.(i) > !min then min' := a.(i); - if a.(i) > !max' && a.(i) < !max then max' := a.(i) - done; - (!min', median a, !max') - -(****) - -let swap_coord polys = List.map (fun (x, y) -> (y, x)) polys - -let ring_bbox (x, y) : float * float * float * float = - let xmin = ref x.(0) in - let ymin = ref y.(0) in - let xmax = ref x.(0) in - let ymax = ref y.(0) in - for i = 1 to Array.length x - 1 do - if x.(i) < !xmin then xmin := x.(i); - if y.(i) < !ymin then ymin := y.(i); - if x.(i) > !xmax then xmax := x.(i); - if y.(i) > !ymax then ymax := y.(i) - done; - (!xmin, !ymin, !xmax, !ymax) - -let min' x y : float = if x < y then x else y -let max' x y : float = if x < y then y else x - -let bbox polys = - match polys with - [] -> assert false - | r :: rem -> - List.fold_left - (fun (xmin, ymin, xmax, ymax) r -> - let (xmin', ymin', xmax', ymax') = ring_bbox r in - (min' xmin xmin', min' ymin ymin', - max' xmax xmax', max' ymax ymax')) - (ring_bbox r) rem - -let size polys = List.fold_left (fun n (x, _) -> n + Array.length x) 0 polys -let size polys = - List.fold_left (fun n (x, _) -> max n (Array.length x)) 0 polys - -(****) - (* let surface = Cairo.SVG.create ~fname:"/tmp/bar.svg" ~width:360. ~height:180. @@ -314,224 +29,20 @@ let _ = Cairo.translate ctx 180. 90. let show = false *) -let d = false -let next_perc = ref 0.005 - -let rec split_rec dir t0 t1 t2 polys l = - let sz = size polys in -if d then Format.eprintf "size: %d (%d)@." sz (List.length polys); - if sz > max_polygon_size then begin - let (xmin, ymin, xmax, ymax) = bbox polys in -if d then Format.eprintf "%f %f - %f %f@." xmin ymin xmax ymax; - if xmax -. xmin < ymax -. ymin then - split_rec (not dir) t0 t1 t2 (swap_coord polys) l - else begin - let m' = ((xmax +. xmin) /. 2.) in - let (m1, m, m2) = stats polys in -if d then Format.eprintf "%f %f %f / %f@." m1 m m2 m'; - let m = - if m2 < m' then m2 else - if m1 > m' then m1 else - if m > (xmin +. m') /. 2. && m < (xmax +. m') /. 2. then - m - else - m' - in -if d then Format.eprintf "Cutting at %f@." m; - let (left, right) = Clipping.perform m polys in - let l = split_rec dir t0 t1 ((t1 +. t2) /. 2.) left l in - split_rec dir t0 ((t1 +. t2) /. 2.) t2 right l - end - end else begin -if t2 >= !next_perc then begin -Util.set_msg - (Format.sprintf "splitting polygons: %s %.0f%% eta %.0fs" - (Util.progress_bar t2) (t2 *. 100.) - ((1. -. t2) *. (Unix.gettimeofday () -. t0) /. t2)); -(*Format.eprintf "%.1f%%@." t2;*) -next_perc := !next_perc +. 0.005 -end; - let polys = if dir then polys else swap_coord polys in -(* -if show then begin -List.iter - (fun (x, y) -> - Cairo.move_to ctx (x.(0) /. 10_000_000.) (y.(0) /. 10_000_000.); - for i = 1 to Array.length x - 2 do - Cairo.line_to ctx (x.(i) /. 10_000_000.) (y.(i) /. 10_000_000.) - done; - Cairo.Path.close ctx) - polys; -Cairo.fill ctx; -if t2 > 0.2 then -(Cairo.Surface.finish surface; exit 0) -end; -*) - polys @ l - end - -(****) - -module Bbox = Rtree.Bbox - -let int_of_sint i = if i >= 0 then 2 * i else - 2 * i - 1 - -let rec write_varint a p v = - if v < 128 then begin - a.[p] <- Char.chr (v land 127); - p + 1 - end else begin - a.[p] <- Char.chr ((v land 127) + 128); - write_varint a (p + 1) (v lsr 7) - end - -let write_signed_varint a p v = write_varint a p (int_of_sint v) - -let output_int_2 ch v = - output_byte ch (v land 0xff); - output_byte ch (v lsr 8) - -let open_rtree name ratio = - let (nm, st) = Rtree.open_out name in - let ch = open_out (Column.file_in_database (Filename.concat name "ratio")) in - Printf.fprintf ch "%d\n" ratio; - close_out ch; - let ch = open_out nm in - let lengths = Array.make (leaf_size / 4) 0 in - let n = ref 0 in - let buf = String.make (2 * leaf_size) '\000' in - let pos = ref 0 in - let bbox = ref Bbox.empty in - let last_lat = ref 0 in - let last_lon = ref 0 in - let flush_ways () = - output_int_2 ch !n; - for i = 0 to !n - 1 do - output_int_2 ch lengths.(i); - done; - Rtree.append st !bbox; - output ch buf 0 (leaf_size - !n * 2 - 2); - n := 0; - pos := 0; - bbox := Bbox.empty; - last_lat := 0; - last_lon := 0 - in - let rec add_polygon (lon, lat) = - let pos' = ref !pos in - for i = 0 to Array.length lat - 2 do - pos' := write_signed_varint buf !pos' (lat.(i) - !last_lat); - last_lat := lat.(i); - pos' := write_signed_varint buf !pos' (lon.(i) - !last_lon); - last_lon := lon.(i) - done; - let pos' = !pos' in - if (!n + 1) * 2 + 2 + pos' > leaf_size then begin - assert (!n > 0); - flush_ways (); - add_polygon (lon, lat) - end else begin - lengths.(!n) <- Array.length lat - 1; - incr n; - pos := pos'; - let max x y : int = if x > y then x else y in - let min x y : int = if x > y then y else x in - let bbox' = - { Bbox. - min_lat = Array.fold_left min max_int lat; - max_lat = Array.fold_left max min_int lat; - min_lon = Array.fold_left min max_int lon; - max_lon = Array.fold_left max min_int lon } - in - bbox := Bbox.union !bbox bbox'; - if !n * 2 + 2 + pos' > leaf_size then flush_ways () - end - in - let close () = - if !n > 0 then flush_ways (); - Rtree.close_out st; - close_out ch - in - (add_polygon, close) - -(****) - -let simplify level polys = - Format.eprintf "==== Simplifying.@."; - (* We multiply by 4. to get good approximations near the poles *) - let scale = 256. /. 360. *. 2. ** level *. 4. in - let small_area = (10_000_000. /. scale) ** 2. in - let ratio = float (truncate (10_000_000. /. scale /. 2.)) in - (truncate ratio, - List.fold_left - (fun lst (lon, lat) -> - let (lon, lat) = Douglas_peucker.perform ratio lon lat in - if - abs_float (Geometry.polygon_area_float lon lat) <= small_area - then - lst - else - (lon, lat) :: lst) - [] polys) - -(****) - -let rescale_ring ratio (lon, lat) = - let l = Array.length lat in - let delta = ratio / 2 - 1 in - let lat' = Array.make l 0 in - let lon' = Array.make l 0 in - for i = 0 to l - 1 do - lat'.(i) <- (truncate lat.(i) + delta) / ratio; - lon'.(i) <- (truncate lon.(i) + delta) / ratio - done; - (lon', lat') - -let build_rtree name ratio polys = - Format.eprintf "==== Splitting.@."; - next_perc := 0.005; - let l = split_rec true (Unix.gettimeofday ()) 0. 1. polys [] in - Util.set_msg ""; -(* -Format.eprintf "==> %d@." (List.length l); -*) - Format.eprintf "==== Sorting.@."; - let a = Array.of_list l in - let a = - Array.map - (fun p -> - let (lon_min, lat_min, lon_max, lat_max) = ring_bbox p in - let lat = truncate ((lat_min +. lat_max) /. 2.) + 90_0000000 in - let lon = truncate ((lon_min +. lon_max) /. 2.) + 180_0000000 in - (Geometry.hilbert_coordinate lat lon, p)) - a - in - Array.sort (fun (x, _) (x', _) -> compare x x') a; - - Format.eprintf "==== Writing B-tree.@."; - let (add_polygon, close) = - open_rtree ("coastline/rtrees/" ^ name) (*"small"*) ratio in - Array.iter - (fun (_, (lon, lat)) -> add_polygon (rescale_ring ratio (lon, lat))) - a; - close () - -(****) let _ = if Array.length Sys.argv <> 2 then Util.fail (Format.sprintf "no file specified"); - let polys = load () in + let polys = Osm_coastline.load () in Gc.compact (); - build_rtree "small" 50 polys; + Osm_coastline.build_rtree "small" 50 polys; Gc.compact (); let build_rtree_simpl level name = - let (ratio, polys') = simplify level polys in - build_rtree name ratio polys'; + let (ratio, polys') = Osm_coastline.simplify level polys in + Osm_coastline.build_rtree name ratio polys'; in build_rtree_simpl 8. "8"; build_rtree_simpl 6. "6"; build_rtree_simpl 4. "4"; build_rtree_simpl 2. "2" - diff --git a/osm/contraction.ml b/osm/contraction.ml index 70ec5e2..486fa36 100644 --- a/osm/contraction.ml +++ b/osm/contraction.ml @@ -16,543 +16,12 @@ * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *) -(* -TODO: -- use memory mappings -- implement search -- oneway edges - - -- Dynamic graph: add/remove edges - -Dynamic graph -- bucket based -- compress buckets when mostly filled -- - - -Operations: -- reserve k edges -- remove an edge -- add an edge -- (remove all edges) - -Keep track of the usage of each page: ---> too small -> compact/merge ---> too large -> split - - -Page full --> balance with next page -Next page also full --> split page - - -Mapping page |-> set of node on this page - -Compress: -- copy to a temp location then blit back -- update nodes - -*) - let _ = Printexc.record_backtrace true - -let debug = ref false - -module Array1 = Bigarray.Array1 - let _ = Column.set_database "/tmp/osm" -let way_ids = Column.open_in (Column.named "highway" "way/id") -let way_id k = Column.get way_ids k - -type t = - { (* Nodes *) - first_edge : (int, Bigarray.int_elt, Bigarray.c_layout) Array1.t; - edge_count : (int, Bigarray.int_elt, Bigarray.c_layout) Array1.t; - priority : (float, Bigarray.float64_elt, Bigarray.c_layout) Array1.t; - depth : (int, Bigarray.int_elt, Bigarray.c_layout) Array1.t; - dist : (int, Bigarray.int_elt, Bigarray.c_layout) Array1.t; - remaining : (int, Bigarray.int_elt, Bigarray.c_layout) Array1.t; - node_count : int; - mutable remaining_count : int; - (* Edges *) - target : (int, Bigarray.int_elt, Bigarray.c_layout) Array1.t; - weight : (int, Bigarray.int_elt, Bigarray.c_layout) Array1.t; - flags : (int, Bigarray.int8_unsigned_elt, Bigarray.c_layout) Array1.t; - original_edges : (int, Bigarray.int_elt, Bigarray.c_layout) Array1.t; - shortcut_node : (int, Bigarray.int_elt, Bigarray.c_layout) Array1.t; - id : (int, Bigarray.int_elt, Bigarray.c_layout) Array1.t; - (* Output *) - out_source : Column.output_stream; - out_target : Column.output_stream; - out_weight : Column.output_stream; - out_flags : Column.output_stream; - out_shortcut : Column.output_stream } - -let first_edge st i = st.first_edge.{i} -let last_edge st i = st.first_edge.{i} + st.edge_count.{i} - 1 - -let forward_edge f = f land 2 <> 0 -let backward_edge f = f land 1 <> 0 -let reverse_edge f = ((f lsl 1) lor (f lsr 1)) land 3 - -let show_path st f i = - let rec follow i = - let d = st.dist.{i} in - if d > 0 then begin - let j = ref (-1) in - for k = first_edge st i to last_edge st i do - if - forward_edge st.flags.{k} && - st.dist.{st.target.{k}} + st.weight.{k} = d - then - j := k - done; -if (!j = -1) then begin - for k = first_edge st i to last_edge st i do - Format.eprintf "XXX %d: %d %d %d@." st.target.{k} st.dist.{st.target.{k}} st.weight.{k} d - done; -end; - Format.fprintf f "--%d(%d)--> %d" (way_id st.id.{!j}) st.weight.{!j} st.target.{!j}; - follow st.target.{!j} - end - in - Format.fprintf f "%d " i; - follow i - -let int_array len = Array1.create Bigarray.int Bigarray.c_layout len -let int8_array len = Array1.create Bigarray.int8_unsigned Bigarray.c_layout len - -let load () = - let c_nodes = Column.open_in (Column.named "highway" "node/idx") in - let c_in_node = Column.open_in (Column.named "highway" "edge/1bis") in - let c_out_node = Column.open_in (Column.named "highway" "edge/2bis") in - let c_weight = Column.open_in (Column.named "highway" "edge/weight_bis") in - let c_flags = Column.open_in (Column.named "highway" "edge/flags_bis") in - let c_id = Column.open_in (Column.named "highway" "edge/idx_bis") in - - let first_edge = int_array (Column.length c_nodes) in - let edge_count = int_array (Column.length c_nodes) in - let priority = - Array1.create Bigarray.float64 Bigarray.c_layout (Column.length c_nodes) in - let depth = int_array (Column.length c_nodes) in - let dist = int_array (Column.length c_nodes) in - let remaining = int_array (Column.length c_nodes) in - Array1.fill depth 0; - Array1.fill dist max_int; - for i = 0 to Array1.dim remaining - 1 do - remaining.{i} <- i - done; - - let len = Column.length c_in_node + Column.length c_nodes + 10 in - let target = int_array len in - let weight = int_array len in - let flags = int8_array len in - let original_edges = int_array len in - let shortcut_node = int_array len in - let id = int_array len in - Array1.fill original_edges 1; - - let s_in_node = Column.stream c_in_node in - let s_out_node = Column.stream c_out_node in - let s_weight = Column.stream c_weight in - let s_flags = Column.stream c_flags in - let s_id = Column.stream c_id in - let rec loop i j = - let n = Column.read s_in_node in - let t = Column.read s_out_node in - let w = Column.read s_weight in - let f = Column.read s_flags in - let way_id = Column.read s_id in - if n <> max_int then begin - let rec loop' i j = - if i < n then begin - let i = i + 1 in - if i > 0 then edge_count.{i - 1} <- j - first_edge.{i - 1}; - first_edge.{i} <- j + 1; -(*Format.eprintf "%d (%d) %d@." i n j;*) - loop' i (j + 1) - end else - j - in - let j = loop' i j in - if t <> n then begin - target.{j} <- t; - weight.{j} <- w; - flags.{j} <- f; - id.{j} <- way_id; - loop n (j + 1) - end else - loop n j - end else begin - edge_count.{i} <- j - first_edge.{i}; - i + 1 - end - in - let node_count = loop (-1) 0 in - let open_column nm = - Column.open_out (Column.named "highway/routing/unordered" nm) in - let out_source = open_column "source" in - let out_target = open_column "target" in - let out_weight = open_column "weight" in - let out_flags = open_column "flags" in - let out_shortcut = open_column "shortcut" in - { first_edge; edge_count; priority; depth; dist; remaining; - node_count; remaining_count = Column.length c_nodes; - target; weight; flags; original_edges; shortcut_node; id; - out_source; out_target; out_weight; out_flags; out_shortcut } - -let move_edge st i j = - st.target.{j} <- st.target.{i}; - st.weight.{j} <- st.weight.{i}; - st.flags.{j} <- st.flags.{i}; - st.original_edges.{j} <- st.original_edges.{i}; - st.shortcut_node.{j} <- st.shortcut_node.{i}; - st.id.{j} <- st.id.{i} - -let rec reserve st i n = -(* -Format.eprintf "reserve: %d %d@." i n; -*) - if i < st.node_count - 1 then begin - let next = first_edge st (i + 1) in - let available = next - last_edge st i - 1 in -(* -Format.eprintf "==> %d@." available; -*) -if available < 0 then Format.eprintf "%d %d %d@." i next (last_edge st i); - assert (available >= 0); - if available < n then begin - let delta = n - available in - reserve st (i + 1) delta; - for j = last_edge st (i + 1) downto first_edge st (i + 1) do - move_edge st j (j + delta) - done; - st.first_edge.{i + 1} <- st.first_edge.{i + 1} + delta - end - end else begin - assert (i + n <= Array1.dim st.target) - (*XXX*) - end - -module Heap = Binary_heap - -type dij = - { heap : Heap.t; - mutable modif : int array; - mutable pos : int } - -let init_dijkstra () = - { heap = Heap.make (); modif = Array.make 4096 0; pos = 0 } - -let push_modif state i = - let p = state.pos in - let l = Array.length state.modif in - if p = l then begin - let a = Array.make (2 * l) 0 in - Array.blit state.modif 0 a 0 l; - state.modif <- a - end; - state.modif.(p) <- i; - state.pos <- p + 1 - -let trace = ref 0 - -let dijkstra st state src ignored_node max_dist id i = - (* Restore distance *) - for i = 0 to state.pos - 1 do - st.dist.{state.modif.(i)} <- max_int - done; - state.pos <- 0; - let heap = state.heap in - Heap.insert heap src 0; -let steps = ref 0 in -let steps' = ref 0 in -if !debug then Format.eprintf ">> max %d@." max_dist; - let rec loop () = -incr steps; - let d = Heap.min_weight heap in - let n = Heap.delete_min heap in - if n <> -1 then begin -if !debug then Format.eprintf "< %d %d (%d)@." n d st.dist.{n}; - let prev_dist = st.dist.{n} in - if prev_dist > d then begin - st.dist.{n} <- d; -if prev_dist = max_int then incr steps'; - if prev_dist = max_int then push_modif state n; -if !debug then Format.eprintf "> %d children@." (st.edge_count.{n}); - for j = first_edge st n to last_edge st n do - if forward_edge st.flags.{j} then begin - let d' = d + st.weight.{j} in - if !debug then Format.eprintf "!%d %d@." d' (st.dist.{st.target.{j}}); - if d' <= max_dist then begin - let t = st.target.{j} in - if t <> ignored_node && st.dist.{t} > d' then - ( - if !debug then Format.eprintf "> %d %d [%d]@." t d' st.id.{j}; - Heap.insert heap t d' -) - end - end - done - end; - loop () - end - in - loop (); - (!steps, !steps') -(* -Format.eprintf "%d %d (%d) %d - %d -- %d@." !steps !steps' target_count (way_id id) max_dist i; -*) -(* -;(*if !trace = 2 then*) Format.eprintf "%d %d %d %d@." !steps !steps' target_count max_dist -*) -(*;print_int !steps; print_char '\n'*) - -let delete_edges st i j = - let n = ref (last_edge st i) in - let k = ref (first_edge st i) in - while !k <= !n do - if st.target.{!k} = j then begin - move_edge st !n !k; - decr n - end else - incr k - done; -(* -Format.eprintf "deleted: %d <- %d : %d %d @." i j (!n + 1 - st.first_edge.{i}) st.edge_count.{i}; -*) - st.edge_count.{i} <- !n + 1 - st.first_edge.{i} - - -let iter_remaining st f = - for i = 0 to st.remaining_count - 1 do - f i st.remaining.{i} - done - -let redundant_edge st i j dir = - let res = ref false in - let target = st.target.{j} in - for k = first_edge st i to last_edge st i do - res := - !res || - (k <> j && st.target.{k} = target && st.flags.{k} land dir <> 0 && - (st.weight.{k} < st.weight.{j} || - (st.weight.{k} = st.weight.{j} && k < j))) - done; - !res - -let prepare_node_contraction dij_state st i = - let deleted_edges = st.edge_count.{i} in - let original_deleted = ref 0 in - let added_edges = ref 0 in - let original_added_edges = ref 0 in - let shortcuts = ref [] in - for j = first_edge st i to last_edge st i do -(* Format.eprintf "%d %d@." i j;*) - original_deleted := !original_deleted + st.original_edges.{j}; - if backward_edge st.flags.{j} && not (redundant_edge st i j 1) then begin - let max_dist = ref (-1) in - for k = first_edge st i to last_edge st i do -(* -print_char '/'; -print_char ' '; -print_int i; -print_char ' '; -print_char '/'; -print_char ' '; -print_int (min st.target.{j} st.target.{k}); -print_char ' '; -print_char '/'; -print_char ' '; -print_int (max st.target.{j} st.target.{k}); -print_char '\n'; -*) - if - k <> j && forward_edge st.flags.{k} && not (redundant_edge st i k 2) - then - max_dist := max !max_dist (st.weight.{j} + st.weight.{k}) - done; - if !max_dist >= 0 then begin -let (steps, steps') = - dijkstra st dij_state (st.target.{j}) i !max_dist (st.id.{j}) i; -in -let success = ref false in - for k = first_edge st i to last_edge st i do - if - k = j || not (forward_edge st.flags.{k}) || redundant_edge st i k 2 - then () else - if st.dist.{st.target.{k}} <= st.weight.{j} + st.weight.{k} then begin -success := true -(* - Format.eprintf " %d -> %d (%d) %d@." k st.dist.{st.target.{k}} - (st.id.{k}) i; - Format.eprintf "ignored node %d <--%d(%d)-- %d --%d(%d)--> %d@." st.target.{j} (way_id st.id.{j}) st.weight.{j} i (way_id st.id.{k}) st.weight.{k} st.target.{k}; - Format.eprintf "path: %a@." (show_path st) st.target.{k} -*) - end else begin - added_edges := !added_edges + 2; - original_added_edges := - !original_added_edges + - 2 * (st.original_edges.{j} + st.original_edges.{k}); - shortcuts := - (st.target.{j}, st.target.{k}, st.weight.{j} + st.weight.{k}, - 2, st.original_edges.{j} + st.original_edges.{k}) :: !shortcuts - end -(*;Format.printf "%d / %d / %b@." steps steps' !success*) - done - end - end - done; - let p = float st.depth.{i} in - let p = - if deleted_edges > 0 then - p +. 2. *. float !added_edges /. float deleted_edges - else - p - in - let p = - if deleted_edges > 0 then - p +. 4. *. float !original_added_edges /. float !original_deleted - else - p - in -(* -Format.eprintf "%d: %f (%d)@." i p deleted_edges; -*) - (p, !shortcuts) - -let add_edge st i j weight flags original_edges shortcut_node id = -(* -Format.eprintf "add %d -> %d (%d)@." i j weight; -*) - reserve st i 1; - st.edge_count.{i} <- st.edge_count.{i} + 1; - let k = last_edge st i in - st.target.{k} <- j; - st.weight.{k} <- weight; - st.flags.{k} <- flags; - st.original_edges.{k} <- original_edges; - st.shortcut_node.{k} <- shortcut_node; - st.id.{k} <- 0 - -let contract_node dij_state st i = - let (_, shortcuts) = prepare_node_contraction dij_state st i in - for j = first_edge st i to last_edge st i do - let k = st.target.{j} in - delete_edges st k i - done; - (*XXX Could be done at a latter stage? *) - for j = first_edge st i to last_edge st i do - Column.append st.out_source i; - Column.append st.out_target st.target.{j}; - Column.append st.out_weight st.weight.{j}; - Column.append st.out_flags st.flags.{j}; - Column.append st.out_shortcut - (if st.original_edges.{j} > 1 then st.shortcut_node.{j} else -1) - done; - let shortcuts1 = Array.of_list shortcuts in -(* - Array.sort - (fun (j, k, _, _, _) (j', k', _, _, _) -> - let c = compare j j' in if c <> 0 then c else compare k k') - shortcuts1; - for l = 1 to Array.length shortcuts1 - 1 do - let (j, k, w, _, _) = shortcuts1.(l -1) in - let (j', k', w', _, _) = shortcuts1.(l) in - assert (j <> j' || k <> k') - done; -*) - let shortcuts2 = - Array.map (fun (j, k, w, _, orig_edges) -> (k, j, w, 1, orig_edges)) - shortcuts1 - in - let shortcuts = Array.append shortcuts1 shortcuts2 in - Array.sort - (fun (j, k, _, _, _) (j', k', _, _, _) -> - let c = compare j j' in if c <> 0 then c else compare k k') - shortcuts; - let len = Array.length shortcuts in -(* -Format.eprintf "%d@." len; -*) - let l = ref 0 in - while !l < len do - let (j, k, w, f, orig_edges) = shortcuts.(!l) in - if - !l + 1 < len && - let (j', k', w', _, _) = shortcuts.(!l + 1) in -(* -Format.eprintf "%d %d/%d %d/%d %d/%d@." !l j' j k' k w' w; -*) - j' = j && k' = k && w' = w - then begin - add_edge st j k w 3 orig_edges i 0; - incr l - end else - add_edge st j k w f orig_edges i 0; - incr l - done - -let update_neighbours dij_state st i = - let a = Array.make st.edge_count.{i} max_int in - let j0 = first_edge st i in - let depth = 1 + st.depth.{i} in - for j = j0 to last_edge st i do - let k = st.target.{j} in - st.depth.{k} <- max st.depth.{k} depth; - a.(j - j0) <- k - done; - Array.sort compare a; - for j = 0 to Array.length a - 1 do - let k = a.(j) in - if j = 0 || k <> a.(j - 1) then begin - st.priority.{k} <- fst (prepare_node_contraction dij_state st k) -(*;if !trace = 2 then Format.eprintf "%d / %f@." k st.priority.{k}*) - end - done - -let max_depth = ref 0 - -let update_neighbour_depths st i = - let depth = 1 + st.depth.{i} in - for j = first_edge st i to last_edge st i do - if depth > !max_depth then begin max_depth := depth; Format.eprintf "depth: %d@." depth end; - let k = st.target.{j} in - st.depth.{k} <- max st.depth.{k} depth - done - -let update_neighbour_priorities updated dij_state st i = - for j = first_edge st i to last_edge st i do - let k = st.target.{j} in - if not (Bitvect.test updated k) then begin - Bitvect.set updated k; - st.priority.{k} <- fst (prepare_node_contraction dij_state st k) -(* -;if !trace = 2 then Format.eprintf "%d / %f@." k st.priority.{k} -*) - end - done - -let check_edges st = - let rem = Bitvect.make (Array1.dim st.first_edge) in - iter_remaining st (fun _ i -> Bitvect.set rem i); - iter_remaining st - (fun _ i -> - for j = first_edge st i to last_edge st i do - let t = st.target.{j} in - if not (Bitvect.test rem t) then - Format.eprintf "DANGLING EDGE %d -> %d@." i t - else begin - let rev = ref false in - for k = first_edge st t to last_edge st t do - if st.target.{k} = i then rev := true - done; - if not !rev then - Format.eprintf "MISSING REVERSE EDGE %d -> %d@." i t - end - done) +(****) +open Osm_contraction let _ = let st = load () in diff --git a/osm/display.ml b/osm/display.ml index b9fd4cf..48bda40 100644 --- a/osm/display.ml +++ b/osm/display.ml @@ -55,1709 +55,12 @@ ===> share between ways! *) -let async_zoom = true -let async_zoom_in = true -let async_delay = 50 (*ms*) - -let debug_time = true - let _ = Printexc.record_backtrace true - let _ = Column.set_database "/tmp/osm" -let (>>) x f = f x - -(****) - -(*XXX Duplicated code...*) - -let sint_of_int i = let i' = i lsr 1 in if i land 1 = 1 then (-i' - 1) else i' - -let rec read_varint_rec a p v offs = - let i = !p in - let c = Char.code a.[i] in - incr p; - if c >= 0x80 then - read_varint_rec a p (v lor ((c land 0x7f) lsl offs)) (offs + 7) - else - v lor (c lsl offs) - -let read_varint a p = read_varint_rec a p 0 0 - -let read_signed_varint a p = sint_of_int (read_varint a p) - -let read_int_2 s pos = Char.code s.[pos] lor (Char.code s.[pos + 1] lsl 8) - -(****) - -let rec log2 x = if x <= 1 then 0 else 1 + log2 (x lsr 1) -let log2_tbl = Array.init 256 log2 -let log2_16 x = - let x' = x lsr 8 in - if x' = 0 then - Array.unsafe_get log2_tbl x - else - 8 + Array.unsafe_get log2_tbl x' -let log2_32 x = - let x' = x lsr 16 in - if x' = 0 then log2_16 x else 16 + log2_16 x' -let log2 x = - let x' = x lsr 32 in - if x' = 0 then log2_32 x else 32 + log2_32 x' - -(****) - -module Surface = Category.Make (struct - type t = - [ `Water | `Forest | `Grass | `Heath | `Rock | `Sand | `Glacier - | `Farmland | `Residential | `Commercial - | `Industrial | `Park | `Cemetery | `Parking | `Building - | `Highway_residential | `Highway_unclassified | `Highway_living_street - | `Highway_service | `Highway_pedestrian | `Highway_track - | `Highway_footway | `Highway_path ] - let list = - [ `Water; `Forest; `Grass; `Heath; `Rock; `Sand; `Glacier; - `Farmland; `Residential; `Commercial; - `Industrial; `Park; `Cemetery; `Parking; `Building; - `Highway_residential; `Highway_unclassified; `Highway_living_street; - `Highway_service; `Highway_pedestrian; `Highway_track; - `Highway_footway; `Highway_path ] -end) - -module Linear_feature = Category.Make (struct - type t = - [ `Motorway | `Trunk | `Primary | `Secondary | `Tertiary - | `Motorway_link | `Trunk_link | `Primary_link - | `Secondary_link | `Tertiary_link - | `Residential | `Unclassified | `Living_street | `Road | `Service - | `Pedestrian | `Track | `Cycleway | `Bridleway | `Footway | `Path | `Steps - | `River | `Canal | `Stream - | `Runway | `Taxiway - | `Rail | `Tram | `Subway ] - let list = - [ `Motorway; `Trunk; `Primary; `Secondary; `Tertiary; - `Motorway_link; `Trunk_link; `Primary_link; - `Secondary_link; `Tertiary_link; - `Residential; `Unclassified; `Living_street; `Road; `Service; - `Pedestrian; `Track; `Cycleway; `Bridleway; `Footway; `Path; `Steps; - `River; `Canal; `Stream; - `Runway; `Taxiway; - `Rail; `Tram; `Subway ] -end) - -let pedestrian_surface = Cairo.PNG.create "images/pedestrian.png" - -(****) - -module Bbox = Rtree.Bbox - -type bbox = Bbox.t = - { min_lat : int; max_lat : int; min_lon : int; max_lon : int } - -let bounding_box ratio x_min y_min x_max y_max = - let ratio = float ratio in - let to_lat y = truncate (Geometry.y_to_lat (y *. 10_000_000.) /. ratio) in - let to_lon x = truncate (x *. 10_000_000. /. ratio) in - { min_lat = to_lat y_min; max_lat = to_lat y_max; - min_lon = to_lon x_min; max_lon = to_lon x_max } - -(****) - -(* Old edge R-tree, now used only to find nodes for routing - (should be eventually superseded by the R-tree containing - linear features). *) - -let leaf_size = 1024 - -let decode_leaf leaves i = - let buf = String.create leaf_size in - seek_in leaves (i * leaf_size); - really_input leaves buf 0 leaf_size; - let node_len = read_int_2 buf 0 in - let edge_len = read_int_2 buf 2 in - let node_lat = Array.make (node_len / 2) 0 in - let node_lon = Array.make (node_len / 2) 0 in - let node_idx = Array.make (node_len / 2) 0 in - let edges = Array.make edge_len 0 in - let i = ref 0 in - let pos = ref 4 in - let lat = ref 0 in - let lon = ref 0 in - let idx = ref 0 in - while !pos < node_len + 4 do - let v = !lat + read_signed_varint buf pos in - node_lat.(!i) <- v; - lat := v; - let v = !lon + read_signed_varint buf pos in - node_lon.(!i) <- v; - lon := v; - let v = !idx + read_signed_varint buf pos in - node_idx.(!i) <- v; - idx := v; - incr i - done; - let node_lat = Array.sub node_lat 0 !i in - let node_lon = Array.sub node_lon 0 !i in - let node_idx = Array.sub node_idx 0 !i in - let i = ref 0 in - let node = ref 0 in - while !pos < edge_len + node_len + 4 do - let v = !node + read_signed_varint buf pos in - edges.(!i) <- v; - node := v; - incr i - done; - let edges = Array.sub edges 0 !i in - (node_lat, node_lon, node_idx, edges) - -let clamp x min max = if x < min then min else if x > max then max else x - -let to_deg x = x * 50 - -let distance_to_box bbox lat lon = - let lat' = clamp lat (to_deg bbox.min_lat) (to_deg bbox.max_lat) in - let lon' = clamp lon (to_deg bbox.min_lon) (to_deg bbox.max_lon) in - Geometry.distance lat lon lat' lon' - -let nearest_point = - let (leaves, routine_nodes) = - Rtree.open_in (Column.file_in_database "highway/r_tree") in - let leaves = open_in leaves in - Rtree.find_nearest_point routine_nodes distance_to_box - (fun j lat lon -> - let (node_lat, node_lon, node_idx, _) = decode_leaf leaves j in - let i0 = ref 0 in - let dist = ref max_int in - for i = 0 to Array.length node_lat - 1 do - let d = - Geometry.distance lat lon (to_deg node_lat.(i)) (to_deg node_lon.(i)) - in - if d < !dist then begin - dist := d; - i0 := i - end - done; - if !dist = max_int then - None - else - Some (!dist, (node_idx.(!i0), node_lat.(!i0), node_lon.(!i0)))) - -(****) - -(*XXX Temp code: rebuild ways *) - -let debug = false - -type t = (bool ref * int) list - -let initialize g = - let g' = Array.make (Array.length g) [] in - Array.iteri - (fun i l -> - List.iter - (fun j -> - let m = ref false in - g'.(i) <- (m, j) :: g'.(i); - g'.(j) <- (m, i) :: g'.(j)) - l) - g; - g' - -let odd_degree g = (ref 0, Array.map (fun l -> (List.length l) land 1 = 1) g) - -let rec next_node g i = - match g.(i) with - [] -> - -1 - | (m, j) :: r -> - g.(i) <- r; - if !m then - next_node g i - else begin - m := true; - j - end - -let rec next_odd (i, o) = - if !i = Array.length o then - -1 - else if o.(!i) then begin - o.(!i) <- false; - !i - end else begin - incr i; - next_odd (i, o) - end - -let rec find graph odd_deg i j cont lst = - let k = next_node graph j in -if debug then Format.eprintf "Going from %d to %d@." j k; - if k = -1 then begin - if i = j then begin -if debug then Format.eprintf "Loop through %d@." i; - (i :: cont, lst) - end else begin -if debug then Format.eprintf "Path %d --> %d@." i j; - (snd odd_deg).(j) <- false; - let k = next_odd odd_deg in -if debug then Format.eprintf "Continuing from %d@." k; - if k = -1 then - ([j], lst) - else begin - let (path, lst) = find graph odd_deg i k cont lst in - ([j], path :: lst) - end - end - end else begin - let (path, lst) = find graph odd_deg i k cont lst in -if debug then Format.eprintf "--@."; - find graph odd_deg j j (path) lst - end - -let rec find_circuits graph odd_deg i lst = - if i = Array.length graph then - lst - else begin - let j = next_node graph i in -if debug then Format.eprintf "Circuit through %d-%d@." i j; - if j = -1 then - find_circuits graph odd_deg (i + 1) lst - else - let (path, lst) = find graph odd_deg i j [] lst in - find_circuits graph odd_deg (i + 1) ((i :: path) :: lst) - end - -let rec find_unclosed_paths graph odd_deg lst = - let i = next_odd odd_deg in - if i <> -1 then begin - let (path, lst) = find graph odd_deg i i [] lst in - find_unclosed_paths graph odd_deg (path :: lst) - end else - lst - -let find_paths graph = - let graph = initialize graph in - let odd_deg = odd_degree graph in -if debug then Format.eprintf "Finding unclosed paths@."; - let lst = find_unclosed_paths graph odd_deg [] in -if debug then Format.eprintf "Finding circuits@."; - find_circuits graph odd_deg 0 lst - -let v = false - -let num = ref 0 -let miss = ref 0 -let find_node tbl i = - incr num; - if not tbl.(i) then begin incr miss; tbl.(i) <- true end - -let chunk tbl edges i l = - let (_, _, _, _, (node_lat, node_lon)) = edges.(i) in - let g = Array.make (Array.length node_lat) [] in - for j = i to i + l - 1 do - let (n1, n2, _, _, _) = edges.(j) in - g.(n1) <- n2 :: g.(n1) - done; - let l = find_paths g in -List.iter (fun p -> List.iter (fun i -> find_node tbl i) p) l; -if v then -List.iter - (fun p -> - Format.eprintf "+ "; - List.iter (fun n -> Format.eprintf "%d " n) p; - Format.eprintf "@.") - l; - l - -let build_paths edges = - Array.sort - (fun (_, _, cat, lay, _) (_, _, cat', lay', _) -> - let c = compare cat cat' in - if c <> 0 then c else - compare lay lay') - edges; - let len = Array.length edges in - let i0 = ref 0 in - let prev_cat = ref (-1) in - let prev_lay = ref (-1) in - let (_, _, _, _, (x, y)) = edges.(0) in - let tbl = Array.make (Array.length x) false in - let paths = ref [] in - for i = 0 to len - 1 do - let (_, _, cat, lay, _) = edges.(i) in - if (cat <> !prev_cat || lay <> !prev_lay) && i > !i0 then begin - paths := - (!prev_cat, !prev_lay, chunk tbl edges !i0 (i - !i0)) :: !paths; - i0 := i - end; - prev_cat := cat; - prev_lay := lay - done; - paths := - (!prev_cat, !prev_lay, chunk tbl edges !i0 (len - !i0)) :: !paths; -if v then Format.eprintf "----@."; - for i = 0 to Array.length tbl - 1 do - assert (tbl.(i)) - done; - Array.of_list - (List.map - (fun (cat, lay_info, l) -> - (* Extract layer and perform sign extension *) - let lay = (lay_info lsr 2) lxor 8 - 8 in - ((cat, lay, lay_info land 1 = 1, lay_info land 2 = 2), - List.map - (fun p -> - Array.of_list (List.map (fun i -> x.(i)) p), - Array.of_list (List.map (fun i -> y.(i)) p)) - l)) - !paths) - -(* R-tree containing linear features *) - -let linear_leaf_read = ref 0 - -let linear_ratio = 50 - -let decode_leaf ratio leaves i = - incr linear_leaf_read; - let leaf_size = 2048 in - let buf = String.create leaf_size in - seek_in leaves (i * leaf_size); - really_input leaves buf 0 leaf_size; - - let node_len = read_int_2 buf 0 in - let x = Array.make (node_len / 2) 0. in - let y = Array.make (node_len / 2) 0. in - let i = ref 0 in - let pos = ref 4 in - let lat = ref 0 in - let lon = ref 0 in - while !pos < node_len + 4 do - let v = !lat + read_signed_varint buf pos in - y.(!i) <- Geometry.lat_to_y (float (v * linear_ratio)); - lat := v; - let v = !lon + read_signed_varint buf pos in - x.(!i) <- float (v * linear_ratio); - lon := v; - incr i - done; - let x = Array.sub x 0 !i in - let y = Array.sub y 0 !i in - let nodes = (x, y) in - - let edge_len = read_int_2 buf 2 in - let edges = Array.make edge_len (0, 0, 0, 0, nodes) in - let i = ref 0 in - let node = ref 0 in - while !pos < edge_len + node_len + 4 do - let n1 = !node + read_signed_varint buf pos in - node := n1; - let n2 = !node + read_signed_varint buf pos in - node := n2; - let cat = Char.code buf.[!pos] in - let layer = Char.code buf.[!pos + 1] in - pos := !pos + 2; - edges.(!i) <- (n1, n2, cat, layer, nodes); - incr i - done; - Array.sub edges 0 !i - -let cache = Lru_cache.make 1000 - -let decode_leaf ratio leaves = - Lru_cache.funct cache - (fun i -> build_paths (decode_leaf linear_ratio leaves i)) - -let open_tree name = - let ratio = linear_ratio in - let (leaves, tree) = Rtree.open_in (Column.file_in_database name) in - let leaves = open_in leaves in - (linear_ratio, decode_leaf ratio leaves, tree) - -let rtrees = - [((-1., 11.5), open_tree "linear/rtrees/large_3"); - ((11.5, 12.5), open_tree "linear/rtrees/large_2"); - ((12.5, 13.5), open_tree "linear/rtrees/large_1"); - ((13.5, 30.), open_tree "linear/rtrees/all")] - -let find_linear_features level x_min y_min x_max y_max = - let lst = ref [] in - List.iter - (fun ((min_level, max_level), (ratio, decode, tree)) -> - if level > min_level && level <= max_level then begin - let bbox = bounding_box ratio x_min y_min x_max y_max in - Rtree.find tree bbox (fun i -> lst := decode i :: !lst) - end) - rtrees; - Array.concat !lst - (****) -let coastline_leaf_size = 2048 - -let decode_coastline ratio leaves i = - let buf = String.create coastline_leaf_size in - seek_in leaves (i * coastline_leaf_size); - really_input leaves buf 0 coastline_leaf_size; - let n = read_int_2 buf 0 in - let pos = ref (2 + 2 * n) in - let lat = ref 0 in - let lon = ref 0 in - let ways = ref [] in - for i = 0 to n - 1 do - let l = read_int_2 buf (2 + 2 * i) in - let x = Array.make (l + 1) 0. in - let y = Array.make (l + 1) 0. in - for j = 0 to l - 1 do - lat := !lat + read_signed_varint buf pos; - lon := !lon + read_signed_varint buf pos; -(*if j = 0 then Format.eprintf "%d %d@." !lon !lat;*) - x.(j) <- float (!lon * ratio); - y.(j) <- Geometry.lat_to_y (float (!lat * ratio)); - done; - x.(l) <- x.(0); - y.(l) <- y.(0); - ways := (x, y) :: !ways - done; - Array.of_list !ways - -let decode_coastline ratio leaves = - Lru_cache.funct cache - (fun i -> decode_coastline ratio leaves i) - -let open_tree name = - let ch = open_in (Column.file_in_database (Filename.concat name "ratio")) in - let ratio = int_of_string (input_line ch) in - close_in ch; - let (leaves, tree) = Rtree.open_in (Column.file_in_database name) in - let leaves = open_in leaves in - (ratio, decode_coastline ratio leaves, tree) - -let rtrees = - if Sys.file_exists (Column.file_in_database "coastline/rtrees/2") then - [((-1., 2.), open_tree "coastline/rtrees/2"); - ((2., 4.), open_tree "coastline/rtrees/4"); - ((4., 6.), open_tree "coastline/rtrees/6"); - ((6., 8.), open_tree "coastline/rtrees/8"); - ((8., 30.), open_tree "coastline/rtrees/small")] - else - [] - -let find_coastline level x_min y_min x_max y_max = - let lst = ref [] in - List.iter - (fun ((min_level, max_level), (ratio, decode, tree)) -> - if level > min_level && level <= max_level then begin - let bbox = bounding_box ratio x_min y_min x_max y_max in - Rtree.find tree bbox - (fun i -> lst := decode i :: !lst) - end) - rtrees; - Array.concat !lst - -(****) - -(* R-tree containing surfaces *) - -let surface_leaf_size = 2048 - -let surface_leaf_read = ref 0 - -let decode_surfaces ratio leaves i = - let buf = String.create surface_leaf_size in - seek_in leaves (i * surface_leaf_size); - really_input leaves buf 0 surface_leaf_size; - let len = read_int_2 buf 0 in - surface_leaf_read := !surface_leaf_read + len; - let buf = - if len > 1 then begin - let buf' = String.create (surface_leaf_size * len) in - String.blit buf 0 buf' 0 surface_leaf_size; - really_input leaves buf' surface_leaf_size - ((len - 1) * surface_leaf_size); - buf' - end else - buf - in - let n = read_int_2 buf 2 in - let pos = ref (4 + 4 * n) in - let lat = ref 0 in - let lon = ref 0 in - let ways = ref [] in - let category = ref 0 in - let layer = ref 0 in - let lst = ref [] in - for i = 0 to n - 1 do - let l = read_int_2 buf (4 + 4 * i) in - let cat = Char.code buf.[4 + 4 * i + 2] in - let lay = Char.code buf.[4 + 4 * i + 3] - 128 in - if cat <> 0 then begin - if !ways <> [] then lst := (!category, !layer, List.rev !ways) :: !lst; - category := cat; - layer := lay; - ways := [] - end; - let x = Array.make (l + 1) 0. in - let y = Array.make (l + 1) 0. in - for j = 0 to l - 1 do - lat := !lat + read_signed_varint buf pos; - lon := !lon + read_signed_varint buf pos; - x.(j) <- float (!lon * ratio); - y.(j) <- Geometry.lat_to_y (float (!lat * ratio)); - done; - x.(l) <- x.(0); - y.(l) <- y.(0); - ways := (x, y) :: !ways - done; - if !ways <> [] then lst := (!category, !layer, List.rev !ways) :: !lst; - !lst - -let prepare_surfaces lst = - Array.of_list - (List.map - (fun (cat, layer, ways) -> - let area = - List.fold_left - (fun a (x, y) -> a +. Geometry.polygon_area_float x y) 0. ways in - (layer, truncate (area +. 0.5), cat, ways)) - lst) - -let decode_surfaces ratio leaves = - Lru_cache.funct cache - (fun i -> prepare_surfaces (decode_surfaces ratio leaves i)) - -let open_tree name = - let ch = open_in (Column.file_in_database (Filename.concat name "ratio")) in - let ratio = int_of_string (input_line ch) in - close_in ch; - let (leaves, tree) = Rtree.open_in (Column.file_in_database name) in - let leaves = open_in leaves in - (ratio, decode_surfaces ratio leaves, tree) - -let large_surfaces = open_tree "surfaces/rtrees/large" - -let rtrees = - [((-1., 6.), open_tree "surfaces/rtrees/06"); - ((6., 7.), open_tree "surfaces/rtrees/07"); - ((7., 8.), open_tree "surfaces/rtrees/08"); - ((8., 9.), open_tree "surfaces/rtrees/09"); - ((9., 10.), open_tree "surfaces/rtrees/10"); - ((10., 12.), open_tree "surfaces/rtrees/12"); - ((12., 30.), large_surfaces); - ((15.5, 30.), open_tree "surfaces/rtrees/small")] - -let find_surfaces level x_min y_min x_max y_max = - let lst = ref [] in - List.iter - (fun ((min_level, max_level), (ratio, decode, tree)) -> - if level > min_level && level <= max_level then begin - let bbox = bounding_box ratio x_min y_min x_max y_max in - Rtree.find tree bbox - (fun i -> lst := decode i :: !lst) - end) - rtrees; - Array.concat !lst - -(**** Pixmap ***) - -type rectangle = { x : int; y : int; width : int; height: int } - -let empty_rectangle = {x = 0; y = 0; width = 0; height = 0} -let rectangle_is_empty r = r.width = 0 || r.height = 0 - -type surface = - { mutable surface : Cairo.Surface.t option; - mutable p_width : int; mutable p_height : int; - mutable valid_rect : rectangle } - -let make_surface () = - { surface = None; p_width = 0; p_height = 0; - valid_rect = empty_rectangle } - -let invalidate_surface p = p.valid_rect <- empty_rectangle - -let grow_surface pm window width height = - let width = max width pm.p_width in - let height = max height pm.p_height in - if width > pm.p_width || height > pm.p_height then begin - let old_p = pm.surface in -(* - let p = GDraw.pixmap ~width ~height ~window () in -*) - let p = - Cairo.Surface.create_similar - (Cairo.get_target (Cairo_gtk.create window#misc#window)) - Cairo.COLOR_ALPHA ~width ~height - in - let r = pm.valid_rect in - begin match old_p with - Some old_p -> - let ctx = Cairo.create p in - Cairo.set_source_surface ctx old_p 0. 0.; - Cairo.rectangle ctx 0. 0. (float r.width) (float r.height); - Cairo.set_operator ctx Cairo.SOURCE; - Cairo.fill ctx -(* - p#put_pixmap ~x:0 ~y:0 ~xsrc:0 ~ysrc:0 - ~width:r.width ~height:r.height old_p#pixmap -*) - | None -> - () - end; - pm.surface <- Some p; - pm.p_width <- width; - pm.p_height <- height - end - -let get_surface pm = match pm.surface with Some p -> p | None -> assert false - -(**** Global state ***) - -type state = - { mutable rect : rectangle; - mutable level : float; - mutable prev_rect : rectangle; - mutable prev_level : float; - mutable active : bool; - mutable timeout : Glib.Timeout.id option; - surface : surface; - mutable marker1 : (int * float * float) option; - mutable marker2 : (int * float * float) option; - mutable path : (int * int) list } - -let compute_scale st = 256. /. 360. *. 2. ** st.level - -(**** Routing ****) - -let find_marker st x y = - let scale = 256. /. 360. *. 2. ** st.level in - let x' = (float st.rect.x +. x) /. scale in - let y' = -. (float st.rect.y +. y) /. scale in - let lat = truncate (Geometry.y_to_lat (y' *. 10_000_000.)) in - let lon = truncate ((x' *. 10_000_000.)) in - Format.eprintf "%d %d@." lat lon; - let (d, (i, lat, lon)) = nearest_point lat lon in - let lat = float lat /. 200000. in - let lon = float lon /. 200000. in - Format.eprintf "%d: %f - %f %f@." i (float d /. 1000.) lat lon; - Some (i, lat, lon) - -let routing = Routing.init () -let node_lat = Column.open_in (Column.named "highway/sorted_node" "lat") -let node_lon = Column.open_in (Column.named "highway/sorted_node" "lon") - -let update_route st = - begin match st.marker1, st.marker2 with - Some (i1, _, _), Some (i2, _, _) -> - let l = Routing.find routing i1 i2 in - st.path <- - List.map (fun n -> (Column.get node_lat n, Column.get node_lon n)) l - | _ -> - () - end - -(****) - -let draw_linear_features st ctx pred stroke i = - let scale = compute_scale st in - let prev_info = ref (-1, 0, false, false) in - let count = ref 0 in - i (fun (info, ways) -> - if pred info then begin - if info <> !prev_info then - if !count > 0 then stroke !prev_info; - prev_info := info; - List.iter - (fun (x, y) -> - let len = Array.length x in - if !count > 0 && !count + len > 10000 then begin - stroke !prev_info; - count := 0 - end; - let coeff = scale /. 10_000_000. in - if st.level >= 15. && Array.length x > 2 then begin - (*XXX This could be precomputed when decoding the path *) - let ((x, y), (x1, y1), (x2, y2)) = - Line_smoothing.perform x y in - let len = Array.length x in - count := !count + 3 * len; - Cairo.move_to ctx (x.(0) *. coeff) (y.(0) *. coeff); - for k = 1 to len - 1 do - Cairo.curve_to ctx - (x1.(k - 1) *. coeff) (y1.(k - 1) *. coeff) - (x2.(k - 1) *. coeff) (y2.(k - 1) *. coeff) - (x.(k) *. coeff) (y.(k) *. coeff) - done - end else begin - count := !count + len; - Cairo.move_to ctx (x.(0) *. coeff) (y.(0) *. coeff); - for k = 1 to len - 1 do - Cairo.line_to ctx (x.(k) *. coeff) (y.(k) *. coeff) - done - end) - ways - end); - if !count > 0 then stroke !prev_info - -let draw_surfaces st ctx pred fill i = - let scale = compute_scale st in - let prev_cat = ref (-1) in - let count = ref 0 in - i - (fun ((_, area, cat, ways) as info) -> - if pred info then begin - if (cat <> !prev_cat || !count > 10000) && !prev_cat <> -1 then begin - fill !prev_cat; - count := 0; - end; - let coeff = scale /. 10_000_000. in - List.iter - (fun (x, y) -> - Cairo.move_to ctx (x.(0) *. coeff) (y.(0) *. coeff); - for i = 1 to Array.length x - 2 do - Cairo.line_to ctx (x.(i) *. coeff) (y.(i) *. coeff) - done; - Cairo.Path.close ctx; - count := !count + Array.length x) - ways; - if cat <> !prev_cat then prev_cat := cat - end); - if !prev_cat <> -1 then fill !prev_cat - -let set_surface_color ctx cat = - match cat with - `Water -> - Cairo.set_source_rgb ctx 0.52 0.94 0.94 - | `Forest -> - Cairo.set_source_rgb ctx 0.1 0.7 0.2 - | `Grass -> - Cairo.set_source_rgb ctx 0.3 0.9 0.3 - | `Heath -> - Cairo.set_source_rgb ctx 0.59 0.74 0.42 - | `Rock -> - Cairo.set_source_rgb ctx 0.37 0.42 0.49 - | `Sand -> - Cairo.set_source_rgb ctx 0.94 0.93 0.22 - | `Glacier -> - Cairo.set_source_rgb ctx 0.80 0.94 0.87 - | `Farmland -> - Cairo.set_source_rgb ctx 0.69 0.94 0.27 - | `Park -> - Cairo.set_source_rgb ctx 0.6 1.0 0.6 - | `Residential -> - Cairo.set_source_rgb ctx 0.91 0.94 0.94 - | `Commercial -> - Cairo.set_source_rgb ctx 0.94 0.78 0.78 - | `Industrial -> - Cairo.set_source_rgb ctx 0.87 0.82 0.85 - | `Parking -> - Cairo.set_source_rgb ctx 0.97 0.94 0.72 - | `Cemetery -> - Cairo.set_source_rgb ctx 0.67 0.8 0.69 - | `Building -> - Cairo.set_source_rgb ctx 0.7 0.7 0.7 - | `Highway_pedestrian | `Highway_track - | `Highway_footway | `Highway_path -> - Cairo.set_source_surface ctx pedestrian_surface 0. 0.; - Cairo.Pattern.set_extend (Cairo.get_source ctx) Cairo.Pattern.REPEAT - | `Highway_residential | `Highway_unclassified - | `Highway_living_street | `Highway_service -> - Cairo.set_source_rgb ctx 0.8 0.8 0.8 - -(****) - -let draw_coastline st ctx coastline x_min y_min x_max y_max = - let scale = compute_scale st in - let coeff = scale /. 10_000_000. in - Array.iter - (fun (x, y) -> - Cairo.move_to ctx (x.(0) *. coeff) (y.(0) *. coeff); - for i = 1 to Array.length x - 2 do - Cairo.line_to ctx (x.(i) *. coeff) (y.(i) *. coeff) - done; - Cairo.Path.close ctx) - coastline; - Cairo.set_source_rgb ctx 0.52 0.94 0.94; - Cairo.fill ctx - -(****) - -let draw_map_high_levels st ctx surfaces linear_features = - let scale = compute_scale st in - let module SP = Surface.Partition in - let partition = SP.make () in - let landuse = - SP.add_group partition - [`Forest; `Grass; `Heath; `Rock; `Sand; `Glacier; - `Farmland; `Residential; `Commercial; - `Industrial; `Park; `Cemetery; `Parking ] - in - let water = SP.add_group partition [`Water] in - let building_or_pedestrian = - SP.add_group partition - [`Building; - `Highway_pedestrian; `Highway_track; `Highway_footway; `Highway_path ] - in - let highway_areas = - SP.add_group partition - [`Highway_residential; `Highway_unclassified; - `Highway_living_street; `Highway_service ] - in -let t = Unix.gettimeofday () in - let area = Array.map (fun (_, a, _, _) -> - log2 a) surfaces in - let layer = Array.map (fun (l, _, _, _) -> l) surfaces in -let surfaces' = surfaces in - let surfaces = - SP.apply partition surfaces (fun (_, _, cat, _) -> cat) - >> SP.order_totally - >> SP.order_by area - >> SP.order_by layer - >> SP.order_by_group - >> SP.select - in -if debug_time then -Format.eprintf "Sorting surfaces (%d elements): %.3f@." (Array.length surfaces') (Unix.gettimeofday () -. t); - -let t = Unix.gettimeofday () in - let module LP = Linear_feature.Partition in - let partition = LP.make () in - let waterway = LP.add_group partition [`River; `Canal; `Stream] in - let ways = - LP.add_group partition - [ `Runway; `Taxiway; - `Rail; `Tram; `Subway; - `Pedestrian; `Track; `Cycleway; `Bridleway; `Footway; `Path; `Steps; - `Residential; `Unclassified; `Living_street; `Road; `Service; - `Tertiary_link; `Secondary_link; `Primary_link; - `Tertiary; `Secondary; `Primary; - `Trunk_link; `Motorway_link; - `Trunk; `Motorway ] - in - let layer = - Array.map - (fun ((cat, lay, is_bridge, is_tunnel), _) -> - 16 + - lay * 3 + (if is_bridge then 1 else 0) - - (if is_tunnel then 1 else 0)) - linear_features - in -let linear_features' = linear_features in - let linear_features = - LP.apply partition linear_features (fun ((cat, _, _, _), _) -> cat) - >> LP.order_totally - >> LP.order_by layer - >> LP.order_by_group - >> LP.select - in -if debug_time then -Format.eprintf "Sorting lines (%d elements): %.3f@." (Array.length linear_features') (Unix.gettimeofday () -. t); -if debug_time then begin -let n = ref 0 in -Array.iter (fun (_, _, _, ways) -> List.iter (fun (x, _) -> n := !n + Array.length x) ways) surfaces'; -Format.eprintf "Surfaces: %d nodes@." !n; - let small_area = truncate (16. *. (*64. *.*) (10_000_000. /. scale) ** 2.) in - let filter (_, area, cat, ways) = - st.level >= 15.5 || - (area > small_area && - (area > 50_000_000 || Surface.of_id cat <> `Building)) - in -let n = ref 0 in -let m = ref 0 in -Array.iter (fun ((_, _, _, ways) as info) -> if filter info then begin incr m; List.iter (fun (x, _) -> n := !n + Array.length x) ways end) surfaces'; -Format.eprintf "Surfaces: %d nodes, %d elements (%.2f%%)@." !n !m (100. *. float !m /. float (Array.length surfaces')); -let n = ref 0 in -Array.iter (fun (_, ways) -> List.iter (fun (x, _) -> n := !n + Array.length x) ways) linear_features'; -Format.eprintf "Lines: %d nodes@." !n -end; - let draw_water_lines layer = - Cairo.set_line_join ctx Cairo.JOIN_ROUND; - begin match layer with - `Tunnel -> - Cairo.set_dash ctx [|2.; 4.|] - | `Bridge -> - Cairo.set_line_cap ctx Cairo.BUTT; - Cairo.set_source_rgb ctx 0. 0. 0.; - Cairo.set_line_width ctx 4.; - linear_features >> LP.with_group waterway >> LP.iter - >> draw_linear_features st ctx - (fun (_, _, is_bridge, is_tunnel) -> is_bridge) - (fun _ -> Cairo.stroke ctx) - | `Ground -> - () - end; - Cairo.set_line_cap ctx Cairo.ROUND; - Cairo.set_source_rgb ctx 0.52 0.94 0.94; - Cairo.set_line_width ctx 2.; - linear_features >> LP.with_group waterway >> LP.iter - >> draw_linear_features st ctx - (fun (_, _, is_bridge, is_tunnel) -> - match layer with - `Tunnel -> is_tunnel - | `Bridge -> is_bridge - | `Ground -> not (is_tunnel || is_bridge)) - (fun _ -> Cairo.stroke ctx); - Cairo.set_dash ctx [||] - in - let draw_casing round i = - Cairo.set_line_cap ctx (if round then Cairo.ROUND else Cairo.BUTT); - Cairo.set_line_join ctx Cairo.JOIN_ROUND; - Cairo.set_dash ctx [||]; - let stroke (cat, _, _, is_tunnel) = - let cat = Linear_feature.of_id cat in - begin match cat with - `Pedestrian | `Track | `Cycleway | `Bridleway | `Footway - | `Path | `Steps | `Rail | `Tram | `Subway -> - Cairo.set_source_rgb ctx 1. 1. 1. - | _ -> - Cairo.set_source_rgb ctx 0. 0. 0. - end; - let line_width = - match cat with - `Trunk | `Motorway -> 11. - | `Trunk_link | `Motorway_link -> 6.5 - | `Tertiary | `Secondary | `Primary -> 8. - | `Tertiary_link | `Secondary_link | `Primary_link -> 5.5 - | `Residential | `Unclassified | `Living_street | `Road -> 6. - | `Service -> 4. - | `Pedestrian | `Track | `Cycleway - | `Bridleway | `Footway | `Path | `Steps -> 4. - | `Rail -> 5. - | `Tram | `Subway -> 6. - | _ -> assert false - in - let line_width = if round then line_width else line_width -. 0.5 in - Cairo.set_line_width ctx line_width; - Cairo.stroke ctx; - in - draw_linear_features st ctx - (fun (cat, _, is_bridge, _) -> - match Linear_feature.of_id cat with - `Motorway | `Trunk | `Primary | `Secondary | `Tertiary - | `Motorway_link | `Trunk_link | `Primary_link - | `Secondary_link | `Tertiary_link - | `Residential | `Unclassified | `Living_street | `Road - | `Service -> - true - | `Pedestrian | `Track | `Cycleway | `Bridleway | `Footway - | `Path | `Steps | `Rail | `Tram | `Subway when is_bridge -> - true - | _ -> - false) - stroke i - in - let draw_inline i = - Cairo.set_line_join ctx Cairo.JOIN_ROUND; - Cairo.set_line_cap ctx Cairo.ROUND; - let stroke (cat, _, _, is_tunnel) = - let cat = Linear_feature.of_id cat in - if cat = `Rail then begin - Cairo.set_line_cap ctx Cairo.BUTT; - Cairo.set_dash ctx [|1.; 3.|]; - Cairo.set_source_rgb ctx 0.27 0.27 0.27; - Cairo.set_line_width ctx 3.; - Cairo.stroke_preserve ctx; - Cairo.set_line_cap ctx Cairo.ROUND; - Cairo.set_dash ctx [||] - end; - begin match cat with - `Motorway | `Trunk | `Motorway_link | `Trunk_link -> - Cairo.set_source_rgb ctx 1. 0.8 0. - | `Pedestrian | `Track | `Cycleway | `Bridleway - | `Footway | `Path | `Steps -> - Cairo.set_source_rgb ctx 0. 0. 0.; - Cairo.set_line_cap ctx Cairo.BUTT; - if cat = `Steps then - Cairo.set_dash ctx [|1.; 2.|] - else - Cairo.set_dash ctx [|2.; 3.|] - | `Residential | `Unclassified | `Living_street | `Road | `Service -> - Cairo.set_source_rgb ctx 0.8 0.8 0.8 - | `Rail | `Tram -> - Cairo.set_source_rgb ctx 0.27 0.27 0.27 - | `Subway -> - Cairo.set_source_rgb ctx 0.6 0.6 0.6 - | `Runway | `Taxiway -> - Cairo.set_source_rgb ctx 0.73 0.73 0.8 - | _ -> - Cairo.set_source_rgb ctx 1. 1. 1. - end; - let line_width = - match cat with - `Motorway | `Trunk -> 6. - | `Motorway_link | `Trunk_link -> 3. - | `Primary | `Secondary | `Tertiary -> 5. - | `Primary_link | `Secondary_link | `Tertiary_link -> 2.5 - | `Residential | `Unclassified | `Living_street | `Road -> 4. - | `Service -> 2.5 - | `Rail -> 1. - | `Tram | `Subway -> 2. - | `Steps -> 3. - | `Pedestrian | `Track | `Cycleway - | `Bridleway | `Footway | `Path -> 1. - | `Runway -> max (2.5e-4 *. scale) 2. - | _ -> 2. - in - Cairo.set_line_width ctx line_width; - Cairo.stroke ctx; - begin match cat with - `Pedestrian | `Track | `Cycleway | `Bridleway - | `Footway | `Path | `Steps -> - Cairo.set_line_cap ctx Cairo.ROUND; - Cairo.set_dash ctx [||] - | _ -> - () - end - in - draw_linear_features st ctx - (fun (cat, _, _, is_tunnel) -> - not is_tunnel || Linear_feature.of_id cat <> `Subway) - stroke i - in - -let t = Unix.gettimeofday () in - (* Draw surfaces *) - let fill_surface cat = - let cat = Surface.of_id cat in - set_surface_color ctx cat; - if st.level >= 17. && cat = `Building then begin - Cairo.set_source_rgb ctx 0.71 0.71 0.71; - Cairo.fill_preserve ctx; - Cairo.set_source_rgb ctx 0.6 0.6 0.6; - Cairo.set_line_width ctx 1.; - Cairo.stroke ctx - end else - Cairo.fill ctx - in - let small_area = truncate (32. *. (10_000_000. /. scale) ** 2.) in - let filter (_, area, cat, ways) = - st.level >= 15.5 || - (area > small_area && - (area > 50_000_000 || Surface.of_id cat <> `Building)) - in - surfaces >> SP.with_group landuse >> SP.iter - >> draw_surfaces st ctx filter fill_surface; - (* Draw water lines (below buildings) *) - draw_water_lines `Ground; - surfaces >> SP.with_group water >> SP.iter - >> draw_surfaces st ctx filter fill_surface; - surfaces >> SP.with_group building_or_pedestrian >> SP.iter - >> draw_surfaces st ctx filter fill_surface; - -if debug_time then -Format.eprintf "Drawing surfaces: %.3f@." (Unix.gettimeofday () -. t); - -let t = Unix.gettimeofday () in - (* Draw linear features *) - let ways = linear_features >> LP.with_group ways in - (* Underground features *) - Cairo.Group.push ctx; - draw_water_lines `Tunnel; - ways >> LP.iter_by_key - (fun layer i -> if layer mod 3 = 0 then draw_casing true (LP.iter i)); - ways >> LP.iter_by_key - (fun layer i -> - if layer mod 3 = 0 then begin - draw_casing false (LP.iter i); draw_inline (LP.iter i) - end); - Cairo.Group.pop_to_source ctx; - Cairo.paint ~alpha:0.4 ctx; - (* Outline *) - ways >> LP.iter_by_key - (fun layer i -> - if layer mod 3 <> 0 then draw_casing true (LP.iter i)); - (* Casing/inline *) - surfaces >> SP.with_group highway_areas >> SP.iter - >> draw_surfaces st ctx filter fill_surface; - ways >> LP.iter_by_key - (fun layer i -> - if layer mod 3 <> 0 then begin - if layer mod 3 = 2 then draw_casing false (LP.iter i); - draw_inline (LP.iter i) - end); - draw_water_lines `Bridge; -if debug_time then -Format.eprintf "Drawing lines: %.3f@." (Unix.gettimeofday () -. t) - -let draw_map_medium_levels st ctx surfaces linear_features = - let scale = compute_scale st in - let module SP = Surface.Partition in - let partition = SP.make () in - let landuse = - SP.add_group partition - [`Forest; `Grass; `Heath; `Rock; `Sand; `Glacier; - `Farmland; `Residential; `Commercial; - `Industrial; `Park; `Cemetery; `Parking ] - in - let water = SP.add_group partition [`Water] in - let building_or_pedestrian = - SP.add_group partition - [`Building; - `Highway_pedestrian; `Highway_track; `Highway_footway; `Highway_path ] - in - let highway_areas = - SP.add_group partition - [`Highway_residential; `Highway_unclassified; - `Highway_living_street; `Highway_service ] - in -let t = Unix.gettimeofday () in - let area = Array.map (fun (_, a, _, _) -> - log2 a) surfaces in - let layer = Array.map (fun (l, _, _, _) -> l) surfaces in -let surfaces' = surfaces in - let surfaces = - SP.apply partition surfaces (fun (_, _, cat, _) -> cat) - >> SP.order_totally - >> SP.order_by area - >> SP.order_by layer - >> SP.order_by_group - >> SP.select - in -if debug_time then -Format.eprintf "Sorting surfaces (%d elements): %.3f@." (Array.length surfaces') (Unix.gettimeofday () -. t); - -let t = Unix.gettimeofday () in - let module LP = Linear_feature.Partition in - let partition = LP.make () in - let waterway = LP.add_group partition [`River; `Canal] in - let ignored_ways = - LP.add_group partition [ `Footway; `Steps; `Service; `Stream ] in - let minor_roads = - LP.add_group partition - [ `Runway; `Taxiway; - `Rail; `Tram; `Subway; - `Pedestrian; `Track; `Cycleway; `Bridleway; `Path; - `Residential; `Unclassified; `Living_street; `Road ] - in - let major_road_links = - LP.add_group partition - [ `Tertiary_link; `Secondary_link; `Primary_link ] in - let major_roads = - LP.add_group partition [ `Tertiary; `Secondary; `Primary ] in - let highway_links = - LP.add_group partition [ `Trunk_link; `Motorway_link ] in - let highways = LP.add_group partition [ `Trunk; `Motorway ] in -let linear_features' = linear_features in - let linear_features = - LP.apply partition linear_features (fun ((cat, _, _, _), _) -> cat) - >> LP.order_totally - >> LP.order_by_group - >> LP.select - in -if debug_time then -Format.eprintf "Sorting lines (%d elements): %.3f@." (Array.length linear_features') (Unix.gettimeofday () -. t); -if debug_time then begin -let n = ref 0 in -Array.iter (fun (_, _, _, ways) -> List.iter (fun (x, _) -> n := !n + Array.length x) ways) surfaces'; -Format.eprintf "Surfaces: %d nodes@." !n; - let small_area = truncate (32. *. (10_000_000. /. scale) ** 2.) in - let filter (_, area, cat, ways) = - (area > small_area && - (area > 50_000_000 || Surface.of_id cat <> `Building)) - in -let n = ref 0 in -let m = ref 0 in -Array.iter (fun ((_, _, _, ways) as info) -> if filter info then begin incr m; List.iter (fun (x, _) -> n := !n + Array.length x) ways end) surfaces'; -Format.eprintf "Surfaces: %d nodes, %d elements (%.2f%%)@." !n !m (100. *. float !m /. float (Array.length surfaces')); -let n = ref 0 in -Array.iter (fun (_, ways) -> List.iter (fun (x, _) -> n := !n + Array.length x) ways) linear_features'; -Format.eprintf "Lines: %d nodes@." !n -end; - let draw_water_lines () = - Cairo.set_line_cap ctx Cairo.BUTT; - Cairo.set_line_join ctx Cairo.JOIN_MITER; - Cairo.set_miter_limit ctx 1.4; - Cairo.set_source_rgb ctx 0.52 0.94 0.94; - Cairo.set_line_width ctx 2.; - linear_features >> LP.with_group waterway >> LP.iter - >> draw_linear_features st ctx - (fun (_, _, is_bridge, is_tunnel) -> not is_tunnel) - (fun _ -> Cairo.stroke ctx); - Cairo.set_dash ctx [||] - in - let draw_casing i = - Cairo.set_line_cap ctx Cairo.BUTT; - Cairo.set_line_join ctx Cairo.JOIN_MITER; - Cairo.set_miter_limit ctx 1.4; - Cairo.set_dash ctx [||]; - Cairo.set_source_rgb ctx 0. 0. 0.; - let stroke (cat, _, _, is_tunnel) = - let cat = Linear_feature.of_id cat in - let line_width = - match cat with - `Trunk | `Motorway -> 8. - | `Motorway_link | `Trunk_link - | `Primary_link | `Secondary_link | `Tertiary_link -> 3. - | _ -> 5. - in - Cairo.set_line_width ctx line_width; - Cairo.stroke ctx; - in - draw_linear_features st ctx - (fun (cat, _, is_bridge, _) -> - match Linear_feature.of_id cat with - `Motorway | `Trunk | `Primary | `Secondary | `Tertiary - | `Motorway_link | `Trunk_link - | `Primary_link | `Secondary_link | `Tertiary_link - | `Residential | `Unclassified | `Living_street | `Road -> - true - | _ -> - false) - stroke i - in - let draw_inline i = - Cairo.set_line_cap ctx Cairo.BUTT; - Cairo.set_line_join ctx Cairo.JOIN_MITER; - Cairo.set_miter_limit ctx 1.4; - let stroke (cat, _, _, is_tunnel) = - let cat = Linear_feature.of_id cat in - begin match cat with - `Motorway | `Trunk | `Motorway_link | `Trunk_link -> - Cairo.set_source_rgb ctx 1. 0.8 0. - | `Pedestrian | `Track | `Path | `Cycleway | `Bridleway -> - Cairo.set_source_rgb ctx 0. 0. 0.; - Cairo.set_line_cap ctx Cairo.BUTT; - Cairo.set_dash ctx [|2.; 3.|] - | `Residential | `Unclassified | `Living_street | `Road -> - Cairo.set_source_rgb ctx 0.7 0.7 0.7 - | `Rail | `Tram -> - Cairo.set_source_rgb ctx 0.27 0.27 0.27 - | `Subway -> - Cairo.set_source_rgb ctx 0.6 0.6 0.6 - | `Runway | `Taxiway -> - Cairo.set_source_rgb ctx 0.73 0.73 0.8 - | _ -> - Cairo.set_source_rgb ctx 1. 1. 1. - end; - let line_width = - match cat with - `Motorway | `Trunk -> 4. - | `Motorway_link | `Trunk_link - | `Primary_link | `Secondary_link | `Tertiary_link -> 1. - | `Primary | `Secondary | `Tertiary -> 3. - | `Residential | `Unclassified | `Living_street | `Road -> - 1.5 - | `Rail -> 1. - | `Tram | `Subway -> 2. - | `Pedestrian | `Track | `Cycleway | `Path | `Bridleway -> 1. - | `Runway -> max (2.5e-4 *. scale) 2. - | _ -> 1. - in - Cairo.set_line_width ctx line_width; - Cairo.stroke ctx; - begin match cat with - `Pedestrian | `Track | `Path | `Cycleway | `Bridleway -> - Cairo.set_dash ctx [||] - | _ -> - () - end - in - draw_linear_features st ctx - (fun (cat, _, _, is_tunnel) -> - not is_tunnel || Linear_feature.of_id cat <> `Subway) - stroke i - in - -let t = Unix.gettimeofday () in - (* Draw surfaces *) - let fill_surface cat = - let cat = Surface.of_id cat in - set_surface_color ctx cat; - Cairo.fill ctx - in - let small_area = truncate ((*64.*)16. *. (10_000_000. /. scale) ** 2.) in - let filter (_, area, cat, ways) = - (area > small_area && - (area > 50_000_000 || Surface.of_id cat <> `Building)) - in - surfaces >> SP.with_group landuse >> SP.iter - >> draw_surfaces st ctx filter fill_surface; - (* Draw water lines (below buildings) *) - draw_water_lines (); - surfaces >> SP.with_group water >> SP.iter - >> draw_surfaces st ctx filter fill_surface; - surfaces >> SP.with_group building_or_pedestrian >> SP.iter - >> draw_surfaces st ctx filter fill_surface; - -if debug_time then -Format.eprintf "Drawing surfaces: %.3f@." (Unix.gettimeofday () -. t); - -let t = Unix.gettimeofday () in - (* Draw linear features *) - let l = - [minor_roads; major_road_links; major_roads; highway_links; highways] in - List.iter - (fun g -> - if g <> minor_roads then - linear_features >> LP.with_group g >> LP.iter >> draw_casing; - linear_features >> LP.with_group g >> LP.iter >> draw_inline) - l; -if debug_time then -Format.eprintf "Drawing lines: %.3f@." (Unix.gettimeofday () -. t) - -let partition_high = - ([ `Footway; `Steps; `Service; - `Tram; `Subway; `Taxiway; - `Pedestrian; `Track; `Cycleway; `Bridleway; `Path; - `Residential; `Unclassified; `Living_street; `Road; - `Stream ], - [ `River; `Canal ], - [ `Runway; `Rail; `Tertiary_link; `Tertiary ], - [ `Primary_link; `Primary; `Secondary_link; `Secondary ], - [ `Trunk_link; `Motorway_link; `Trunk; `Motorway ]) - -let partition_medium = - ([ `Footway; `Steps; `Service; - `Tram; `Subway; `Taxiway; `Runway; - `Pedestrian; `Track; `Cycleway; `Bridleway; `Path; - `Residential; `Unclassified; `Living_street; `Road; - `Tertiary_link; `Tertiary; - `Canal; `Stream ], - [ `River ], - [ `Rail; `Secondary_link; `Secondary ], - [ `Primary_link; `Primary ], - [ `Trunk_link; `Motorway_link; `Trunk; `Motorway ]) - -let partition_low = - ([ `Footway; `Steps; `Service; - `Tram; `Subway; `Taxiway; `Runway; - `Pedestrian; `Track; `Cycleway; `Bridleway; `Path; - `Residential; `Unclassified; `Living_street; `Road; - `Tertiary_link; `Tertiary; `Secondary_link; `Secondary; - `Canal; `Stream ], - [ `River ], - [ `Rail ], - [ `Primary_link; `Primary ], - [ `Trunk_link; `Motorway_link; `Trunk; `Motorway ]) - -let draw_map_intermediate_levels - st ctx road_partition surfaces linear_features = - let scale = compute_scale st in - let module SP = Surface.Partition in - let partition = SP.make () in - let landuse = - SP.add_group partition - [`Forest; `Grass; `Heath; `Rock; `Sand; `Glacier; - `Farmland; `Residential; `Commercial; - `Industrial; `Park; `Cemetery; `Parking ] - in - let water = SP.add_group partition [`Water] in - let building_or_pedestrian = - SP.add_group partition - [`Building; - `Highway_pedestrian; `Highway_track; `Highway_footway; `Highway_path ] - in - let highway_areas = - SP.add_group partition - [`Highway_residential; `Highway_unclassified; - `Highway_living_street; `Highway_service ] - in -let t = Unix.gettimeofday () in - let area = Array.map (fun (_, a, _, _) -> - log2 a) surfaces in - let layer = Array.map (fun (l, _, _, _) -> l) surfaces in -let surfaces' = surfaces in - let surfaces = - SP.apply partition surfaces (fun (_, _, cat, _) -> cat) - >> SP.order_totally - >> SP.order_by area - >> SP.order_by layer - >> SP.order_by_group - >> SP.select - in -if debug_time then -Format.eprintf "Sorting surfaces (%d elements): %.3f@." (Array.length surfaces') (Unix.gettimeofday () -. t); - -let t = Unix.gettimeofday () in - let module LP = Linear_feature.Partition in - let partition = LP.make () in - let (ignored_way_cats, waterway_cats, - minor_road_cats, major_road_cats, highway_cats) = - road_partition in - let waterway = LP.add_group partition waterway_cats in - let ignored_ways = LP.add_group partition ignored_way_cats in - let minor_roads = LP.add_group partition minor_road_cats in - let major_roads = LP.add_group partition major_road_cats in - let highways = LP.add_group partition highway_cats in - -let linear_features = - let eps = (10_000_000. /. scale) in - Array.map - (fun (((cat, _, _, _) as info, ways) as item) -> - if List.memq (Linear_feature.of_id cat) ignored_way_cats then - item - else - (info, List.filter (fun (x, y) -> Array.length x > 0) (List.map (fun (x, y) -> Douglas_peucker.perform eps x y) ways))) - linear_features -in -(**) -let linear_features' = linear_features in - let linear_features = - LP.apply partition linear_features (fun ((cat, _, _, _), _) -> cat) - >> LP.order_totally - >> LP.order_by_group - >> LP.select - in -if debug_time then -Format.eprintf "Sorting lines (%d elements): %.3f@." (Array.length linear_features') (Unix.gettimeofday () -. t); -if debug_time then begin -let n = ref 0 in -Array.iter (fun (_, _, _, ways) -> List.iter (fun (x, _) -> n := !n + Array.length x) ways) surfaces'; -Format.eprintf "Surfaces: %d nodes@." !n; - let small_area = truncate (16. *. (*64. *.*) (10_000_000. /. scale) ** 2.) in - let filter (_, area, cat, ways) = - (area > small_area && - (area > 50_000_000 || Surface.of_id cat <> `Building)) - in -let n = ref 0 in -let m = ref 0 in -Array.iter (fun ((_, _, _, ways) as info) -> if filter info then begin incr m; List.iter (fun (x, _) -> n := !n + Array.length x) ways end) surfaces'; -Format.eprintf "Surfaces: %d nodes, %d elements (%.2f%%)@." !n !m (100. *. float !m /. float (Array.length surfaces')); -let n = ref 0 in -Array.iter (fun (_, ways) -> List.iter (fun (x, _) -> n := !n + Array.length x) ways) linear_features'; -Format.eprintf "Lines: %d nodes@." !n; -let n = ref 0 in -let n' = ref 0 in -let m = ref 0 in -Array.iter (fun ((cat, _, _, _), ways) -> if not (List.memq (Linear_feature.of_id cat) ignored_way_cats) then begin incr m; List.iter (fun (x, y) -> n := !n + Array.length x; let eps = (10_000_000. /. scale) in let (x', y') = Douglas_peucker.perform eps x y in n' := !n' + Array.length x') ways end) linear_features'; -Format.eprintf "Lines: %d nodes (%d), %d elements (%.2f%%)@." !n !n' !m (100. *. float !m /. float (Array.length linear_features')); -end; - - let draw_water_lines () = - Cairo.set_line_cap ctx Cairo.BUTT; - Cairo.set_line_join ctx Cairo.JOIN_MITER; - Cairo.set_miter_limit ctx 1.4; - Cairo.set_source_rgb ctx 0.52 0.94 0.94; - Cairo.set_line_width ctx 2.; - linear_features >> LP.with_group waterway >> LP.iter - >> draw_linear_features st ctx - (fun (_, _, is_bridge, is_tunnel) -> not is_tunnel) - (fun _ -> Cairo.stroke ctx); - Cairo.set_dash ctx [||] - in - let emph_secondary = List.mem `Secondary major_road_cats in - let draw_casing i = - Cairo.set_line_cap ctx Cairo.BUTT; - Cairo.set_line_join ctx Cairo.JOIN_MITER; - Cairo.set_miter_limit ctx 1.4; - Cairo.set_dash ctx [||]; - let stroke (cat, _, _, is_tunnel) = - let cat = Linear_feature.of_id cat in - Cairo.set_source_rgb ctx 0. 0. 0.; - let line_width = - match cat with - `Trunk | `Motorway -> 8. - | `Trunk_link | `Motorway_link -> 5. - | `Primary -> 5. - | `Primary_link -> 4. - | `Secondary | `Secondary_link -> 3.5 - | _ -> assert false - in - Cairo.set_line_width ctx line_width; - Cairo.stroke ctx; - in - draw_linear_features st ctx - (fun (cat, _, is_bridge, _) -> - let cat = Linear_feature.of_id cat in - match cat with - `Motorway | `Trunk | `Primary | `Secondary - | `Motorway_link | `Trunk_link | `Primary_link | `Secondary_link -> - true - | `Secondary | `Secondary_link -> - emph_secondary - | _ -> - false) - stroke i - in - let draw_inline i = - Cairo.set_line_cap ctx Cairo.BUTT; - Cairo.set_line_join ctx Cairo.JOIN_MITER; - Cairo.set_miter_limit ctx 1.4; - let stroke (cat, _, _, is_tunnel) = - let cat = Linear_feature.of_id cat in - begin match cat with - `Motorway | `Trunk | `Motorway_link | `Trunk_link -> - Cairo.set_source_rgb ctx 1. 0.8 0. - | `Tertiary | `Tertiary_link -> - Cairo.set_source_rgb ctx 0.7 0.7 0.7 - | `Secondary | `Secondary_link when not emph_secondary -> - Cairo.set_source_rgb ctx 0.7 0.7 0.7 - | `Rail -> - Cairo.set_source_rgb ctx 0.27 0.27 0.27 - | `Runway -> - Cairo.set_source_rgb ctx 0.73 0.73 0.8 - | _ -> - Cairo.set_source_rgb ctx 1. 1. 1. - end; - let line_width = - match cat with - `Motorway | `Trunk -> 4. - | `Motorway_link | `Trunk_link -> 1. - | `Primary -> 3. - | `Primary_link -> 2. - | `Rail -> 1. - | `Secondary -> 1.5 - | `Secondary_link | `Tertiary_link | `Secondary | `Tertiary -> 1.5 - | `Runway -> max (2.5e-4 *. scale) 1. - | `Taxiway -> 1. - | _ -> assert false - in - Cairo.set_line_width ctx line_width; - Cairo.stroke ctx - in - draw_linear_features st ctx - (fun (cat, _, _, is_tunnel) -> true) - stroke i - in - -let t = Unix.gettimeofday () in - (* Draw surfaces *) - let fill_surface cat = - let cat = Surface.of_id cat in - set_surface_color ctx cat; - Cairo.fill ctx - in - let small_area = truncate ((*64.*)16. *. (10_000_000. /. scale) ** 2.) in - let filter (_, area, cat, ways) = - (area > small_area && - (area > 50_000_000 || Surface.of_id cat <> `Building)) - in - surfaces >> SP.with_group landuse >> SP.iter - >> draw_surfaces st ctx filter fill_surface; - (* Draw water lines (below buildings) *) - draw_water_lines (); - surfaces >> SP.with_group water >> SP.iter - >> draw_surfaces st ctx filter fill_surface; - surfaces >> SP.with_group building_or_pedestrian >> SP.iter - >> draw_surfaces st ctx filter fill_surface; - -if debug_time then -Format.eprintf "Drawing surfaces: %.3f@." (Unix.gettimeofday () -. t); - -let t = Unix.gettimeofday () in - (* Draw linear features *) - let l = [minor_roads; major_roads; highways] in - List.iter - (fun g -> - if g <> minor_roads then - linear_features >> LP.with_group g >> LP.iter >> draw_casing; - linear_features >> LP.with_group g >> LP.iter >> draw_inline) - l; -if debug_time then -Format.eprintf "Drawing lines: %.3f@." (Unix.gettimeofday () -. t) - -(****) - -let draw_map st pm x y width height = -(* -Format.eprintf "map: %d %d %d %d@." x y width height; -*) - let ctx = Cairo.create pm in - - (* Background *) - Cairo.rectangle ctx (float x) (float y) (float width) (float height); - Cairo.clip_preserve ctx; - Cairo.set_source_rgb ctx 1. 1. 1.; - Cairo.fill ctx; - - let x = st.rect.x + x in - let y = st.rect.y + y in - - Cairo.scale ctx 1. (-1.); - Cairo.translate ctx (-. float st.rect.x) (float st.rect.y); - - let extra = 7 in - - let scale = compute_scale st in - let lon_min = float (x - extra) /. scale in - let lon_max = float (x + width + extra) /. scale in - let lat_min = -. float (y + height + extra) /. scale in - let lat_max = -. float (y - extra) /. scale in - - (* Load surfaces *) -let t = Unix.gettimeofday () in - let surfaces = find_surfaces st.level lon_min lat_min lon_max lat_max in -if debug_time then -Format.eprintf "Loading surfaces: %.3f@." (Unix.gettimeofday () -. t); - - (* Load linear features *) -let t = Unix.gettimeofday () in - let linear_features = - find_linear_features st.level lon_min lat_min lon_max lat_max in -if debug_time then -Format.eprintf "Loading lines: %.3f@." (Unix.gettimeofday () -. t); - - (* Load coastline *) -let t = Unix.gettimeofday () in - let coastline = - find_coastline st.level lon_min lat_min lon_max lat_max in -if debug_time then -Format.eprintf "Loading coastline: %.3f@." (Unix.gettimeofday () -. t); - -if debug_time then -Format.eprintf "Surfaces: %d / lines: %d@." (2 * !surface_leaf_read) (2 * !linear_leaf_read); - - draw_coastline st ctx coastline lon_min lat_min lon_max lat_max; - if st.level >= 14.5 then - draw_map_high_levels st ctx surfaces linear_features - else if st.level >= 13.5 then - draw_map_medium_levels st ctx surfaces linear_features - else if st.level >= 12.5 then - draw_map_intermediate_levels - st ctx partition_high surfaces linear_features - else if st.level >= 11.5 then - draw_map_intermediate_levels - st ctx partition_medium surfaces linear_features - else - draw_map_intermediate_levels - st ctx partition_low surfaces linear_features - -(****) - -let draw_route st ctx = - let scale = compute_scale st in - Cairo.scale ctx 1. (-1.); - Cairo.translate ctx (-. float st.rect.x) (float st.rect.y); - let path = - List.map - (fun (lat, lon) -> - let approx x = - float ((x + linear_ratio / 2 - 1) / linear_ratio * linear_ratio) - in - (approx lon /. 10_000_000., - Geometry.lat_to_y (approx lat) /. 10_000_000.)) - st.path - in - begin match path with - [] -> - () - | (x, y) :: rem -> - Cairo.move_to ctx (x *. scale) (y *. scale); - List.iter - (fun (x, y) -> Cairo.line_to ctx (x *. scale) (y *. scale)) - rem; - Cairo.set_line_width ctx 2.; - Cairo.set_source_rgb ctx 0. 0. 1.; - Cairo.stroke ctx - end; - - begin match st.marker1 with - Some (_, lat, lon) -> - let x = lon in - let y = Geometry.lat_to_y (lat *. 10_000_000.) /. 10_000_000. in - Cairo.arc ctx (x *. scale) (y *. scale) 4. 0. (2. *. Geometry.pi); - Cairo.set_source_rgb ctx 1. 0. 0.; - Cairo.fill ctx - | None -> - () - end; - begin match st.marker2 with - Some (_, lat, lon) -> - let x = lon in - let y = Geometry.lat_to_y (lat *. 10_000_000.) /. 10_000_000. in - Cairo.arc ctx (x *. scale) (y *. scale) 4. 0. (2. *. Geometry.pi); - Cairo.set_source_rgb ctx 0. 0.6 0.; - Cairo.fill ctx - | None -> - () - end +open Osm_display let _ = let width = 512 in @@ -1773,7 +76,7 @@ let st = let lat = ref 48.850 in let lon = ref 2.350 in begin - let (ratio, _, tree) = large_surfaces in + let (ratio, _, tree) = Lazy.force large_surfaces in let bbox = Rtree.bounding_box tree in let c x = truncate (x *. 10_000_000. /. float ratio +. 0.5) in Format.eprintf "%a %d %d@." Bbox.print bbox (c !lat) (c !lon); @@ -1787,7 +90,7 @@ let scale = compute_scale st in st.rect <- { st.rect with x = truncate (!lon *. scale) - width / 2; - y = - truncate (Geometry.lat_to_y (!lat *. 10_000_000.) /. 10_000_000. *. scale) - height / 2 }; + y = - truncate (Osm_geometry.lat_to_y (!lat *. 10_000_000.) /. 10_000_000. *. scale) - height / 2 }; Format.eprintf "%d %d@." st.rect.x st.rect.y; ignore (GMain.Main.init ()); @@ -1923,7 +226,7 @@ let t = Unix.gettimeofday () in Cairo.save ctx; (* Workaround for a Cairo bug (in ATI Catalyst drivers?): *) if st.active then Cairo.set_operator ctx Cairo.SOURCE; - Cairo.fill_preserve ctx; + Cairo.fill_preserve ctx; Cairo.restore ctx; Cairo.clip ctx; draw_route st ctx; diff --git a/osm/highway.ml b/osm/highway.ml index 619de19..69bd87c 100644 --- a/osm/highway.ml +++ b/osm/highway.ml @@ -31,7 +31,7 @@ let compute_order out latitude longitude = if lat <> max_int then begin let lat = lat + 90_0000000 in let lon = lon + 180_0000000 in - Column.append o (Geometry.hilbert_coordinate lat lon); + Column.append o (Osm_geometry.hilbert_coordinate lat lon); loop () end in @@ -211,7 +211,7 @@ List.iter (fun (k,v) -> Format.eprintf "%s(%d)=%s(%d)@." (Dictionary.get dict k) if w = w' && info.P.speed > 0. then begin if info.P.backward_speed <= 0. then info.P.backward_speed <- info.P.speed; - let dist = Geometry.distance lat lon lat' lon' in + let dist = Osm_geometry.distance lat lon lat' lon' in (* Format.eprintf "%d %d %d %d %d@." lon lat lon' lat' dist; *) diff --git a/osm/category.ml b/osm/lib/osm_category.ml similarity index 100% rename from osm/category.ml rename to osm/lib/osm_category.ml diff --git a/osm/category.mli b/osm/lib/osm_category.mli similarity index 100% rename from osm/category.mli rename to osm/lib/osm_category.mli diff --git a/osm/clipping.ml b/osm/lib/osm_clipping.ml similarity index 100% rename from osm/clipping.ml rename to osm/lib/osm_clipping.ml diff --git a/osm/clipping.mli b/osm/lib/osm_clipping.mli similarity index 100% rename from osm/clipping.mli rename to osm/lib/osm_clipping.mli diff --git a/osm/lib/osm_coastline.ml b/osm/lib/osm_coastline.ml new file mode 100644 index 0000000..d3a4e03 --- /dev/null +++ b/osm/lib/osm_coastline.ml @@ -0,0 +1,504 @@ +(* OSM tools + * Copyright (C) 2013 Jérôme Vouillon + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, with linking exception; + * either version 2.1 of the License, or (at your option) any later version. + * + * 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 Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + *) + +(* +TODO +==== +- simplification should depend on latitude + (should be done after projection) +*) + +let leaf_size = 2048 +let max_polygon_size = 200 + +(****) + +let read_big_int ch = + let i3 = input_byte ch in + let i2 = input_byte ch in + let i1 = input_byte ch in + let i0 = input_byte ch in + i0 lor (i1 lsl 8) lor (i2 lsl 16) lor (i3 lsl 24) + +let read_lit_int ch = + let i0 = input_byte ch in + let i1 = input_byte ch in + let i2 = input_byte ch in + let i3 = input_byte ch in + i0 lor (i1 lsl 8) lor (i2 lsl 16) lor (i3 lsl 24) + +let read_float ch = + let i0 = Int64.of_int (read_lit_int ch) in + let i1 = Int64.of_int (read_lit_int ch) in + Int64.float_of_bits (Int64.logor i0 (Int64.shift_left i1 32)) + +let build_rings l = + let (lst, rem) = + List.partition + (fun (x, y) -> + let l = Array.length x in x.(0) = x.(l - 1) && y.(0) = y.(l - 1)) + l + in + let succ = Hashtbl.create 17 in + let pred = Hashtbl.create 17 in + let rem = Array.of_list rem in + Array.iteri + (fun i (x, y) -> + let l = Array.length x in + Hashtbl.add succ (x.(0), y.(0)) i; + Hashtbl.add pred (x.(l - 1), y.(l - 1)) i) + rem; + let lst = ref lst in + let rec follow i l = + let (x, y) = rem.(i) in + let len = Array.length x in + try + let k = Hashtbl.find pred (x.(0), y.(0)) in + let l = (Array.sub x 1 (len - 1), Array.sub y 1 (len - 1)) :: l in + follow k l + with Not_found -> + (x, y) :: l + in + for i = 0 to Array.length rem - 1 do + let (x, y) = rem.(i) in + let len = Array.length x in + if not (Hashtbl.mem succ (x.(len - 1), y.(len - 1))) then begin + let l = follow i [] in + let x = Array.concat (List.map fst l) in + let y = Array.concat (List.map snd l) in + let len = Array.length x in + let (x, y) = + if x.(0) < -1799999000. && x.(len - 1) < -1799999000. then begin + x.(0) <- -1800000000.; x.(len - 1) <- -1800000000.; + (Array.append x [|x.(0)|], Array.append y [|y.(0)|]) + end else if + x.(0) > 1799999000. && x.(len - 1) > 1799999000. + then begin + x.(0) <- 1800000000.; x.(len - 1) <- 1800000000.; + (Array.append x [|x.(0)|], Array.append y [|y.(0)|]) + end else if + x.(0) = -1800000000. && x.(len -1) = 1800000000. + then begin + (* Antartica *) + (Array.append x [|1800000000.; -1800000000.; x.(0)|], + Array.append y [| 850000000.; 850000000.; y.(0)|]) + end else begin +Format.eprintf "%f %f - %f %f@." x.(0) y.(0) x.(len - 1) y.(len - 1); + assert false + end + in + lst := (x, y) :: !lst + end + done; + !lst + +let read_record ch = + let num = read_big_int ch in + let len = read_big_int ch in + ignore (num, len); + let typ = read_lit_int ch in + assert (typ = 3); (* polyline *) + + for i = 0 to 3 do ignore (read_float ch) done; (* bbox *) + + let num_parts = read_lit_int ch in + assert (num_parts = 1); + + let num_points = read_lit_int ch in + + ignore (read_lit_int ch); (* first and unique part *) + + let lon = Array.make (num_points - 1) 0. in + let lat = Array.make (num_points - 1) 0. in + + for i = 0 to num_points - 2 do + lon.(i) <- 10_000_000. *. read_float ch; + lat.(i) <- 10_000_000. *. read_float ch + done; + let lon0 = 10_000_000. *. read_float ch in + let lat0 = 10_000_000. *. read_float ch in + ((lon, lat), (lon0, lat0)) + +let open_file () = + let ch = open_in Sys.argv.(1) in + let magic = read_big_int ch in + if magic = 9994 then begin + seek_in ch 0; + ch + end else if magic = 0x504b0304 then begin + close_in ch; + let z = Zip.open_in Sys.argv.(1) in + let e = + try + Zip.find_entry z "coastlines-split-4326/lines.shp" + with Not_found -> + Util.fail "Zip entry 'coastlines-split-4326/lines.shp' not found" + in + let (r, w) = Unix.pipe () in + match Unix.fork () with + 0 -> + Unix.close r; + let ch = Unix.out_channel_of_descr w in + Zip.copy_entry_to_channel z e ch; + flush ch; + exit 0 + | _ -> + Unix.close w; + Unix.in_channel_of_descr r + end else + Util.fail (Format.sprintf "bad file format") + +let load () = + Format.eprintf "==== Loading.@."; + + let ch = open_file () in + + for i = 0 to 24 do + ignore (read_big_int ch) + done; + +(* + let a = read_float ch in + Format.eprintf "%f@." a; + let a = read_float ch in + Format.eprintf "%f@." a; + let a = read_float ch in + Format.eprintf "%f@." a; + let a = read_float ch in + Format.eprintf "%f@." a; +*) + + let lst = ref [] in + let cur = ref [] in + let lon0' = ref 0. in + let lat0' = ref 0. in + begin try + while true do + let ((lon, lat), (lon0, lat0)) = read_record ch in + if !cur <> [] && (lat.(0) <> !lat0' || lon.(0) <> !lon0') then begin + let path = List.rev (([|!lon0'|], [|!lat0'|]) :: !cur) in + let lon'' = Array.concat (List.map fst path) in + let lat'' = Array.concat (List.map snd path) in + lst := (lon'', lat'') :: !lst; + cur := [] + end; + cur := (lon, lat) :: !cur; + lon0' := lon0; + lat0' := lat0 + done + with End_of_file -> () end; + let polys = build_rings !lst in +(* +Format.eprintf "%d@." (List.length polys); +*) + polys + +(**** Quickselect ****) + +let swap (tbl : float array) i j = + let v = tbl.(i) in tbl.(i) <- tbl.(j); tbl.(j) <- v + +let partition (tbl : float array) left right pivot = + let p = tbl.(pivot) in + swap tbl pivot right; + let j = ref left in + for i = left to right - 1 do + if tbl.(i) < p then begin + swap tbl i !j; + incr j + end + done; + swap tbl right !j; + !j + +let rec select tbl left right k = + if left = right then + tbl.(left) + else begin + let pivot = Random.int (right - left + 1) + left in + let pivot = partition tbl left right pivot in + if pivot = k then + tbl.(pivot) + else if k < pivot then + select tbl left (pivot - 1) k + else + select tbl (pivot + 1) right k + end + +let median tbl = + let len = Array.length tbl in + select tbl 0 (len - 1) (len / 2) + +let poly_median polys = + let l = List.rev_map fst polys in (* Tail recursive! *) + median (Array.concat l) + +let stats polys = + let a = Array.concat (List.rev_map fst polys) in (* Tail recursive! *) + let min = ref max_float in + let max = ref (-. max_float) in + for i = 0 to Array.length a - 1 do + if a.(i) < !min then min := a.(i); + if a.(i) > !max then max := a.(i) + done; + let min' = ref max_float in + let max' = ref (-. max_float) in + for i = 0 to Array.length a - 1 do + if a.(i) < !min' && a.(i) > !min then min' := a.(i); + if a.(i) > !max' && a.(i) < !max then max' := a.(i) + done; + (!min', median a, !max') + +(****) + +let swap_coord polys = List.map (fun (x, y) -> (y, x)) polys + +let ring_bbox (x, y) : float * float * float * float = + let xmin = ref x.(0) in + let ymin = ref y.(0) in + let xmax = ref x.(0) in + let ymax = ref y.(0) in + for i = 1 to Array.length x - 1 do + if x.(i) < !xmin then xmin := x.(i); + if y.(i) < !ymin then ymin := y.(i); + if x.(i) > !xmax then xmax := x.(i); + if y.(i) > !ymax then ymax := y.(i) + done; + (!xmin, !ymin, !xmax, !ymax) + +let min' x y : float = if x < y then x else y +let max' x y : float = if x < y then y else x + +let bbox polys = + match polys with + [] -> assert false + | r :: rem -> + List.fold_left + (fun (xmin, ymin, xmax, ymax) r -> + let (xmin', ymin', xmax', ymax') = ring_bbox r in + (min' xmin xmin', min' ymin ymin', + max' xmax xmax', max' ymax ymax')) + (ring_bbox r) rem + +let size polys = List.fold_left (fun n (x, _) -> n + Array.length x) 0 polys +let size polys = + List.fold_left (fun n (x, _) -> max n (Array.length x)) 0 polys + +(****) + +let d = false +let next_perc = ref 0.005 + +let rec split_rec dir t0 t1 t2 polys l = + let sz = size polys in +if d then Format.eprintf "size: %d (%d)@." sz (List.length polys); + if sz > max_polygon_size then begin + let (xmin, ymin, xmax, ymax) = bbox polys in +if d then Format.eprintf "%f %f - %f %f@." xmin ymin xmax ymax; + if xmax -. xmin < ymax -. ymin then + split_rec (not dir) t0 t1 t2 (swap_coord polys) l + else begin + let m' = ((xmax +. xmin) /. 2.) in + let (m1, m, m2) = stats polys in +if d then Format.eprintf "%f %f %f / %f@." m1 m m2 m'; + let m = + if m2 < m' then m2 else + if m1 > m' then m1 else + if m > (xmin +. m') /. 2. && m < (xmax +. m') /. 2. then + m + else + m' + in +if d then Format.eprintf "Cutting at %f@." m; + let (left, right) = Osm_clipping.perform m polys in + let l = split_rec dir t0 t1 ((t1 +. t2) /. 2.) left l in + split_rec dir t0 ((t1 +. t2) /. 2.) t2 right l + end + end else begin +if t2 >= !next_perc then begin +Util.set_msg + (Format.sprintf "splitting polygons: %s %.0f%% eta %.0fs" + (Util.progress_bar t2) (t2 *. 100.) + ((1. -. t2) *. (Unix.gettimeofday () -. t0) /. t2)); +(*Format.eprintf "%.1f%%@." t2;*) +next_perc := !next_perc +. 0.005 +end; + let polys = if dir then polys else swap_coord polys in +(* +if show then begin +List.iter + (fun (x, y) -> + Cairo.move_to ctx (x.(0) /. 10_000_000.) (y.(0) /. 10_000_000.); + for i = 1 to Array.length x - 2 do + Cairo.line_to ctx (x.(i) /. 10_000_000.) (y.(i) /. 10_000_000.) + done; + Cairo.Path.close ctx) + polys; +Cairo.fill ctx; +if t2 > 0.2 then +(Cairo.Surface.finish surface; exit 0) +end; +*) + polys @ l + end + +(****) + +module Bbox = Rtree.Bbox + +let int_of_sint i = if i >= 0 then 2 * i else - 2 * i - 1 + +let rec write_varint a p v = + if v < 128 then begin + a.[p] <- Char.chr (v land 127); + p + 1 + end else begin + a.[p] <- Char.chr ((v land 127) + 128); + write_varint a (p + 1) (v lsr 7) + end + +let write_signed_varint a p v = write_varint a p (int_of_sint v) + +let output_int_2 ch v = + output_byte ch (v land 0xff); + output_byte ch (v lsr 8) + +let open_rtree name ratio = + let (nm, st) = Rtree.open_out name in + let ch = open_out (Column.file_in_database (Filename.concat name "ratio")) in + Printf.fprintf ch "%d\n" ratio; + close_out ch; + let ch = open_out nm in + let lengths = Array.make (leaf_size / 4) 0 in + let n = ref 0 in + let buf = String.make (2 * leaf_size) '\000' in + let pos = ref 0 in + let bbox = ref Bbox.empty in + let last_lat = ref 0 in + let last_lon = ref 0 in + let flush_ways () = + output_int_2 ch !n; + for i = 0 to !n - 1 do + output_int_2 ch lengths.(i); + done; + Rtree.append st !bbox; + output ch buf 0 (leaf_size - !n * 2 - 2); + n := 0; + pos := 0; + bbox := Bbox.empty; + last_lat := 0; + last_lon := 0 + in + let rec add_polygon (lon, lat) = + let pos' = ref !pos in + for i = 0 to Array.length lat - 2 do + pos' := write_signed_varint buf !pos' (lat.(i) - !last_lat); + last_lat := lat.(i); + pos' := write_signed_varint buf !pos' (lon.(i) - !last_lon); + last_lon := lon.(i) + done; + let pos' = !pos' in + if (!n + 1) * 2 + 2 + pos' > leaf_size then begin + assert (!n > 0); + flush_ways (); + add_polygon (lon, lat) + end else begin + lengths.(!n) <- Array.length lat - 1; + incr n; + pos := pos'; + let max x y : int = if x > y then x else y in + let min x y : int = if x > y then y else x in + let bbox' = + { Bbox. + min_lat = Array.fold_left min max_int lat; + max_lat = Array.fold_left max min_int lat; + min_lon = Array.fold_left min max_int lon; + max_lon = Array.fold_left max min_int lon } + in + bbox := Bbox.union !bbox bbox'; + if !n * 2 + 2 + pos' > leaf_size then flush_ways () + end + in + let close () = + if !n > 0 then flush_ways (); + Rtree.close_out st; + close_out ch + in + (add_polygon, close) + +(****) + +let simplify level polys = + Format.eprintf "==== Simplifying.@."; + (* We multiply by 4. to get good approximations near the poles *) + let scale = 256. /. 360. *. 2. ** level *. 4. in + let small_area = (10_000_000. /. scale) ** 2. in + let ratio = float (truncate (10_000_000. /. scale /. 2.)) in + (truncate ratio, + List.fold_left + (fun lst (lon, lat) -> + let (lon, lat) = Osm_douglas_peucker.perform ratio lon lat in + if + abs_float (Osm_geometry.polygon_area_float lon lat) <= small_area + then + lst + else + (lon, lat) :: lst) + [] polys) + +(****) + +let rescale_ring ratio (lon, lat) = + let l = Array.length lat in + let delta = ratio / 2 - 1 in + let lat' = Array.make l 0 in + let lon' = Array.make l 0 in + for i = 0 to l - 1 do + lat'.(i) <- (truncate lat.(i) + delta) / ratio; + lon'.(i) <- (truncate lon.(i) + delta) / ratio + done; + (lon', lat') + +let build_rtree name ratio polys = + Format.eprintf "==== Splitting.@."; + next_perc := 0.005; + let l = split_rec true (Unix.gettimeofday ()) 0. 1. polys [] in + Util.set_msg ""; +(* +Format.eprintf "==> %d@." (List.length l); +*) + Format.eprintf "==== Sorting.@."; + let a = Array.of_list l in + let a = + Array.map + (fun p -> + let (lon_min, lat_min, lon_max, lat_max) = ring_bbox p in + let lat = truncate ((lat_min +. lat_max) /. 2.) + 90_0000000 in + let lon = truncate ((lon_min +. lon_max) /. 2.) + 180_0000000 in + (Osm_geometry.hilbert_coordinate lat lon, p)) + a + in + Array.sort (fun (x, _) (x', _) -> compare x x') a; + + Format.eprintf "==== Writing B-tree.@."; + let (add_polygon, close) = + open_rtree ("coastline/rtrees/" ^ name) (*"small"*) ratio in + Array.iter + (fun (_, (lon, lat)) -> add_polygon (rescale_ring ratio (lon, lat))) + a; + close () diff --git a/osm/lib/osm_contraction.ml b/osm/lib/osm_contraction.ml new file mode 100644 index 0000000..f75e7e4 --- /dev/null +++ b/osm/lib/osm_contraction.ml @@ -0,0 +1,550 @@ +(* OSM tools + * Copyright (C) 2013 Jérôme Vouillon + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, with linking exception; + * either version 2.1 of the License, or (at your option) any later version. + * + * 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 Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + *) + +(* +TODO: +- use memory mappings +- implement search +- oneway edges + + +- Dynamic graph: add/remove edges + +Dynamic graph +- bucket based +- compress buckets when mostly filled +- + + +Operations: +- reserve k edges +- remove an edge +- add an edge +- (remove all edges) + +Keep track of the usage of each page: +--> too small -> compact/merge +--> too large -> split + + +Page full --> balance with next page +Next page also full --> split page + + +Mapping page |-> set of node on this page + +Compress: +- copy to a temp location then blit back +- update nodes + +*) + +let debug = ref false + +module Array1 = Bigarray.Array1 + +let way_ids = lazy (Column.open_in (Column.named "highway" "way/id")) +let way_id k = Column.get (Lazy.force way_ids) k + +type t = + { (* Nodes *) + first_edge : (int, Bigarray.int_elt, Bigarray.c_layout) Array1.t; + edge_count : (int, Bigarray.int_elt, Bigarray.c_layout) Array1.t; + priority : (float, Bigarray.float64_elt, Bigarray.c_layout) Array1.t; + depth : (int, Bigarray.int_elt, Bigarray.c_layout) Array1.t; + dist : (int, Bigarray.int_elt, Bigarray.c_layout) Array1.t; + remaining : (int, Bigarray.int_elt, Bigarray.c_layout) Array1.t; + node_count : int; + mutable remaining_count : int; + (* Edges *) + target : (int, Bigarray.int_elt, Bigarray.c_layout) Array1.t; + weight : (int, Bigarray.int_elt, Bigarray.c_layout) Array1.t; + flags : (int, Bigarray.int8_unsigned_elt, Bigarray.c_layout) Array1.t; + original_edges : (int, Bigarray.int_elt, Bigarray.c_layout) Array1.t; + shortcut_node : (int, Bigarray.int_elt, Bigarray.c_layout) Array1.t; + id : (int, Bigarray.int_elt, Bigarray.c_layout) Array1.t; + (* Output *) + out_source : Column.output_stream; + out_target : Column.output_stream; + out_weight : Column.output_stream; + out_flags : Column.output_stream; + out_shortcut : Column.output_stream } + +let first_edge st i = st.first_edge.{i} +let last_edge st i = st.first_edge.{i} + st.edge_count.{i} - 1 + +let forward_edge f = f land 2 <> 0 +let backward_edge f = f land 1 <> 0 +let reverse_edge f = ((f lsl 1) lor (f lsr 1)) land 3 + +let show_path st f i = + let rec follow i = + let d = st.dist.{i} in + if d > 0 then begin + let j = ref (-1) in + for k = first_edge st i to last_edge st i do + if + forward_edge st.flags.{k} && + st.dist.{st.target.{k}} + st.weight.{k} = d + then + j := k + done; +if (!j = -1) then begin + for k = first_edge st i to last_edge st i do + Format.eprintf "XXX %d: %d %d %d@." st.target.{k} st.dist.{st.target.{k}} st.weight.{k} d + done; +end; + Format.fprintf f "--%d(%d)--> %d" (way_id st.id.{!j}) st.weight.{!j} st.target.{!j}; + follow st.target.{!j} + end + in + Format.fprintf f "%d " i; + follow i + +let int_array len = Array1.create Bigarray.int Bigarray.c_layout len +let int8_array len = Array1.create Bigarray.int8_unsigned Bigarray.c_layout len + +let load () = + let c_nodes = Column.open_in (Column.named "highway" "node/idx") in + let c_in_node = Column.open_in (Column.named "highway" "edge/1bis") in + let c_out_node = Column.open_in (Column.named "highway" "edge/2bis") in + let c_weight = Column.open_in (Column.named "highway" "edge/weight_bis") in + let c_flags = Column.open_in (Column.named "highway" "edge/flags_bis") in + let c_id = Column.open_in (Column.named "highway" "edge/idx_bis") in + + let first_edge = int_array (Column.length c_nodes) in + let edge_count = int_array (Column.length c_nodes) in + let priority = + Array1.create Bigarray.float64 Bigarray.c_layout (Column.length c_nodes) in + let depth = int_array (Column.length c_nodes) in + let dist = int_array (Column.length c_nodes) in + let remaining = int_array (Column.length c_nodes) in + Array1.fill depth 0; + Array1.fill dist max_int; + for i = 0 to Array1.dim remaining - 1 do + remaining.{i} <- i + done; + + let len = Column.length c_in_node + Column.length c_nodes + 10 in + let target = int_array len in + let weight = int_array len in + let flags = int8_array len in + let original_edges = int_array len in + let shortcut_node = int_array len in + let id = int_array len in + Array1.fill original_edges 1; + + let s_in_node = Column.stream c_in_node in + let s_out_node = Column.stream c_out_node in + let s_weight = Column.stream c_weight in + let s_flags = Column.stream c_flags in + let s_id = Column.stream c_id in + let rec loop i j = + let n = Column.read s_in_node in + let t = Column.read s_out_node in + let w = Column.read s_weight in + let f = Column.read s_flags in + let way_id = Column.read s_id in + if n <> max_int then begin + let rec loop' i j = + if i < n then begin + let i = i + 1 in + if i > 0 then edge_count.{i - 1} <- j - first_edge.{i - 1}; + first_edge.{i} <- j + 1; +(*Format.eprintf "%d (%d) %d@." i n j;*) + loop' i (j + 1) + end else + j + in + let j = loop' i j in + if t <> n then begin + target.{j} <- t; + weight.{j} <- w; + flags.{j} <- f; + id.{j} <- way_id; + loop n (j + 1) + end else + loop n j + end else begin + edge_count.{i} <- j - first_edge.{i}; + i + 1 + end + in + let node_count = loop (-1) 0 in + let open_column nm = + Column.open_out (Column.named "highway/routing/unordered" nm) in + let out_source = open_column "source" in + let out_target = open_column "target" in + let out_weight = open_column "weight" in + let out_flags = open_column "flags" in + let out_shortcut = open_column "shortcut" in + { first_edge; edge_count; priority; depth; dist; remaining; + node_count; remaining_count = Column.length c_nodes; + target; weight; flags; original_edges; shortcut_node; id; + out_source; out_target; out_weight; out_flags; out_shortcut } + +let move_edge st i j = + st.target.{j} <- st.target.{i}; + st.weight.{j} <- st.weight.{i}; + st.flags.{j} <- st.flags.{i}; + st.original_edges.{j} <- st.original_edges.{i}; + st.shortcut_node.{j} <- st.shortcut_node.{i}; + st.id.{j} <- st.id.{i} + +let rec reserve st i n = +(* +Format.eprintf "reserve: %d %d@." i n; +*) + if i < st.node_count - 1 then begin + let next = first_edge st (i + 1) in + let available = next - last_edge st i - 1 in +(* +Format.eprintf "==> %d@." available; +*) +if available < 0 then Format.eprintf "%d %d %d@." i next (last_edge st i); + assert (available >= 0); + if available < n then begin + let delta = n - available in + reserve st (i + 1) delta; + for j = last_edge st (i + 1) downto first_edge st (i + 1) do + move_edge st j (j + delta) + done; + st.first_edge.{i + 1} <- st.first_edge.{i + 1} + delta + end + end else begin + assert (i + n <= Array1.dim st.target) + (*XXX*) + end + +module Heap = Binary_heap + +type dij = + { heap : Heap.t; + mutable modif : int array; + mutable pos : int } + +let init_dijkstra () = + { heap = Heap.make (); modif = Array.make 4096 0; pos = 0 } + +let push_modif state i = + let p = state.pos in + let l = Array.length state.modif in + if p = l then begin + let a = Array.make (2 * l) 0 in + Array.blit state.modif 0 a 0 l; + state.modif <- a + end; + state.modif.(p) <- i; + state.pos <- p + 1 + +let trace = ref 0 + +let dijkstra st state src ignored_node max_dist id i = + (* Restore distance *) + for i = 0 to state.pos - 1 do + st.dist.{state.modif.(i)} <- max_int + done; + state.pos <- 0; + let heap = state.heap in + Heap.insert heap src 0; +let steps = ref 0 in +let steps' = ref 0 in +if !debug then Format.eprintf ">> max %d@." max_dist; + let rec loop () = +incr steps; + let d = Heap.min_weight heap in + let n = Heap.delete_min heap in + if n <> -1 then begin +if !debug then Format.eprintf "< %d %d (%d)@." n d st.dist.{n}; + let prev_dist = st.dist.{n} in + if prev_dist > d then begin + st.dist.{n} <- d; +if prev_dist = max_int then incr steps'; + if prev_dist = max_int then push_modif state n; +if !debug then Format.eprintf "> %d children@." (st.edge_count.{n}); + for j = first_edge st n to last_edge st n do + if forward_edge st.flags.{j} then begin + let d' = d + st.weight.{j} in + if !debug then Format.eprintf "!%d %d@." d' (st.dist.{st.target.{j}}); + if d' <= max_dist then begin + let t = st.target.{j} in + if t <> ignored_node && st.dist.{t} > d' then + ( + if !debug then Format.eprintf "> %d %d [%d]@." t d' st.id.{j}; + Heap.insert heap t d' +) + end + end + done + end; + loop () + end + in + loop (); + (!steps, !steps') +(* +Format.eprintf "%d %d (%d) %d - %d -- %d@." !steps !steps' target_count (way_id id) max_dist i; +*) +(* +;(*if !trace = 2 then*) Format.eprintf "%d %d %d %d@." !steps !steps' target_count max_dist +*) +(*;print_int !steps; print_char '\n'*) + +let delete_edges st i j = + let n = ref (last_edge st i) in + let k = ref (first_edge st i) in + while !k <= !n do + if st.target.{!k} = j then begin + move_edge st !n !k; + decr n + end else + incr k + done; +(* +Format.eprintf "deleted: %d <- %d : %d %d @." i j (!n + 1 - st.first_edge.{i}) st.edge_count.{i}; +*) + st.edge_count.{i} <- !n + 1 - st.first_edge.{i} + + +let iter_remaining st f = + for i = 0 to st.remaining_count - 1 do + f i st.remaining.{i} + done + +let redundant_edge st i j dir = + let res = ref false in + let target = st.target.{j} in + for k = first_edge st i to last_edge st i do + res := + !res || + (k <> j && st.target.{k} = target && st.flags.{k} land dir <> 0 && + (st.weight.{k} < st.weight.{j} || + (st.weight.{k} = st.weight.{j} && k < j))) + done; + !res + +let prepare_node_contraction dij_state st i = + let deleted_edges = st.edge_count.{i} in + let original_deleted = ref 0 in + let added_edges = ref 0 in + let original_added_edges = ref 0 in + let shortcuts = ref [] in + for j = first_edge st i to last_edge st i do +(* Format.eprintf "%d %d@." i j;*) + original_deleted := !original_deleted + st.original_edges.{j}; + if backward_edge st.flags.{j} && not (redundant_edge st i j 1) then begin + let max_dist = ref (-1) in + for k = first_edge st i to last_edge st i do +(* +print_char '/'; +print_char ' '; +print_int i; +print_char ' '; +print_char '/'; +print_char ' '; +print_int (min st.target.{j} st.target.{k}); +print_char ' '; +print_char '/'; +print_char ' '; +print_int (max st.target.{j} st.target.{k}); +print_char '\n'; +*) + if + k <> j && forward_edge st.flags.{k} && not (redundant_edge st i k 2) + then + max_dist := max !max_dist (st.weight.{j} + st.weight.{k}) + done; + if !max_dist >= 0 then begin +let (steps, steps') = + dijkstra st dij_state (st.target.{j}) i !max_dist (st.id.{j}) i; +in +let success = ref false in + for k = first_edge st i to last_edge st i do + if + k = j || not (forward_edge st.flags.{k}) || redundant_edge st i k 2 + then () else + if st.dist.{st.target.{k}} <= st.weight.{j} + st.weight.{k} then begin +success := true +(* + Format.eprintf " %d -> %d (%d) %d@." k st.dist.{st.target.{k}} + (st.id.{k}) i; + Format.eprintf "ignored node %d <--%d(%d)-- %d --%d(%d)--> %d@." st.target.{j} (way_id st.id.{j}) st.weight.{j} i (way_id st.id.{k}) st.weight.{k} st.target.{k}; + Format.eprintf "path: %a@." (show_path st) st.target.{k} +*) + end else begin + added_edges := !added_edges + 2; + original_added_edges := + !original_added_edges + + 2 * (st.original_edges.{j} + st.original_edges.{k}); + shortcuts := + (st.target.{j}, st.target.{k}, st.weight.{j} + st.weight.{k}, + 2, st.original_edges.{j} + st.original_edges.{k}) :: !shortcuts + end +(*;Format.printf "%d / %d / %b@." steps steps' !success*) + done + end + end + done; + let p = float st.depth.{i} in + let p = + if deleted_edges > 0 then + p +. 2. *. float !added_edges /. float deleted_edges + else + p + in + let p = + if deleted_edges > 0 then + p +. 4. *. float !original_added_edges /. float !original_deleted + else + p + in +(* +Format.eprintf "%d: %f (%d)@." i p deleted_edges; +*) + (p, !shortcuts) + +let add_edge st i j weight flags original_edges shortcut_node id = +(* +Format.eprintf "add %d -> %d (%d)@." i j weight; +*) + reserve st i 1; + st.edge_count.{i} <- st.edge_count.{i} + 1; + let k = last_edge st i in + st.target.{k} <- j; + st.weight.{k} <- weight; + st.flags.{k} <- flags; + st.original_edges.{k} <- original_edges; + st.shortcut_node.{k} <- shortcut_node; + st.id.{k} <- 0 + +let contract_node dij_state st i = + let (_, shortcuts) = prepare_node_contraction dij_state st i in + for j = first_edge st i to last_edge st i do + let k = st.target.{j} in + delete_edges st k i + done; + (*XXX Could be done at a latter stage? *) + for j = first_edge st i to last_edge st i do + Column.append st.out_source i; + Column.append st.out_target st.target.{j}; + Column.append st.out_weight st.weight.{j}; + Column.append st.out_flags st.flags.{j}; + Column.append st.out_shortcut + (if st.original_edges.{j} > 1 then st.shortcut_node.{j} else -1) + done; + let shortcuts1 = Array.of_list shortcuts in +(* + Array.sort + (fun (j, k, _, _, _) (j', k', _, _, _) -> + let c = compare j j' in if c <> 0 then c else compare k k') + shortcuts1; + for l = 1 to Array.length shortcuts1 - 1 do + let (j, k, w, _, _) = shortcuts1.(l -1) in + let (j', k', w', _, _) = shortcuts1.(l) in + assert (j <> j' || k <> k') + done; +*) + let shortcuts2 = + Array.map (fun (j, k, w, _, orig_edges) -> (k, j, w, 1, orig_edges)) + shortcuts1 + in + let shortcuts = Array.append shortcuts1 shortcuts2 in + Array.sort + (fun (j, k, _, _, _) (j', k', _, _, _) -> + let c = compare j j' in if c <> 0 then c else compare k k') + shortcuts; + let len = Array.length shortcuts in +(* +Format.eprintf "%d@." len; +*) + let l = ref 0 in + while !l < len do + let (j, k, w, f, orig_edges) = shortcuts.(!l) in + if + !l + 1 < len && + let (j', k', w', _, _) = shortcuts.(!l + 1) in +(* +Format.eprintf "%d %d/%d %d/%d %d/%d@." !l j' j k' k w' w; +*) + j' = j && k' = k && w' = w + then begin + add_edge st j k w 3 orig_edges i 0; + incr l + end else + add_edge st j k w f orig_edges i 0; + incr l + done + +let update_neighbours dij_state st i = + let a = Array.make st.edge_count.{i} max_int in + let j0 = first_edge st i in + let depth = 1 + st.depth.{i} in + for j = j0 to last_edge st i do + let k = st.target.{j} in + st.depth.{k} <- max st.depth.{k} depth; + a.(j - j0) <- k + done; + Array.sort compare a; + for j = 0 to Array.length a - 1 do + let k = a.(j) in + if j = 0 || k <> a.(j - 1) then begin + st.priority.{k} <- fst (prepare_node_contraction dij_state st k) +(*;if !trace = 2 then Format.eprintf "%d / %f@." k st.priority.{k}*) + end + done + +let max_depth = ref 0 + +let update_neighbour_depths st i = + let depth = 1 + st.depth.{i} in + for j = first_edge st i to last_edge st i do + if depth > !max_depth then begin max_depth := depth; Format.eprintf "depth: %d@." depth end; + let k = st.target.{j} in + st.depth.{k} <- max st.depth.{k} depth + done + +let update_neighbour_priorities updated dij_state st i = + for j = first_edge st i to last_edge st i do + let k = st.target.{j} in + if not (Bitvect.test updated k) then begin + Bitvect.set updated k; + st.priority.{k} <- fst (prepare_node_contraction dij_state st k) +(* +;if !trace = 2 then Format.eprintf "%d / %f@." k st.priority.{k} +*) + end + done + +let check_edges st = + let rem = Bitvect.make (Array1.dim st.first_edge) in + iter_remaining st (fun _ i -> Bitvect.set rem i); + iter_remaining st + (fun _ i -> + for j = first_edge st i to last_edge st i do + let t = st.target.{j} in + if not (Bitvect.test rem t) then + Format.eprintf "DANGLING EDGE %d -> %d@." i t + else begin + let rev = ref false in + for k = first_edge st t to last_edge st t do + if st.target.{k} = i then rev := true + done; + if not !rev then + Format.eprintf "MISSING REVERSE EDGE %d -> %d@." i t + end + done) diff --git a/osm/lib/osm_display.ml b/osm/lib/osm_display.ml new file mode 100644 index 0000000..1fa95d0 --- /dev/null +++ b/osm/lib/osm_display.ml @@ -0,0 +1,1760 @@ +(* OSM tools + * Copyright (C) 2013 Jérôme Vouillon + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, with linking exception; + * either version 2.1 of the License, or (at your option) any later version. + * + * 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 Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + *) + +(* +* aeroway=runway can be a surface + +* When a leaf bounding box is very large, perform clipping of the + objects it contains, to avoid rendering artifacts (Cairo uses + integers internally, +/- 8 millions, bug #20091) + +* Redesign R-tree of linear features + + mostly constant: + layer 4 bits + bridge/tunnel 2 bit + + category 5 bits + oneway/access 5 bits (car 2, bikes 2, pedestrian 1) + ==> 2 bytes => 6 bits remaining for the number of ways + +* Rendering performance: + * do not render railway=rail;service=* at low zoom levels + * multiple surface R-trees: + R-tree of small surfaces; R-tree of other surfaces; + R-tree with about the 1% largest surfaces; ... + +* Rendering fixes + ==> Highway with area=yes + ==> do not draw them as ways as well + ==> Outline of highway surfaces + ==> Notre Dame de Paris and Sacré Coeur are missing! + ==> use same algorithm as osm2pgsql to deal with multipolygon tags + +* We could improve the rendering of tunnels, ... by classifying nodes: + - do not use round linecap at tunnel extremities + - draw additional circles when a path extremities do not agree + - detect nodes with level mismatch + +* One-way arrows (foot, bicycle, car), accessibility (foot, bicycle, car) + ===> share between ways! +*) + +let async_zoom = true +let async_zoom_in = true +let async_delay = 50 (*ms*) + +let debug_time = true + +let (>>) x f = f x + +(****) + +(*XXX Duplicated code...*) + +let sint_of_int i = let i' = i lsr 1 in if i land 1 = 1 then (-i' - 1) else i' + +let rec read_varint_rec a p v offs = + let i = !p in + let c = Char.code a.[i] in + incr p; + if c >= 0x80 then + read_varint_rec a p (v lor ((c land 0x7f) lsl offs)) (offs + 7) + else + v lor (c lsl offs) + +let read_varint a p = read_varint_rec a p 0 0 + +let read_signed_varint a p = sint_of_int (read_varint a p) + +let read_int_2 s pos = Char.code s.[pos] lor (Char.code s.[pos + 1] lsl 8) + +(****) + +let rec log2 x = if x <= 1 then 0 else 1 + log2 (x lsr 1) +let log2_tbl = Array.init 256 log2 +let log2_16 x = + let x' = x lsr 8 in + if x' = 0 then + Array.unsafe_get log2_tbl x + else + 8 + Array.unsafe_get log2_tbl x' +let log2_32 x = + let x' = x lsr 16 in + if x' = 0 then log2_16 x else 16 + log2_16 x' +let log2 x = + let x' = x lsr 32 in + if x' = 0 then log2_32 x else 32 + log2_32 x' + +(****) + +module Surface = Osm_category.Make (struct + type t = + [ `Water | `Forest | `Grass | `Heath | `Rock | `Sand | `Glacier + | `Farmland | `Residential | `Commercial + | `Industrial | `Park | `Cemetery | `Parking | `Building + | `Highway_residential | `Highway_unclassified | `Highway_living_street + | `Highway_service | `Highway_pedestrian | `Highway_track + | `Highway_footway | `Highway_path ] + let list = + [ `Water; `Forest; `Grass; `Heath; `Rock; `Sand; `Glacier; + `Farmland; `Residential; `Commercial; + `Industrial; `Park; `Cemetery; `Parking; `Building; + `Highway_residential; `Highway_unclassified; `Highway_living_street; + `Highway_service; `Highway_pedestrian; `Highway_track; + `Highway_footway; `Highway_path ] +end) + +module Linear_feature = Osm_category.Make (struct + type t = + [ `Motorway | `Trunk | `Primary | `Secondary | `Tertiary + | `Motorway_link | `Trunk_link | `Primary_link + | `Secondary_link | `Tertiary_link + | `Residential | `Unclassified | `Living_street | `Road | `Service + | `Pedestrian | `Track | `Cycleway | `Bridleway | `Footway | `Path | `Steps + | `River | `Canal | `Stream + | `Runway | `Taxiway + | `Rail | `Tram | `Subway ] + let list = + [ `Motorway; `Trunk; `Primary; `Secondary; `Tertiary; + `Motorway_link; `Trunk_link; `Primary_link; + `Secondary_link; `Tertiary_link; + `Residential; `Unclassified; `Living_street; `Road; `Service; + `Pedestrian; `Track; `Cycleway; `Bridleway; `Footway; `Path; `Steps; + `River; `Canal; `Stream; + `Runway; `Taxiway; + `Rail; `Tram; `Subway ] +end) + +let pedestrian_surface = lazy (Cairo.PNG.create "images/pedestrian.png") + +(****) + +module Bbox = Rtree.Bbox + +type bbox = Bbox.t = + { min_lat : int; max_lat : int; min_lon : int; max_lon : int } + +let bounding_box ratio x_min y_min x_max y_max = + let ratio = float ratio in + let to_lat y = truncate (Osm_geometry.y_to_lat (y *. 10_000_000.) /. ratio) in + let to_lon x = truncate (x *. 10_000_000. /. ratio) in + { min_lat = to_lat y_min; max_lat = to_lat y_max; + min_lon = to_lon x_min; max_lon = to_lon x_max } + +(****) + +(* Old edge R-tree, now used only to find nodes for routing + (should be eventually superseded by the R-tree containing + linear features). *) + +let leaf_size = 1024 + +let decode_leaf leaves i = + let buf = String.create leaf_size in + seek_in leaves (i * leaf_size); + really_input leaves buf 0 leaf_size; + let node_len = read_int_2 buf 0 in + let edge_len = read_int_2 buf 2 in + let node_lat = Array.make (node_len / 2) 0 in + let node_lon = Array.make (node_len / 2) 0 in + let node_idx = Array.make (node_len / 2) 0 in + let edges = Array.make edge_len 0 in + let i = ref 0 in + let pos = ref 4 in + let lat = ref 0 in + let lon = ref 0 in + let idx = ref 0 in + while !pos < node_len + 4 do + let v = !lat + read_signed_varint buf pos in + node_lat.(!i) <- v; + lat := v; + let v = !lon + read_signed_varint buf pos in + node_lon.(!i) <- v; + lon := v; + let v = !idx + read_signed_varint buf pos in + node_idx.(!i) <- v; + idx := v; + incr i + done; + let node_lat = Array.sub node_lat 0 !i in + let node_lon = Array.sub node_lon 0 !i in + let node_idx = Array.sub node_idx 0 !i in + let i = ref 0 in + let node = ref 0 in + while !pos < edge_len + node_len + 4 do + let v = !node + read_signed_varint buf pos in + edges.(!i) <- v; + node := v; + incr i + done; + let edges = Array.sub edges 0 !i in + (node_lat, node_lon, node_idx, edges) + +let clamp x min max = if x < min then min else if x > max then max else x + +let to_deg x = x * 50 + +let distance_to_box bbox lat lon = + let lat' = clamp lat (to_deg bbox.min_lat) (to_deg bbox.max_lat) in + let lon' = clamp lon (to_deg bbox.min_lon) (to_deg bbox.max_lon) in + Osm_geometry.distance lat lon lat' lon' + +let nearest_point = lazy begin + let (leaves, routine_nodes) = + Rtree.open_in (Column.file_in_database "highway/r_tree") in + let leaves = open_in leaves in + Rtree.find_nearest_point routine_nodes distance_to_box + (fun j lat lon -> + let (node_lat, node_lon, node_idx, _) = decode_leaf leaves j in + let i0 = ref 0 in + let dist = ref max_int in + for i = 0 to Array.length node_lat - 1 do + let d = + Osm_geometry.distance lat lon (to_deg node_lat.(i)) (to_deg node_lon.(i)) + in + if d < !dist then begin + dist := d; + i0 := i + end + done; + if !dist = max_int then + None + else + Some (!dist, (node_idx.(!i0), node_lat.(!i0), node_lon.(!i0)))) +end + +(****) + +(*XXX Temp code: rebuild ways *) + +let debug = false + +type t = (bool ref * int) list + +let initialize g = + let g' = Array.make (Array.length g) [] in + Array.iteri + (fun i l -> + List.iter + (fun j -> + let m = ref false in + g'.(i) <- (m, j) :: g'.(i); + g'.(j) <- (m, i) :: g'.(j)) + l) + g; + g' + +let odd_degree g = (ref 0, Array.map (fun l -> (List.length l) land 1 = 1) g) + +let rec next_node g i = + match g.(i) with + [] -> + -1 + | (m, j) :: r -> + g.(i) <- r; + if !m then + next_node g i + else begin + m := true; + j + end + +let rec next_odd (i, o) = + if !i = Array.length o then + -1 + else if o.(!i) then begin + o.(!i) <- false; + !i + end else begin + incr i; + next_odd (i, o) + end + +let rec find graph odd_deg i j cont lst = + let k = next_node graph j in +if debug then Format.eprintf "Going from %d to %d@." j k; + if k = -1 then begin + if i = j then begin +if debug then Format.eprintf "Loop through %d@." i; + (i :: cont, lst) + end else begin +if debug then Format.eprintf "Path %d --> %d@." i j; + (snd odd_deg).(j) <- false; + let k = next_odd odd_deg in +if debug then Format.eprintf "Continuing from %d@." k; + if k = -1 then + ([j], lst) + else begin + let (path, lst) = find graph odd_deg i k cont lst in + ([j], path :: lst) + end + end + end else begin + let (path, lst) = find graph odd_deg i k cont lst in +if debug then Format.eprintf "--@."; + find graph odd_deg j j (path) lst + end + +let rec find_circuits graph odd_deg i lst = + if i = Array.length graph then + lst + else begin + let j = next_node graph i in +if debug then Format.eprintf "Circuit through %d-%d@." i j; + if j = -1 then + find_circuits graph odd_deg (i + 1) lst + else + let (path, lst) = find graph odd_deg i j [] lst in + find_circuits graph odd_deg (i + 1) ((i :: path) :: lst) + end + +let rec find_unclosed_paths graph odd_deg lst = + let i = next_odd odd_deg in + if i <> -1 then begin + let (path, lst) = find graph odd_deg i i [] lst in + find_unclosed_paths graph odd_deg (path :: lst) + end else + lst + +let find_paths graph = + let graph = initialize graph in + let odd_deg = odd_degree graph in +if debug then Format.eprintf "Finding unclosed paths@."; + let lst = find_unclosed_paths graph odd_deg [] in +if debug then Format.eprintf "Finding circuits@."; + find_circuits graph odd_deg 0 lst + +let v = false + +let num = ref 0 +let miss = ref 0 +let find_node tbl i = + incr num; + if not tbl.(i) then begin incr miss; tbl.(i) <- true end + +let chunk tbl edges i l = + let (_, _, _, _, (node_lat, node_lon)) = edges.(i) in + let g = Array.make (Array.length node_lat) [] in + for j = i to i + l - 1 do + let (n1, n2, _, _, _) = edges.(j) in + g.(n1) <- n2 :: g.(n1) + done; + let l = find_paths g in +List.iter (fun p -> List.iter (fun i -> find_node tbl i) p) l; +if v then +List.iter + (fun p -> + Format.eprintf "+ "; + List.iter (fun n -> Format.eprintf "%d " n) p; + Format.eprintf "@.") + l; + l + +let build_paths edges = + Array.sort + (fun (_, _, cat, lay, _) (_, _, cat', lay', _) -> + let c = compare cat cat' in + if c <> 0 then c else + compare lay lay') + edges; + let len = Array.length edges in + let i0 = ref 0 in + let prev_cat = ref (-1) in + let prev_lay = ref (-1) in + let (_, _, _, _, (x, y)) = edges.(0) in + let tbl = Array.make (Array.length x) false in + let paths = ref [] in + for i = 0 to len - 1 do + let (_, _, cat, lay, _) = edges.(i) in + if (cat <> !prev_cat || lay <> !prev_lay) && i > !i0 then begin + paths := + (!prev_cat, !prev_lay, chunk tbl edges !i0 (i - !i0)) :: !paths; + i0 := i + end; + prev_cat := cat; + prev_lay := lay + done; + paths := + (!prev_cat, !prev_lay, chunk tbl edges !i0 (len - !i0)) :: !paths; +if v then Format.eprintf "----@."; + for i = 0 to Array.length tbl - 1 do + assert (tbl.(i)) + done; + Array.of_list + (List.map + (fun (cat, lay_info, l) -> + (* Extract layer and perform sign extension *) + let lay = (lay_info lsr 2) lxor 8 - 8 in + ((cat, lay, lay_info land 1 = 1, lay_info land 2 = 2), + List.map + (fun p -> + Array.of_list (List.map (fun i -> x.(i)) p), + Array.of_list (List.map (fun i -> y.(i)) p)) + l)) + !paths) + +(* R-tree containing linear features *) + +let linear_leaf_read = ref 0 + +let linear_ratio = 50 + +let decode_leaf ratio leaves i = + incr linear_leaf_read; + let leaf_size = 2048 in + let buf = String.create leaf_size in + seek_in leaves (i * leaf_size); + really_input leaves buf 0 leaf_size; + + let node_len = read_int_2 buf 0 in + let x = Array.make (node_len / 2) 0. in + let y = Array.make (node_len / 2) 0. in + let i = ref 0 in + let pos = ref 4 in + let lat = ref 0 in + let lon = ref 0 in + while !pos < node_len + 4 do + let v = !lat + read_signed_varint buf pos in + y.(!i) <- Osm_geometry.lat_to_y (float (v * linear_ratio)); + lat := v; + let v = !lon + read_signed_varint buf pos in + x.(!i) <- float (v * linear_ratio); + lon := v; + incr i + done; + let x = Array.sub x 0 !i in + let y = Array.sub y 0 !i in + let nodes = (x, y) in + + let edge_len = read_int_2 buf 2 in + let edges = Array.make edge_len (0, 0, 0, 0, nodes) in + let i = ref 0 in + let node = ref 0 in + while !pos < edge_len + node_len + 4 do + let n1 = !node + read_signed_varint buf pos in + node := n1; + let n2 = !node + read_signed_varint buf pos in + node := n2; + let cat = Char.code buf.[!pos] in + let layer = Char.code buf.[!pos + 1] in + pos := !pos + 2; + edges.(!i) <- (n1, n2, cat, layer, nodes); + incr i + done; + Array.sub edges 0 !i + +let cache = Lru_cache.make 1000 + +let decode_leaf ratio leaves = + Lru_cache.funct cache + (fun i -> build_paths (decode_leaf linear_ratio leaves i)) + +let open_tree name = + let ratio = linear_ratio in + let (leaves, tree) = Rtree.open_in (Column.file_in_database name) in + let leaves = open_in leaves in + (linear_ratio, decode_leaf ratio leaves, tree) + +let rtrees = lazy begin + [((-1., 11.5), open_tree "linear/rtrees/large_3"); + ((11.5, 12.5), open_tree "linear/rtrees/large_2"); + ((12.5, 13.5), open_tree "linear/rtrees/large_1"); + ((13.5, 30.), open_tree "linear/rtrees/all")] +end + +let find_linear_features level x_min y_min x_max y_max = + let lst = ref [] in + List.iter + (fun ((min_level, max_level), (ratio, decode, tree)) -> + if level > min_level && level <= max_level then begin + let bbox = bounding_box ratio x_min y_min x_max y_max in + Rtree.find tree bbox (fun i -> lst := decode i :: !lst) + end) + (Lazy.force rtrees); + Array.concat !lst + +(****) + +let coastline_leaf_size = 2048 + +let decode_coastline ratio leaves i = + let buf = String.create coastline_leaf_size in + seek_in leaves (i * coastline_leaf_size); + really_input leaves buf 0 coastline_leaf_size; + let n = read_int_2 buf 0 in + let pos = ref (2 + 2 * n) in + let lat = ref 0 in + let lon = ref 0 in + let ways = ref [] in + for i = 0 to n - 1 do + let l = read_int_2 buf (2 + 2 * i) in + let x = Array.make (l + 1) 0. in + let y = Array.make (l + 1) 0. in + for j = 0 to l - 1 do + lat := !lat + read_signed_varint buf pos; + lon := !lon + read_signed_varint buf pos; +(*if j = 0 then Format.eprintf "%d %d@." !lon !lat;*) + x.(j) <- float (!lon * ratio); + y.(j) <- Osm_geometry.lat_to_y (float (!lat * ratio)); + done; + x.(l) <- x.(0); + y.(l) <- y.(0); + ways := (x, y) :: !ways + done; + Array.of_list !ways + +let decode_coastline ratio leaves = + Lru_cache.funct cache + (fun i -> decode_coastline ratio leaves i) + +let open_tree name = + let ch = open_in (Column.file_in_database (Filename.concat name "ratio")) in + let ratio = int_of_string (input_line ch) in + close_in ch; + let (leaves, tree) = Rtree.open_in (Column.file_in_database name) in + let leaves = open_in leaves in + (ratio, decode_coastline ratio leaves, tree) + +let rtrees = lazy begin + if Sys.file_exists (Column.file_in_database "coastline/rtrees/2") then + [((-1., 2.), open_tree "coastline/rtrees/2"); + ((2., 4.), open_tree "coastline/rtrees/4"); + ((4., 6.), open_tree "coastline/rtrees/6"); + ((6., 8.), open_tree "coastline/rtrees/8"); + ((8., 30.), open_tree "coastline/rtrees/small")] + else + [] +end + +let find_coastline level x_min y_min x_max y_max = + let lst = ref [] in + List.iter + (fun ((min_level, max_level), (ratio, decode, tree)) -> + if level > min_level && level <= max_level then begin + let bbox = bounding_box ratio x_min y_min x_max y_max in + Rtree.find tree bbox + (fun i -> lst := decode i :: !lst) + end) + (Lazy.force rtrees); + Array.concat !lst + +(****) + +(* R-tree containing surfaces *) + +let surface_leaf_size = 2048 + +let surface_leaf_read = ref 0 + +let decode_surfaces ratio leaves i = + let buf = String.create surface_leaf_size in + seek_in leaves (i * surface_leaf_size); + really_input leaves buf 0 surface_leaf_size; + let len = read_int_2 buf 0 in + surface_leaf_read := !surface_leaf_read + len; + let buf = + if len > 1 then begin + let buf' = String.create (surface_leaf_size * len) in + String.blit buf 0 buf' 0 surface_leaf_size; + really_input leaves buf' surface_leaf_size + ((len - 1) * surface_leaf_size); + buf' + end else + buf + in + let n = read_int_2 buf 2 in + let pos = ref (4 + 4 * n) in + let lat = ref 0 in + let lon = ref 0 in + let ways = ref [] in + let category = ref 0 in + let layer = ref 0 in + let lst = ref [] in + for i = 0 to n - 1 do + let l = read_int_2 buf (4 + 4 * i) in + let cat = Char.code buf.[4 + 4 * i + 2] in + let lay = Char.code buf.[4 + 4 * i + 3] - 128 in + if cat <> 0 then begin + if !ways <> [] then lst := (!category, !layer, List.rev !ways) :: !lst; + category := cat; + layer := lay; + ways := [] + end; + let x = Array.make (l + 1) 0. in + let y = Array.make (l + 1) 0. in + for j = 0 to l - 1 do + lat := !lat + read_signed_varint buf pos; + lon := !lon + read_signed_varint buf pos; + x.(j) <- float (!lon * ratio); + y.(j) <- Osm_geometry.lat_to_y (float (!lat * ratio)); + done; + x.(l) <- x.(0); + y.(l) <- y.(0); + ways := (x, y) :: !ways + done; + if !ways <> [] then lst := (!category, !layer, List.rev !ways) :: !lst; + !lst + +let prepare_surfaces lst = + Array.of_list + (List.map + (fun (cat, layer, ways) -> + let area = + List.fold_left + (fun a (x, y) -> a +. Osm_geometry.polygon_area_float x y) 0. ways in + (layer, truncate (area +. 0.5), cat, ways)) + lst) + +let decode_surfaces ratio leaves = + Lru_cache.funct cache + (fun i -> prepare_surfaces (decode_surfaces ratio leaves i)) + +let open_tree name = + let ch = open_in (Column.file_in_database (Filename.concat name "ratio")) in + let ratio = int_of_string (input_line ch) in + close_in ch; + let (leaves, tree) = Rtree.open_in (Column.file_in_database name) in + let leaves = open_in leaves in + (ratio, decode_surfaces ratio leaves, tree) + +let large_surfaces = lazy (open_tree "surfaces/rtrees/large") + +let rtrees = lazy begin + [((-1., 6.), open_tree "surfaces/rtrees/06"); + ((6., 7.), open_tree "surfaces/rtrees/07"); + ((7., 8.), open_tree "surfaces/rtrees/08"); + ((8., 9.), open_tree "surfaces/rtrees/09"); + ((9., 10.), open_tree "surfaces/rtrees/10"); + ((10., 12.), open_tree "surfaces/rtrees/12"); + ((12., 30.), Lazy.force large_surfaces); + ((15.5, 30.), open_tree "surfaces/rtrees/small")] +end + +let find_surfaces level x_min y_min x_max y_max = + let lst = ref [] in + List.iter + (fun ((min_level, max_level), (ratio, decode, tree)) -> + if level > min_level && level <= max_level then begin + let bbox = bounding_box ratio x_min y_min x_max y_max in + Rtree.find tree bbox + (fun i -> lst := decode i :: !lst) + end) + (Lazy.force rtrees); + Array.concat !lst + +(**** Pixmap ***) + +type rectangle = { x : int; y : int; width : int; height: int } + +let empty_rectangle = {x = 0; y = 0; width = 0; height = 0} +let rectangle_is_empty r = r.width = 0 || r.height = 0 + +type surface = + { mutable surface : Cairo.Surface.t option; + mutable p_width : int; mutable p_height : int; + mutable valid_rect : rectangle } + +let make_surface () = + { surface = None; p_width = 0; p_height = 0; + valid_rect = empty_rectangle } + +let invalidate_surface p = p.valid_rect <- empty_rectangle + +let grow_surface pm window width height = + let width = max width pm.p_width in + let height = max height pm.p_height in + if width > pm.p_width || height > pm.p_height then begin + let old_p = pm.surface in +(* + let p = GDraw.pixmap ~width ~height ~window () in +*) + let p = + Cairo.Surface.create_similar + (Cairo.get_target (Cairo_gtk.create window#misc#window)) + Cairo.COLOR_ALPHA ~width ~height + in + let r = pm.valid_rect in + begin match old_p with + Some old_p -> + let ctx = Cairo.create p in + Cairo.set_source_surface ctx old_p 0. 0.; + Cairo.rectangle ctx 0. 0. (float r.width) (float r.height); + Cairo.set_operator ctx Cairo.SOURCE; + Cairo.fill ctx +(* + p#put_pixmap ~x:0 ~y:0 ~xsrc:0 ~ysrc:0 + ~width:r.width ~height:r.height old_p#pixmap +*) + | None -> + () + end; + pm.surface <- Some p; + pm.p_width <- width; + pm.p_height <- height + end + +let get_surface pm = match pm.surface with Some p -> p | None -> assert false + +(**** Global state ***) + +type state = + { mutable rect : rectangle; + mutable level : float; + mutable prev_rect : rectangle; + mutable prev_level : float; + mutable active : bool; + mutable timeout : Glib.Timeout.id option; + surface : surface; + mutable marker1 : (int * float * float) option; + mutable marker2 : (int * float * float) option; + mutable path : (int * int) list } + +let compute_scale st = 256. /. 360. *. 2. ** st.level + +(**** Routing ****) + +let find_marker st x y = + let scale = 256. /. 360. *. 2. ** st.level in + let x' = (float st.rect.x +. x) /. scale in + let y' = -. (float st.rect.y +. y) /. scale in + let lat = truncate (Osm_geometry.y_to_lat (y' *. 10_000_000.)) in + let lon = truncate ((x' *. 10_000_000.)) in + Format.eprintf "%d %d@." lat lon; + let (d, (i, lat, lon)) = Lazy.force nearest_point lat lon in + let lat = float lat /. 200000. in + let lon = float lon /. 200000. in + Format.eprintf "%d: %f - %f %f@." i (float d /. 1000.) lat lon; + Some (i, lat, lon) + +let routing = lazy (Routing.init ()) +let node_lat = lazy (Column.open_in (Column.named "highway/sorted_node" "lat")) +let node_lon = lazy (Column.open_in (Column.named "highway/sorted_node" "lon")) + +let update_route st = + begin match st.marker1, st.marker2 with + Some (i1, _, _), Some (i2, _, _) -> + let l = Routing.find (Lazy.force routing) i1 i2 in + st.path <- + List.map (fun n -> (Column.get (Lazy.force node_lat) n, Column.get (Lazy.force node_lon) n)) l + | _ -> + () + end + +(****) + +let draw_linear_features st ctx pred stroke i = + let scale = compute_scale st in + let prev_info = ref (-1, 0, false, false) in + let count = ref 0 in + i (fun (info, ways) -> + if pred info then begin + if info <> !prev_info then + if !count > 0 then stroke !prev_info; + prev_info := info; + List.iter + (fun (x, y) -> + let len = Array.length x in + if !count > 0 && !count + len > 10000 then begin + stroke !prev_info; + count := 0 + end; + let coeff = scale /. 10_000_000. in + if st.level >= 15. && Array.length x > 2 then begin + (*XXX This could be precomputed when decoding the path *) + let ((x, y), (x1, y1), (x2, y2)) = + Line_smoothing.perform x y in + let len = Array.length x in + count := !count + 3 * len; + Cairo.move_to ctx (x.(0) *. coeff) (y.(0) *. coeff); + for k = 1 to len - 1 do + Cairo.curve_to ctx + (x1.(k - 1) *. coeff) (y1.(k - 1) *. coeff) + (x2.(k - 1) *. coeff) (y2.(k - 1) *. coeff) + (x.(k) *. coeff) (y.(k) *. coeff) + done + end else begin + count := !count + len; + Cairo.move_to ctx (x.(0) *. coeff) (y.(0) *. coeff); + for k = 1 to len - 1 do + Cairo.line_to ctx (x.(k) *. coeff) (y.(k) *. coeff) + done + end) + ways + end); + if !count > 0 then stroke !prev_info + +let draw_surfaces st ctx pred fill i = + let scale = compute_scale st in + let prev_cat = ref (-1) in + let count = ref 0 in + i + (fun ((_, area, cat, ways) as info) -> + if pred info then begin + if (cat <> !prev_cat || !count > 10000) && !prev_cat <> -1 then begin + fill !prev_cat; + count := 0; + end; + let coeff = scale /. 10_000_000. in + List.iter + (fun (x, y) -> + Cairo.move_to ctx (x.(0) *. coeff) (y.(0) *. coeff); + for i = 1 to Array.length x - 2 do + Cairo.line_to ctx (x.(i) *. coeff) (y.(i) *. coeff) + done; + Cairo.Path.close ctx; + count := !count + Array.length x) + ways; + if cat <> !prev_cat then prev_cat := cat + end); + if !prev_cat <> -1 then fill !prev_cat + +let set_surface_color ctx cat = + match cat with + `Water -> + Cairo.set_source_rgb ctx 0.52 0.94 0.94 + | `Forest -> + Cairo.set_source_rgb ctx 0.1 0.7 0.2 + | `Grass -> + Cairo.set_source_rgb ctx 0.3 0.9 0.3 + | `Heath -> + Cairo.set_source_rgb ctx 0.59 0.74 0.42 + | `Rock -> + Cairo.set_source_rgb ctx 0.37 0.42 0.49 + | `Sand -> + Cairo.set_source_rgb ctx 0.94 0.93 0.22 + | `Glacier -> + Cairo.set_source_rgb ctx 0.80 0.94 0.87 + | `Farmland -> + Cairo.set_source_rgb ctx 0.69 0.94 0.27 + | `Park -> + Cairo.set_source_rgb ctx 0.6 1.0 0.6 + | `Residential -> + Cairo.set_source_rgb ctx 0.91 0.94 0.94 + | `Commercial -> + Cairo.set_source_rgb ctx 0.94 0.78 0.78 + | `Industrial -> + Cairo.set_source_rgb ctx 0.87 0.82 0.85 + | `Parking -> + Cairo.set_source_rgb ctx 0.97 0.94 0.72 + | `Cemetery -> + Cairo.set_source_rgb ctx 0.67 0.8 0.69 + | `Building -> + Cairo.set_source_rgb ctx 0.7 0.7 0.7 + | `Highway_pedestrian | `Highway_track + | `Highway_footway | `Highway_path -> + Cairo.set_source_surface ctx (Lazy.force pedestrian_surface) 0. 0.; + Cairo.Pattern.set_extend (Cairo.get_source ctx) Cairo.Pattern.REPEAT + | `Highway_residential | `Highway_unclassified + | `Highway_living_street | `Highway_service -> + Cairo.set_source_rgb ctx 0.8 0.8 0.8 + +(****) + +let draw_coastline st ctx coastline x_min y_min x_max y_max = + let scale = compute_scale st in + let coeff = scale /. 10_000_000. in + Array.iter + (fun (x, y) -> + Cairo.move_to ctx (x.(0) *. coeff) (y.(0) *. coeff); + for i = 1 to Array.length x - 2 do + Cairo.line_to ctx (x.(i) *. coeff) (y.(i) *. coeff) + done; + Cairo.Path.close ctx) + coastline; + Cairo.set_source_rgb ctx 0.52 0.94 0.94; + Cairo.fill ctx + +(****) + +let draw_map_high_levels st ctx surfaces linear_features = + let scale = compute_scale st in + let module SP = Surface.Partition in + let partition = SP.make () in + let landuse = + SP.add_group partition + [`Forest; `Grass; `Heath; `Rock; `Sand; `Glacier; + `Farmland; `Residential; `Commercial; + `Industrial; `Park; `Cemetery; `Parking ] + in + let water = SP.add_group partition [`Water] in + let building_or_pedestrian = + SP.add_group partition + [`Building; + `Highway_pedestrian; `Highway_track; `Highway_footway; `Highway_path ] + in + let highway_areas = + SP.add_group partition + [`Highway_residential; `Highway_unclassified; + `Highway_living_street; `Highway_service ] + in +let t = Unix.gettimeofday () in + let area = Array.map (fun (_, a, _, _) -> - log2 a) surfaces in + let layer = Array.map (fun (l, _, _, _) -> l) surfaces in +let surfaces' = surfaces in + let surfaces = + SP.apply partition surfaces (fun (_, _, cat, _) -> cat) + >> SP.order_totally + >> SP.order_by area + >> SP.order_by layer + >> SP.order_by_group + >> SP.select + in +if debug_time then +Format.eprintf "Sorting surfaces (%d elements): %.3f@." (Array.length surfaces') (Unix.gettimeofday () -. t); + +let t = Unix.gettimeofday () in + let module LP = Linear_feature.Partition in + let partition = LP.make () in + let waterway = LP.add_group partition [`River; `Canal; `Stream] in + let ways = + LP.add_group partition + [ `Runway; `Taxiway; + `Rail; `Tram; `Subway; + `Pedestrian; `Track; `Cycleway; `Bridleway; `Footway; `Path; `Steps; + `Residential; `Unclassified; `Living_street; `Road; `Service; + `Tertiary_link; `Secondary_link; `Primary_link; + `Tertiary; `Secondary; `Primary; + `Trunk_link; `Motorway_link; + `Trunk; `Motorway ] + in + let layer = + Array.map + (fun ((cat, lay, is_bridge, is_tunnel), _) -> + 16 + + lay * 3 + (if is_bridge then 1 else 0) + - (if is_tunnel then 1 else 0)) + linear_features + in +let linear_features' = linear_features in + let linear_features = + LP.apply partition linear_features (fun ((cat, _, _, _), _) -> cat) + >> LP.order_totally + >> LP.order_by layer + >> LP.order_by_group + >> LP.select + in +if debug_time then +Format.eprintf "Sorting lines (%d elements): %.3f@." (Array.length linear_features') (Unix.gettimeofday () -. t); +if debug_time then begin +let n = ref 0 in +Array.iter (fun (_, _, _, ways) -> List.iter (fun (x, _) -> n := !n + Array.length x) ways) surfaces'; +Format.eprintf "Surfaces: %d nodes@." !n; + let small_area = truncate (16. *. (*64. *.*) (10_000_000. /. scale) ** 2.) in + let filter (_, area, cat, ways) = + st.level >= 15.5 || + (area > small_area && + (area > 50_000_000 || Surface.of_id cat <> `Building)) + in +let n = ref 0 in +let m = ref 0 in +Array.iter (fun ((_, _, _, ways) as info) -> if filter info then begin incr m; List.iter (fun (x, _) -> n := !n + Array.length x) ways end) surfaces'; +Format.eprintf "Surfaces: %d nodes, %d elements (%.2f%%)@." !n !m (100. *. float !m /. float (Array.length surfaces')); +let n = ref 0 in +Array.iter (fun (_, ways) -> List.iter (fun (x, _) -> n := !n + Array.length x) ways) linear_features'; +Format.eprintf "Lines: %d nodes@." !n +end; + let draw_water_lines layer = + Cairo.set_line_join ctx Cairo.JOIN_ROUND; + begin match layer with + `Tunnel -> + Cairo.set_dash ctx [|2.; 4.|] + | `Bridge -> + Cairo.set_line_cap ctx Cairo.BUTT; + Cairo.set_source_rgb ctx 0. 0. 0.; + Cairo.set_line_width ctx 4.; + linear_features >> LP.with_group waterway >> LP.iter + >> draw_linear_features st ctx + (fun (_, _, is_bridge, is_tunnel) -> is_bridge) + (fun _ -> Cairo.stroke ctx) + | `Ground -> + () + end; + Cairo.set_line_cap ctx Cairo.ROUND; + Cairo.set_source_rgb ctx 0.52 0.94 0.94; + Cairo.set_line_width ctx 2.; + linear_features >> LP.with_group waterway >> LP.iter + >> draw_linear_features st ctx + (fun (_, _, is_bridge, is_tunnel) -> + match layer with + `Tunnel -> is_tunnel + | `Bridge -> is_bridge + | `Ground -> not (is_tunnel || is_bridge)) + (fun _ -> Cairo.stroke ctx); + Cairo.set_dash ctx [||] + in + let draw_casing round i = + Cairo.set_line_cap ctx (if round then Cairo.ROUND else Cairo.BUTT); + Cairo.set_line_join ctx Cairo.JOIN_ROUND; + Cairo.set_dash ctx [||]; + let stroke (cat, _, _, is_tunnel) = + let cat = Linear_feature.of_id cat in + begin match cat with + `Pedestrian | `Track | `Cycleway | `Bridleway | `Footway + | `Path | `Steps | `Rail | `Tram | `Subway -> + Cairo.set_source_rgb ctx 1. 1. 1. + | _ -> + Cairo.set_source_rgb ctx 0. 0. 0. + end; + let line_width = + match cat with + `Trunk | `Motorway -> 11. + | `Trunk_link | `Motorway_link -> 6.5 + | `Tertiary | `Secondary | `Primary -> 8. + | `Tertiary_link | `Secondary_link | `Primary_link -> 5.5 + | `Residential | `Unclassified | `Living_street | `Road -> 6. + | `Service -> 4. + | `Pedestrian | `Track | `Cycleway + | `Bridleway | `Footway | `Path | `Steps -> 4. + | `Rail -> 5. + | `Tram | `Subway -> 6. + | _ -> assert false + in + let line_width = if round then line_width else line_width -. 0.5 in + Cairo.set_line_width ctx line_width; + Cairo.stroke ctx; + in + draw_linear_features st ctx + (fun (cat, _, is_bridge, _) -> + match Linear_feature.of_id cat with + `Motorway | `Trunk | `Primary | `Secondary | `Tertiary + | `Motorway_link | `Trunk_link | `Primary_link + | `Secondary_link | `Tertiary_link + | `Residential | `Unclassified | `Living_street | `Road + | `Service -> + true + | `Pedestrian | `Track | `Cycleway | `Bridleway | `Footway + | `Path | `Steps | `Rail | `Tram | `Subway when is_bridge -> + true + | _ -> + false) + stroke i + in + let draw_inline i = + Cairo.set_line_join ctx Cairo.JOIN_ROUND; + Cairo.set_line_cap ctx Cairo.ROUND; + let stroke (cat, _, _, is_tunnel) = + let cat = Linear_feature.of_id cat in + if cat = `Rail then begin + Cairo.set_line_cap ctx Cairo.BUTT; + Cairo.set_dash ctx [|1.; 3.|]; + Cairo.set_source_rgb ctx 0.27 0.27 0.27; + Cairo.set_line_width ctx 3.; + Cairo.stroke_preserve ctx; + Cairo.set_line_cap ctx Cairo.ROUND; + Cairo.set_dash ctx [||] + end; + begin match cat with + `Motorway | `Trunk | `Motorway_link | `Trunk_link -> + Cairo.set_source_rgb ctx 1. 0.8 0. + | `Pedestrian | `Track | `Cycleway | `Bridleway + | `Footway | `Path | `Steps -> + Cairo.set_source_rgb ctx 0. 0. 0.; + Cairo.set_line_cap ctx Cairo.BUTT; + if cat = `Steps then + Cairo.set_dash ctx [|1.; 2.|] + else + Cairo.set_dash ctx [|2.; 3.|] + | `Residential | `Unclassified | `Living_street | `Road | `Service -> + Cairo.set_source_rgb ctx 0.8 0.8 0.8 + | `Rail | `Tram -> + Cairo.set_source_rgb ctx 0.27 0.27 0.27 + | `Subway -> + Cairo.set_source_rgb ctx 0.6 0.6 0.6 + | `Runway | `Taxiway -> + Cairo.set_source_rgb ctx 0.73 0.73 0.8 + | _ -> + Cairo.set_source_rgb ctx 1. 1. 1. + end; + let line_width = + match cat with + `Motorway | `Trunk -> 6. + | `Motorway_link | `Trunk_link -> 3. + | `Primary | `Secondary | `Tertiary -> 5. + | `Primary_link | `Secondary_link | `Tertiary_link -> 2.5 + | `Residential | `Unclassified | `Living_street | `Road -> 4. + | `Service -> 2.5 + | `Rail -> 1. + | `Tram | `Subway -> 2. + | `Steps -> 3. + | `Pedestrian | `Track | `Cycleway + | `Bridleway | `Footway | `Path -> 1. + | `Runway -> max (2.5e-4 *. scale) 2. + | _ -> 2. + in + Cairo.set_line_width ctx line_width; + Cairo.stroke ctx; + begin match cat with + `Pedestrian | `Track | `Cycleway | `Bridleway + | `Footway | `Path | `Steps -> + Cairo.set_line_cap ctx Cairo.ROUND; + Cairo.set_dash ctx [||] + | _ -> + () + end + in + draw_linear_features st ctx + (fun (cat, _, _, is_tunnel) -> + not is_tunnel || Linear_feature.of_id cat <> `Subway) + stroke i + in + +let t = Unix.gettimeofday () in + (* Draw surfaces *) + let fill_surface cat = + let cat = Surface.of_id cat in + set_surface_color ctx cat; + if st.level >= 17. && cat = `Building then begin + Cairo.set_source_rgb ctx 0.71 0.71 0.71; + Cairo.fill_preserve ctx; + Cairo.set_source_rgb ctx 0.6 0.6 0.6; + Cairo.set_line_width ctx 1.; + Cairo.stroke ctx + end else + Cairo.fill ctx + in + let small_area = truncate (32. *. (10_000_000. /. scale) ** 2.) in + let filter (_, area, cat, ways) = + st.level >= 15.5 || + (area > small_area && + (area > 50_000_000 || Surface.of_id cat <> `Building)) + in + surfaces >> SP.with_group landuse >> SP.iter + >> draw_surfaces st ctx filter fill_surface; + (* Draw water lines (below buildings) *) + draw_water_lines `Ground; + surfaces >> SP.with_group water >> SP.iter + >> draw_surfaces st ctx filter fill_surface; + surfaces >> SP.with_group building_or_pedestrian >> SP.iter + >> draw_surfaces st ctx filter fill_surface; + +if debug_time then +Format.eprintf "Drawing surfaces: %.3f@." (Unix.gettimeofday () -. t); + +let t = Unix.gettimeofday () in + (* Draw linear features *) + let ways = linear_features >> LP.with_group ways in + (* Underground features *) + Cairo.Group.push ctx; + draw_water_lines `Tunnel; + ways >> LP.iter_by_key + (fun layer i -> if layer mod 3 = 0 then draw_casing true (LP.iter i)); + ways >> LP.iter_by_key + (fun layer i -> + if layer mod 3 = 0 then begin + draw_casing false (LP.iter i); draw_inline (LP.iter i) + end); + Cairo.Group.pop_to_source ctx; + Cairo.paint ~alpha:0.4 ctx; + (* Outline *) + ways >> LP.iter_by_key + (fun layer i -> + if layer mod 3 <> 0 then draw_casing true (LP.iter i)); + (* Casing/inline *) + surfaces >> SP.with_group highway_areas >> SP.iter + >> draw_surfaces st ctx filter fill_surface; + ways >> LP.iter_by_key + (fun layer i -> + if layer mod 3 <> 0 then begin + if layer mod 3 = 2 then draw_casing false (LP.iter i); + draw_inline (LP.iter i) + end); + draw_water_lines `Bridge; +if debug_time then +Format.eprintf "Drawing lines: %.3f@." (Unix.gettimeofday () -. t) + +let draw_map_medium_levels st ctx surfaces linear_features = + let scale = compute_scale st in + let module SP = Surface.Partition in + let partition = SP.make () in + let landuse = + SP.add_group partition + [`Forest; `Grass; `Heath; `Rock; `Sand; `Glacier; + `Farmland; `Residential; `Commercial; + `Industrial; `Park; `Cemetery; `Parking ] + in + let water = SP.add_group partition [`Water] in + let building_or_pedestrian = + SP.add_group partition + [`Building; + `Highway_pedestrian; `Highway_track; `Highway_footway; `Highway_path ] + in + let highway_areas = + SP.add_group partition + [`Highway_residential; `Highway_unclassified; + `Highway_living_street; `Highway_service ] + in +let t = Unix.gettimeofday () in + let area = Array.map (fun (_, a, _, _) -> - log2 a) surfaces in + let layer = Array.map (fun (l, _, _, _) -> l) surfaces in +let surfaces' = surfaces in + let surfaces = + SP.apply partition surfaces (fun (_, _, cat, _) -> cat) + >> SP.order_totally + >> SP.order_by area + >> SP.order_by layer + >> SP.order_by_group + >> SP.select + in +if debug_time then +Format.eprintf "Sorting surfaces (%d elements): %.3f@." (Array.length surfaces') (Unix.gettimeofday () -. t); + +let t = Unix.gettimeofday () in + let module LP = Linear_feature.Partition in + let partition = LP.make () in + let waterway = LP.add_group partition [`River; `Canal] in + let ignored_ways = + LP.add_group partition [ `Footway; `Steps; `Service; `Stream ] in + let minor_roads = + LP.add_group partition + [ `Runway; `Taxiway; + `Rail; `Tram; `Subway; + `Pedestrian; `Track; `Cycleway; `Bridleway; `Path; + `Residential; `Unclassified; `Living_street; `Road ] + in + let major_road_links = + LP.add_group partition + [ `Tertiary_link; `Secondary_link; `Primary_link ] in + let major_roads = + LP.add_group partition [ `Tertiary; `Secondary; `Primary ] in + let highway_links = + LP.add_group partition [ `Trunk_link; `Motorway_link ] in + let highways = LP.add_group partition [ `Trunk; `Motorway ] in +let linear_features' = linear_features in + let linear_features = + LP.apply partition linear_features (fun ((cat, _, _, _), _) -> cat) + >> LP.order_totally + >> LP.order_by_group + >> LP.select + in +if debug_time then +Format.eprintf "Sorting lines (%d elements): %.3f@." (Array.length linear_features') (Unix.gettimeofday () -. t); +if debug_time then begin +let n = ref 0 in +Array.iter (fun (_, _, _, ways) -> List.iter (fun (x, _) -> n := !n + Array.length x) ways) surfaces'; +Format.eprintf "Surfaces: %d nodes@." !n; + let small_area = truncate (32. *. (10_000_000. /. scale) ** 2.) in + let filter (_, area, cat, ways) = + (area > small_area && + (area > 50_000_000 || Surface.of_id cat <> `Building)) + in +let n = ref 0 in +let m = ref 0 in +Array.iter (fun ((_, _, _, ways) as info) -> if filter info then begin incr m; List.iter (fun (x, _) -> n := !n + Array.length x) ways end) surfaces'; +Format.eprintf "Surfaces: %d nodes, %d elements (%.2f%%)@." !n !m (100. *. float !m /. float (Array.length surfaces')); +let n = ref 0 in +Array.iter (fun (_, ways) -> List.iter (fun (x, _) -> n := !n + Array.length x) ways) linear_features'; +Format.eprintf "Lines: %d nodes@." !n +end; + let draw_water_lines () = + Cairo.set_line_cap ctx Cairo.BUTT; + Cairo.set_line_join ctx Cairo.JOIN_MITER; + Cairo.set_miter_limit ctx 1.4; + Cairo.set_source_rgb ctx 0.52 0.94 0.94; + Cairo.set_line_width ctx 2.; + linear_features >> LP.with_group waterway >> LP.iter + >> draw_linear_features st ctx + (fun (_, _, is_bridge, is_tunnel) -> not is_tunnel) + (fun _ -> Cairo.stroke ctx); + Cairo.set_dash ctx [||] + in + let draw_casing i = + Cairo.set_line_cap ctx Cairo.BUTT; + Cairo.set_line_join ctx Cairo.JOIN_MITER; + Cairo.set_miter_limit ctx 1.4; + Cairo.set_dash ctx [||]; + Cairo.set_source_rgb ctx 0. 0. 0.; + let stroke (cat, _, _, is_tunnel) = + let cat = Linear_feature.of_id cat in + let line_width = + match cat with + `Trunk | `Motorway -> 8. + | `Motorway_link | `Trunk_link + | `Primary_link | `Secondary_link | `Tertiary_link -> 3. + | _ -> 5. + in + Cairo.set_line_width ctx line_width; + Cairo.stroke ctx; + in + draw_linear_features st ctx + (fun (cat, _, is_bridge, _) -> + match Linear_feature.of_id cat with + `Motorway | `Trunk | `Primary | `Secondary | `Tertiary + | `Motorway_link | `Trunk_link + | `Primary_link | `Secondary_link | `Tertiary_link + | `Residential | `Unclassified | `Living_street | `Road -> + true + | _ -> + false) + stroke i + in + let draw_inline i = + Cairo.set_line_cap ctx Cairo.BUTT; + Cairo.set_line_join ctx Cairo.JOIN_MITER; + Cairo.set_miter_limit ctx 1.4; + let stroke (cat, _, _, is_tunnel) = + let cat = Linear_feature.of_id cat in + begin match cat with + `Motorway | `Trunk | `Motorway_link | `Trunk_link -> + Cairo.set_source_rgb ctx 1. 0.8 0. + | `Pedestrian | `Track | `Path | `Cycleway | `Bridleway -> + Cairo.set_source_rgb ctx 0. 0. 0.; + Cairo.set_line_cap ctx Cairo.BUTT; + Cairo.set_dash ctx [|2.; 3.|] + | `Residential | `Unclassified | `Living_street | `Road -> + Cairo.set_source_rgb ctx 0.7 0.7 0.7 + | `Rail | `Tram -> + Cairo.set_source_rgb ctx 0.27 0.27 0.27 + | `Subway -> + Cairo.set_source_rgb ctx 0.6 0.6 0.6 + | `Runway | `Taxiway -> + Cairo.set_source_rgb ctx 0.73 0.73 0.8 + | _ -> + Cairo.set_source_rgb ctx 1. 1. 1. + end; + let line_width = + match cat with + `Motorway | `Trunk -> 4. + | `Motorway_link | `Trunk_link + | `Primary_link | `Secondary_link | `Tertiary_link -> 1. + | `Primary | `Secondary | `Tertiary -> 3. + | `Residential | `Unclassified | `Living_street | `Road -> + 1.5 + | `Rail -> 1. + | `Tram | `Subway -> 2. + | `Pedestrian | `Track | `Cycleway | `Path | `Bridleway -> 1. + | `Runway -> max (2.5e-4 *. scale) 2. + | _ -> 1. + in + Cairo.set_line_width ctx line_width; + Cairo.stroke ctx; + begin match cat with + `Pedestrian | `Track | `Path | `Cycleway | `Bridleway -> + Cairo.set_dash ctx [||] + | _ -> + () + end + in + draw_linear_features st ctx + (fun (cat, _, _, is_tunnel) -> + not is_tunnel || Linear_feature.of_id cat <> `Subway) + stroke i + in + +let t = Unix.gettimeofday () in + (* Draw surfaces *) + let fill_surface cat = + let cat = Surface.of_id cat in + set_surface_color ctx cat; + Cairo.fill ctx + in + let small_area = truncate ((*64.*)16. *. (10_000_000. /. scale) ** 2.) in + let filter (_, area, cat, ways) = + (area > small_area && + (area > 50_000_000 || Surface.of_id cat <> `Building)) + in + surfaces >> SP.with_group landuse >> SP.iter + >> draw_surfaces st ctx filter fill_surface; + (* Draw water lines (below buildings) *) + draw_water_lines (); + surfaces >> SP.with_group water >> SP.iter + >> draw_surfaces st ctx filter fill_surface; + surfaces >> SP.with_group building_or_pedestrian >> SP.iter + >> draw_surfaces st ctx filter fill_surface; + +if debug_time then +Format.eprintf "Drawing surfaces: %.3f@." (Unix.gettimeofday () -. t); + +let t = Unix.gettimeofday () in + (* Draw linear features *) + let l = + [minor_roads; major_road_links; major_roads; highway_links; highways] in + List.iter + (fun g -> + if g <> minor_roads then + linear_features >> LP.with_group g >> LP.iter >> draw_casing; + linear_features >> LP.with_group g >> LP.iter >> draw_inline) + l; +if debug_time then +Format.eprintf "Drawing lines: %.3f@." (Unix.gettimeofday () -. t) + +let partition_high = + ([ `Footway; `Steps; `Service; + `Tram; `Subway; `Taxiway; + `Pedestrian; `Track; `Cycleway; `Bridleway; `Path; + `Residential; `Unclassified; `Living_street; `Road; + `Stream ], + [ `River; `Canal ], + [ `Runway; `Rail; `Tertiary_link; `Tertiary ], + [ `Primary_link; `Primary; `Secondary_link; `Secondary ], + [ `Trunk_link; `Motorway_link; `Trunk; `Motorway ]) + +let partition_medium = + ([ `Footway; `Steps; `Service; + `Tram; `Subway; `Taxiway; `Runway; + `Pedestrian; `Track; `Cycleway; `Bridleway; `Path; + `Residential; `Unclassified; `Living_street; `Road; + `Tertiary_link; `Tertiary; + `Canal; `Stream ], + [ `River ], + [ `Rail; `Secondary_link; `Secondary ], + [ `Primary_link; `Primary ], + [ `Trunk_link; `Motorway_link; `Trunk; `Motorway ]) + +let partition_low = + ([ `Footway; `Steps; `Service; + `Tram; `Subway; `Taxiway; `Runway; + `Pedestrian; `Track; `Cycleway; `Bridleway; `Path; + `Residential; `Unclassified; `Living_street; `Road; + `Tertiary_link; `Tertiary; `Secondary_link; `Secondary; + `Canal; `Stream ], + [ `River ], + [ `Rail ], + [ `Primary_link; `Primary ], + [ `Trunk_link; `Motorway_link; `Trunk; `Motorway ]) + +let draw_map_intermediate_levels + st ctx road_partition surfaces linear_features = + let scale = compute_scale st in + let module SP = Surface.Partition in + let partition = SP.make () in + let landuse = + SP.add_group partition + [`Forest; `Grass; `Heath; `Rock; `Sand; `Glacier; + `Farmland; `Residential; `Commercial; + `Industrial; `Park; `Cemetery; `Parking ] + in + let water = SP.add_group partition [`Water] in + let building_or_pedestrian = + SP.add_group partition + [`Building; + `Highway_pedestrian; `Highway_track; `Highway_footway; `Highway_path ] + in + let highway_areas = + SP.add_group partition + [`Highway_residential; `Highway_unclassified; + `Highway_living_street; `Highway_service ] + in +let t = Unix.gettimeofday () in + let area = Array.map (fun (_, a, _, _) -> - log2 a) surfaces in + let layer = Array.map (fun (l, _, _, _) -> l) surfaces in +let surfaces' = surfaces in + let surfaces = + SP.apply partition surfaces (fun (_, _, cat, _) -> cat) + >> SP.order_totally + >> SP.order_by area + >> SP.order_by layer + >> SP.order_by_group + >> SP.select + in +if debug_time then +Format.eprintf "Sorting surfaces (%d elements): %.3f@." (Array.length surfaces') (Unix.gettimeofday () -. t); + +let t = Unix.gettimeofday () in + let module LP = Linear_feature.Partition in + let partition = LP.make () in + let (ignored_way_cats, waterway_cats, + minor_road_cats, major_road_cats, highway_cats) = + road_partition in + let waterway = LP.add_group partition waterway_cats in + let ignored_ways = LP.add_group partition ignored_way_cats in + let minor_roads = LP.add_group partition minor_road_cats in + let major_roads = LP.add_group partition major_road_cats in + let highways = LP.add_group partition highway_cats in + +let linear_features = + let eps = (10_000_000. /. scale) in + Array.map + (fun (((cat, _, _, _) as info, ways) as item) -> + if List.memq (Linear_feature.of_id cat) ignored_way_cats then + item + else + (info, List.filter (fun (x, y) -> Array.length x > 0) (List.map (fun (x, y) -> Osm_douglas_peucker.perform eps x y) ways))) + linear_features +in +(**) +let linear_features' = linear_features in + let linear_features = + LP.apply partition linear_features (fun ((cat, _, _, _), _) -> cat) + >> LP.order_totally + >> LP.order_by_group + >> LP.select + in +if debug_time then +Format.eprintf "Sorting lines (%d elements): %.3f@." (Array.length linear_features') (Unix.gettimeofday () -. t); +if debug_time then begin +let n = ref 0 in +Array.iter (fun (_, _, _, ways) -> List.iter (fun (x, _) -> n := !n + Array.length x) ways) surfaces'; +Format.eprintf "Surfaces: %d nodes@." !n; + let small_area = truncate (16. *. (*64. *.*) (10_000_000. /. scale) ** 2.) in + let filter (_, area, cat, ways) = + (area > small_area && + (area > 50_000_000 || Surface.of_id cat <> `Building)) + in +let n = ref 0 in +let m = ref 0 in +Array.iter (fun ((_, _, _, ways) as info) -> if filter info then begin incr m; List.iter (fun (x, _) -> n := !n + Array.length x) ways end) surfaces'; +Format.eprintf "Surfaces: %d nodes, %d elements (%.2f%%)@." !n !m (100. *. float !m /. float (Array.length surfaces')); +let n = ref 0 in +Array.iter (fun (_, ways) -> List.iter (fun (x, _) -> n := !n + Array.length x) ways) linear_features'; +Format.eprintf "Lines: %d nodes@." !n; +let n = ref 0 in +let n' = ref 0 in +let m = ref 0 in +Array.iter (fun ((cat, _, _, _), ways) -> if not (List.memq (Linear_feature.of_id cat) ignored_way_cats) then begin incr m; List.iter (fun (x, y) -> n := !n + Array.length x; let eps = (10_000_000. /. scale) in let (x', y') = Osm_douglas_peucker.perform eps x y in n' := !n' + Array.length x') ways end) linear_features'; +Format.eprintf "Lines: %d nodes (%d), %d elements (%.2f%%)@." !n !n' !m (100. *. float !m /. float (Array.length linear_features')); +end; + + let draw_water_lines () = + Cairo.set_line_cap ctx Cairo.BUTT; + Cairo.set_line_join ctx Cairo.JOIN_MITER; + Cairo.set_miter_limit ctx 1.4; + Cairo.set_source_rgb ctx 0.52 0.94 0.94; + Cairo.set_line_width ctx 2.; + linear_features >> LP.with_group waterway >> LP.iter + >> draw_linear_features st ctx + (fun (_, _, is_bridge, is_tunnel) -> not is_tunnel) + (fun _ -> Cairo.stroke ctx); + Cairo.set_dash ctx [||] + in + let emph_secondary = List.mem `Secondary major_road_cats in + let draw_casing i = + Cairo.set_line_cap ctx Cairo.BUTT; + Cairo.set_line_join ctx Cairo.JOIN_MITER; + Cairo.set_miter_limit ctx 1.4; + Cairo.set_dash ctx [||]; + let stroke (cat, _, _, is_tunnel) = + let cat = Linear_feature.of_id cat in + Cairo.set_source_rgb ctx 0. 0. 0.; + let line_width = + match cat with + `Trunk | `Motorway -> 8. + | `Trunk_link | `Motorway_link -> 5. + | `Primary -> 5. + | `Primary_link -> 4. + | `Secondary | `Secondary_link -> 3.5 + | _ -> assert false + in + Cairo.set_line_width ctx line_width; + Cairo.stroke ctx; + in + draw_linear_features st ctx + (fun (cat, _, is_bridge, _) -> + let cat = Linear_feature.of_id cat in + match cat with + `Motorway | `Trunk | `Primary | `Secondary + | `Motorway_link | `Trunk_link | `Primary_link | `Secondary_link -> + true + | `Secondary | `Secondary_link -> + emph_secondary + | _ -> + false) + stroke i + in + let draw_inline i = + Cairo.set_line_cap ctx Cairo.BUTT; + Cairo.set_line_join ctx Cairo.JOIN_MITER; + Cairo.set_miter_limit ctx 1.4; + let stroke (cat, _, _, is_tunnel) = + let cat = Linear_feature.of_id cat in + begin match cat with + `Motorway | `Trunk | `Motorway_link | `Trunk_link -> + Cairo.set_source_rgb ctx 1. 0.8 0. + | `Tertiary | `Tertiary_link -> + Cairo.set_source_rgb ctx 0.7 0.7 0.7 + | `Secondary | `Secondary_link when not emph_secondary -> + Cairo.set_source_rgb ctx 0.7 0.7 0.7 + | `Rail -> + Cairo.set_source_rgb ctx 0.27 0.27 0.27 + | `Runway -> + Cairo.set_source_rgb ctx 0.73 0.73 0.8 + | _ -> + Cairo.set_source_rgb ctx 1. 1. 1. + end; + let line_width = + match cat with + `Motorway | `Trunk -> 4. + | `Motorway_link | `Trunk_link -> 1. + | `Primary -> 3. + | `Primary_link -> 2. + | `Rail -> 1. + | `Secondary -> 1.5 + | `Secondary_link | `Tertiary_link | `Secondary | `Tertiary -> 1.5 + | `Runway -> max (2.5e-4 *. scale) 1. + | `Taxiway -> 1. + | _ -> assert false + in + Cairo.set_line_width ctx line_width; + Cairo.stroke ctx + in + draw_linear_features st ctx + (fun (cat, _, _, is_tunnel) -> true) + stroke i + in + +let t = Unix.gettimeofday () in + (* Draw surfaces *) + let fill_surface cat = + let cat = Surface.of_id cat in + set_surface_color ctx cat; + Cairo.fill ctx + in + let small_area = truncate ((*64.*)16. *. (10_000_000. /. scale) ** 2.) in + let filter (_, area, cat, ways) = + (area > small_area && + (area > 50_000_000 || Surface.of_id cat <> `Building)) + in + surfaces >> SP.with_group landuse >> SP.iter + >> draw_surfaces st ctx filter fill_surface; + (* Draw water lines (below buildings) *) + draw_water_lines (); + surfaces >> SP.with_group water >> SP.iter + >> draw_surfaces st ctx filter fill_surface; + surfaces >> SP.with_group building_or_pedestrian >> SP.iter + >> draw_surfaces st ctx filter fill_surface; + +if debug_time then +Format.eprintf "Drawing surfaces: %.3f@." (Unix.gettimeofday () -. t); + +let t = Unix.gettimeofday () in + (* Draw linear features *) + let l = [minor_roads; major_roads; highways] in + List.iter + (fun g -> + if g <> minor_roads then + linear_features >> LP.with_group g >> LP.iter >> draw_casing; + linear_features >> LP.with_group g >> LP.iter >> draw_inline) + l; +if debug_time then +Format.eprintf "Drawing lines: %.3f@." (Unix.gettimeofday () -. t) + +(****) + +let draw_map st pm x y width height = +(* +Format.eprintf "map: %d %d %d %d@." x y width height; +*) + let ctx = Cairo.create pm in + + (* Background *) + Cairo.rectangle ctx (float x) (float y) (float width) (float height); + Cairo.clip_preserve ctx; + Cairo.set_source_rgb ctx 1. 1. 1.; + Cairo.fill ctx; + + let x = st.rect.x + x in + let y = st.rect.y + y in + + Cairo.scale ctx 1. (-1.); + Cairo.translate ctx (-. float st.rect.x) (float st.rect.y); + + let extra = 7 in + + let scale = compute_scale st in + let lon_min = float (x - extra) /. scale in + let lon_max = float (x + width + extra) /. scale in + let lat_min = -. float (y + height + extra) /. scale in + let lat_max = -. float (y - extra) /. scale in + + (* Load surfaces *) +let t = Unix.gettimeofday () in + let surfaces = find_surfaces st.level lon_min lat_min lon_max lat_max in +if debug_time then +Format.eprintf "Loading surfaces: %.3f@." (Unix.gettimeofday () -. t); + + (* Load linear features *) +let t = Unix.gettimeofday () in + let linear_features = + find_linear_features st.level lon_min lat_min lon_max lat_max in +if debug_time then +Format.eprintf "Loading lines: %.3f@." (Unix.gettimeofday () -. t); + + (* Load coastline *) +let t = Unix.gettimeofday () in + let coastline = + find_coastline st.level lon_min lat_min lon_max lat_max in +if debug_time then +Format.eprintf "Loading coastline: %.3f@." (Unix.gettimeofday () -. t); + +if debug_time then +Format.eprintf "Surfaces: %d / lines: %d@." (2 * !surface_leaf_read) (2 * !linear_leaf_read); + + draw_coastline st ctx coastline lon_min lat_min lon_max lat_max; + if st.level >= 14.5 then + draw_map_high_levels st ctx surfaces linear_features + else if st.level >= 13.5 then + draw_map_medium_levels st ctx surfaces linear_features + else if st.level >= 12.5 then + draw_map_intermediate_levels + st ctx partition_high surfaces linear_features + else if st.level >= 11.5 then + draw_map_intermediate_levels + st ctx partition_medium surfaces linear_features + else + draw_map_intermediate_levels + st ctx partition_low surfaces linear_features + +(****) + +let draw_route st ctx = + let scale = compute_scale st in + Cairo.scale ctx 1. (-1.); + Cairo.translate ctx (-. float st.rect.x) (float st.rect.y); + let path = + List.map + (fun (lat, lon) -> + let approx x = + float ((x + linear_ratio / 2 - 1) / linear_ratio * linear_ratio) + in + (approx lon /. 10_000_000., + Osm_geometry.lat_to_y (approx lat) /. 10_000_000.)) + st.path + in + begin match path with + [] -> + () + | (x, y) :: rem -> + Cairo.move_to ctx (x *. scale) (y *. scale); + List.iter + (fun (x, y) -> Cairo.line_to ctx (x *. scale) (y *. scale)) + rem; + Cairo.set_line_width ctx 2.; + Cairo.set_source_rgb ctx 0. 0. 1.; + Cairo.stroke ctx + end; + + begin match st.marker1 with + Some (_, lat, lon) -> + let x = lon in + let y = Osm_geometry.lat_to_y (lat *. 10_000_000.) /. 10_000_000. in + Cairo.arc ctx (x *. scale) (y *. scale) 4. 0. (2. *. Osm_geometry.pi); + Cairo.set_source_rgb ctx 1. 0. 0.; + Cairo.fill ctx + | None -> + () + end; + begin match st.marker2 with + Some (_, lat, lon) -> + let x = lon in + let y = Osm_geometry.lat_to_y (lat *. 10_000_000.) /. 10_000_000. in + Cairo.arc ctx (x *. scale) (y *. scale) 4. 0. (2. *. Osm_geometry.pi); + Cairo.set_source_rgb ctx 0. 0.6 0.; + Cairo.fill ctx + | None -> + () + end diff --git a/osm/douglas_peucker.ml b/osm/lib/osm_douglas_peucker.ml similarity index 100% rename from osm/douglas_peucker.ml rename to osm/lib/osm_douglas_peucker.ml diff --git a/osm/douglas_peucker.mli b/osm/lib/osm_douglas_peucker.mli similarity index 100% rename from osm/douglas_peucker.mli rename to osm/lib/osm_douglas_peucker.mli diff --git a/osm/geometry.ml b/osm/lib/osm_geometry.ml similarity index 100% rename from osm/geometry.ml rename to osm/lib/osm_geometry.ml diff --git a/osm/geometry.mli b/osm/lib/osm_geometry.mli similarity index 100% rename from osm/geometry.mli rename to osm/lib/osm_geometry.mli diff --git a/osm/linear.ml b/osm/linear.ml index 5e6d4a2..3d70606 100644 --- a/osm/linear.ml +++ b/osm/linear.ml @@ -46,7 +46,7 @@ let compute_order out latitude longitude = if lat <> max_int then begin let lat = lat + 90_0000000 in let lon = lon + 180_0000000 in - Column.append o (Geometry.hilbert_coordinate lat lon); + Column.append o (Osm_geometry.hilbert_coordinate lat lon); loop () end in @@ -108,7 +108,7 @@ type state = mutable prev_bbox : Bbox.t; mutable bbox : Bbox.t } -module Feature = Category.Make (struct +module Feature = Osm_category.Make (struct type t = [ `Motorway | `Trunk | `Primary | `Secondary | `Tertiary | `Motorway_link | `Trunk_link | `Primary_link diff --git a/osm/multipolygons.ml b/osm/multipolygons.ml index 097d91f..7da8781 100644 --- a/osm/multipolygons.ml +++ b/osm/multipolygons.ml @@ -529,7 +529,7 @@ try if debug () then begin if let (_, (px, py)) = coords.(i) in - not (Geometry.is_simple_polygon px py) + not (Osm_geometry.is_simple_polygon px py) then begin Format.eprintf "NOT SIMPLE (%d)@." i; raise Exit @@ -538,7 +538,7 @@ try for j = 0 to len - 1 do if i <> j then inclusion.(i).(j) <- - Geometry.polygon_in_polygon (snd coords.(i)) (snd coords.(j)) + Osm_geometry.polygon_in_polygon (snd coords.(i)) (snd coords.(j)) done done; if debug () then @@ -552,11 +552,11 @@ try for j = 0 to i - 1 do if inclusion.(i).(j) && inclusion.(j).(i) then begin let mi = - Geometry.polygon_mostly_in_polygon + Osm_geometry.polygon_mostly_in_polygon (snd coords.(i)) (snd coords.(j)) in let mj = - Geometry.polygon_mostly_in_polygon + Osm_geometry.polygon_mostly_in_polygon (snd coords.(j)) (snd coords.(i)) in if debug () then diff --git a/osm/surfaces.ml b/osm/surfaces.ml index 3af3a38..a3678dd 100644 --- a/osm/surfaces.ml +++ b/osm/surfaces.ml @@ -107,7 +107,7 @@ let compute_order out latitude longitude = if lat <> max_int then begin let lat = lat + 90_0000000 in let lon = lon + 180_0000000 in - Column.append o (Geometry.hilbert_coordinate lat lon); + Column.append o (Osm_geometry.hilbert_coordinate lat lon); loop () end in @@ -118,7 +118,7 @@ let compute_order out latitude longitude = let _ = Column.set_database "/tmp/osm" -module Surface = Category.Make (struct +module Surface = Osm_category.Make (struct type t = [ `Water | `Forest | `Grass | `Heath | `Rock | `Sand | `Glacier | `Farmland | `Residential | `Commercial @@ -126,7 +126,7 @@ module Surface = Category.Make (struct | `Highway_residential | `Highway_unclassified | `Highway_living_street | `Highway_service | `Highway_pedestrian | `Highway_track | `Highway_footway | `Highway_path ] - let list = + let list = [ `Water; `Forest; `Grass; `Heath; `Rock; `Sand; `Glacier; `Farmland; `Residential; `Commercial; `Industrial; `Park; `Cemetery; `Parking; `Building; @@ -574,7 +574,7 @@ Printf.fprintf ch "%d\n" (lon.(i) - !last_lon); add_polygon let simplify_way ratio (lat, lon) = - let (lat, lon) = Douglas_peucker.perform_int ratio lat lon in + let (lat, lon) = Osm_douglas_peucker.perform_int ratio lat lon in let delta = ratio / 2 - 1 in let l = Array.length lat in for i = 0 to l - 1 do @@ -661,7 +661,7 @@ let _ = let (outer_way, inner_ways) = ways in let area = List.fold_left - (fun a (lat, lon) -> a + Geometry.polygon_area lon lat) 0 + (fun a (lat, lon) -> a + Osm_geometry.polygon_area lon lat) 0 (outer_way :: inner_ways) in if @@ -758,7 +758,7 @@ let _ = List.iter (fun (role, (lat, lon)) -> let sign = if role = 0 then 1 else -1 in - if Geometry.polygon_area lon lat * sign < 0 then begin + if Osm_geometry.polygon_area lon lat * sign < 0 then begin reverse lat; reverse lon end) ways;