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