shp.c

   1 /*
   2  *      Hic Est Leo -- Reading ESRI Shape Files
   3  *
   4  *      (c) 2014--2015 Martin Mares <mj@ucw.cz>
   5  *
   6  *      FIXME: Currently, this parser handles only the subset
   7  *      of shape file syntax which is used by gdal_contours.
   8  */
   9
  10 #undef LOCAL_DEBUG
  11
  12 #include "leo.h"
  13 #include "osm.h"
  14 #include "map.h"
  15 #include "shp.h"
  16
  17 #include <ucw/conf.h>
  18 #include <ucw/fastbuf.h>
  19 #include <ucw/gary.h>
  20 #include <ucw/unaligned.h>
  21
  22 #include <fcntl.h>
  23 #include <stdio.h>
  24
  25 static osm_id_t shp_id_counter;
  26
  27 static double shp_get_double(byte *p)
  28 {
  29   // FIXME: This is non-portable!
  30   double x = 0;
  31   memcpy(&x, p, 8);
  32   return x;
  33 }
  34
  35 static void shp_record_polyline(byte *buf, u32 len)
  36 {
  37   if (len < 44)
  38     die("Polyline record too short: %u bytes", len);
  39
  40   u32 num_parts = get_u32_le(buf+36);
  41   u32 num_points = get_u32_le(buf+40);
  42   DBG("%u parts, %u points", num_parts, num_points);
  43   if (num_parts != 1)
  44     die("Polylines with multiple parts are not supported yet");
  45
  46   if (len < 44 + 4*num_parts + 16*num_points)
  47     die("Polyline record too short for %u parts and %u points: %u bytes", num_parts, num_points, len);
  48
  49   struct osm_way *w = osm_way_new(shp_id_counter++);
  50   osm_obj_add_tag(&w->o, "shape", "contour");
  51
  52   byte *p = buf + 44 + 4*num_parts;
  53   for (uint i=0; i < num_points; i++)
  54     {
  55       double x = shp_get_double(p);
  56       double y = shp_get_double(p+8);
  57       DBG(">>> x=%.10g y=%.10g", x, y);
  58       struct osm_node *n = osm_node_new(shp_id_counter++);
  59       n->x = x;
  60       n->y = y;
  61       osm_way_add_node(w, n);
  62       p += 16;
  63     }
  64 }
  65
  66 void shp_parse(const char *name)
  67 {
  68   msg(L_INFO, "Loading shape file %s", name);
  69   struct fastbuf *fb = bopen_file(name, O_RDONLY, NULL);
  70   shp_id_counter = 1;
  71
  72   byte *buf;
  73   GARY_INIT(buf, 100);
  74
  75   uns len = bread(fb, buf, 100);
  76   if (len != 100 || get_u32_be(buf) != 9994)
  77     die("Invalid shape file header");
  78
  79   u32 version = get_u32_le(buf+28);
  80   if (version != 1000)
  81     die("Unknown shape file version %u", version);
  82
  83   u32 type = get_u32_le(buf+32);
  84   if (type != 3)
  85     die("Unsupported shape file type %u", type);
  86
  87   u32 last_recno = 0;
  88   for (;;)
  89     {
  90       u64 pos = btell(fb);
  91       uns len = bread(fb, buf, 8);
  92       if (!len)
  93         break;
  94       if (len != 8)
  95         die("Truncated record header at %ju", (uintmax_t) pos);
  96       u32 recno = get_u32_be(buf);
  97       u32 reclen = get_u32_be(buf+4) * 2;
  98       if (recno != ++last_recno)
  99         die("Unexpected record #%u (should be #%u) at %ju", recno, last_recno, (uintmax_t) pos);
 100       if (reclen > GARY_SIZE(buf))
 101         GARY_RESIZE(buf, reclen);
 102       len = bread(fb, buf, reclen);
 103       if (len != reclen)
 104         die("Truncated record: %u < %u at %ju", len, reclen, (uintmax_t) pos);
 105       if (len < 4)
 106         die("Record too short: %u bytes at %ju", reclen, (uintmax_t) pos);
 107
 108       u32 rectype = get_u32_le(buf);
 109       DBG("@%ju: recno=%u len=%u type=%u", (uintmax_t) pos, recno, reclen, rectype);
 110       if (rectype == 3)
 111         shp_record_polyline(buf, len);
 112     }
 113
 114   GARY_FREE(buf);
 115   bclose(fb);
 116 }