2021-08-23 20:03:53 (UTC-03:00)
Marcel Rodrigues <marcelgmr@gmail.com>
map: implement scaled poly cache
diff --git a/map.lua b/map.lua index 1c932cb..c801955 100644 --- a/map.lua +++ b/map.lua @@ -1,3 +1,5 @@ +local bio = require "bio" + local abs = math.abs local rad, deg = math.rad, math.deg local cos, sin, tan = math.cos, math.sin, math.tan @@ -197,6 +199,72 @@ function Frame:mapped(polys) end end +function Frame:near(bb) + -- not strictly correct, since projected geobbox != 2D rectangle + -- just used to limit cache and avoid 16-bit overflows + local ns = {} + local add_ll = function(lon, lat) + local x, y = self:map(lon, lat) + table.insert(ns, abs(x-self.w/2)) + table.insert(ns, abs(y-self.h/2)) + end + add_ll(bb.xmin, bb.ymin) + add_ll(bb.xmax, bb.ymin) + add_ll(bb.xmax, bb.ymax) + add_ll(bb.xmin, bb.ymax) + return math.max(unpack(ns)) < math.max(self.w, self.h)*2 +end + +function Frame:save_cache(fname, sf, filter) + local indices = {} + for i = 1, #sf.tab do + local row = sf.tab[i] + local rec = {} + for j = 1, #sf.fields do + rec[sf.fields[j].name] = row[j] + end + local key = filter(rec) + if key ~= nil then + if self:near(sf:read_bbox(i)) then + table.insert(indices, {i, key:sub(1, 16)}) + end + end + end + local cache = io.open(fname, "w") + for i = 1, #indices do + local index, key = unpack(indices[i]) + cache:write(key, string.rep("\0", 20 - #key)) + end + cache:write(string.rep("\0", 20)) -- end list of entries + for i = 1, #indices do + local index, key = unpack(indices[i]) + local offset = cache:seek() + cache:seek("set", i * 20 - 4) + bio.write_beu32(cache, offset) + cache:seek("set", offset) + local bb, lens, polys = sf:read_polygons(index) + for poly in self:mapped(polys) do + local ox, oy = unpack(poly()) + ox, oy = math.floor(ox + 0.5), math.floor(oy + 0.5) + bio.write_bei16(cache, ox) + bio.write_bei16(cache, oy) + for point in poly do + local x, y = unpack(point) + x, y = math.floor(x + 0.5), math.floor(y + 0.5) + local dx, dy = x-ox, y-oy + if dx ~= 0 or dy ~= 0 then + bio.write_bei16(cache, dx) + bio.write_bei16(cache, dy) + ox, oy = x, y + end + end + bio.write_beu32(cache, 0) -- end list of points + end + bio.write_beu32(cache, 0) -- end list of polys + end + cache:close() +end + function Frame:save(fname) local frm = io.open(fname, "w") local model = self.proj:model() @@ -259,7 +327,60 @@ local function load_frame(fname) return new_frame(proj, bbox) end +local Cache = {} +Cache.__index = Cache + +function Cache:keys() + local cache = self.fp + local offset = 0 + return function() + cache:seek("set", offset) + offset = offset + 20 + local ckey = cache:read(16) + if ckey:byte() ~= 0 then + return ckey:sub(1, ckey:find("\0")-1) + end + end +end + +function Cache:get_polys(key) + local cache = self.fp + cache:seek("set", 0) + local ckey = cache:read(16) + local offset = -1 + while ckey:byte() ~= 0 do + if ckey:sub(1, ckey:find("\0")-1) == key then + offset = bio.read_beu32(cache) + break + end + cache:seek("cur", 4) + ckey = cache:read(16) + end + assert(offset > 0, ("key '%s' not found in cache '%s'"):format(key, fname)) + cache:seek("set", offset) + return function() + local ox, oy = bio.read_bei16(cache), bio.read_bei16(cache) + if ox ~= 0 or oy ~= 0 then + local x, y = 0, 0 + return function() + if x ~= ox or y ~= oy then + x, y = ox, oy + local dx, dy = bio.read_bei16(cache), bio.read_bei16(cache) + ox, oy = ox+dx, oy+dy + return {x, y} + end + end + end + end +end + +local function load_cache(fname) + local self = setmetatable({}, Cache) + self.fp = io.open(fname, "r") + return self +end + return { distance=distance, bbox=bbox, centroid=centroid, Proj=Proj, - new_frame=new_frame, load_frame=load_frame + new_frame=new_frame, load_frame=load_frame, load_cache=load_cache } diff --git a/shp.lua b/shp.lua index d102a5b..d6b3f10 100644 --- a/shp.lua +++ b/shp.lua @@ -1,6 +1,7 @@ -local ffi = require "ffi" local bit = require "bit" +local bio = require "bio" + local function rtrim(s, c) local i = #s while s:sub(i, i) == c do @@ -9,76 +10,30 @@ local function rtrim(s, c) return s:sub(1, i) end -local function read_byte(fp) - return fp:read(1):byte(1) -end - -local function read_leu16(fp) - return read_byte(fp) + bit.lshift(read_byte(fp), 8) -end - -local function read_leu32(fp) - return read_leu16(fp) + bit.lshift(read_leu16(fp), 16) -end - -local function read_beu16(fp) - return bit.lshift(read_byte(fp), 8) + read_byte(fp) -end - -local function read_beu32(fp) - return bit.lshift(read_beu16(fp), 16) + read_beu16(fp) -end - -local function signed(u, nbits) - local p = 2^(nbits-1) - if u >= p then u = u - p - p end - return u -end - -local function read_lei16(fp) - return signed(read_leu16(fp), 16) -end - -local function read_lei32(fp) - return signed(read_leu32(fp), 32) -end - -local function read_bei16(fp) - return signed(read_beu16(fp), 16) -end - -local function read_bei32(fp) - return signed(read_beu32(fp), 32) -end - -local function read_led64(fp) - return ffi.cast("double *", ffi.new("char[8]", fp:read(8)))[0] -end - local SF = {} SF.__index = SF function SF:read_dbf() local fp = io.open(self.path .. ".dbf", "rb") - local version = bit.band(read_byte(fp), 0x07) + local version = bit.band(bio.read_byte(fp), 0x07) if version ~= 3 then error("only DBF version 5 is supported") end - self.year = 1900 + read_byte(fp) - self.month = read_byte(fp) - self.day = read_byte(fp) - self.nrecs = read_leu32(fp) + self.year = 1900 + bio.read_byte(fp) + self.month = bio.read_byte(fp) + self.day = bio.read_byte(fp) + self.nrecs = bio.read_leu32(fp) fp:seek("cur", 24) local reclen = 0 local fields = {} - local byte = read_byte(fp) + local byte = bio.read_byte(fp) while byte ~= 0x0D do local field_name = rtrim(string.char(byte) .. fp:read(10), "\000") local field_type = fp:read(1) fp:seek("cur", 4) - local field_length = read_byte(fp) + local field_length = bio.read_byte(fp) reclen = reclen + field_length - local field_dec_count = read_byte(fp) + local field_dec_count = bio.read_byte(fp) local field = { name=field_name, type=field_type, @@ -87,7 +42,7 @@ function SF:read_dbf() } table.insert(fields, field) fp:seek("cur", 14) - byte = read_byte(fp) + byte = bio.read_byte(fp) end local tab = {} for i = 1, self.nrecs do @@ -164,32 +119,32 @@ end local function read_bbox(fp) local bbox = {} - bbox.xmin = read_led64(fp) - bbox.ymin = read_led64(fp) - bbox.xmax = read_led64(fp) - bbox.ymax = read_led64(fp) + bbox.xmin = bio.read_led64(fp) + bbox.ymin = bio.read_led64(fp) + bbox.xmax = bio.read_led64(fp) + bbox.ymax = bio.read_led64(fp) return bbox end -- SHX and SHP share the same header structure! local function read_sf_header(fp) - local file_code = read_bei32(fp) + local file_code = bio.read_bei32(fp) if file_code ~= 9994 then error("does not seem to be a shapefile (invalid FileCode)") end fp:seek("cur", 20) -- unused, always all-zero - local file_len = read_bei32(fp) - local version = read_lei32(fp) + local file_len = bio.read_bei32(fp) + local version = bio.read_lei32(fp) if version ~= 1000 then error("invalid shapefile version") end local header = {} - header.shape = shape_name[read_lei32(fp)] + header.shape = shape_name[bio.read_lei32(fp)] header.bbox = read_bbox(fp) - header.zmin = read_led64(fp) - header.zmax = read_led64(fp) - header.mmin = read_led64(fp) - header.mmax = read_led64(fp) + header.zmin = bio.read_led64(fp) + header.zmax = bio.read_led64(fp) + header.mmin = bio.read_led64(fp) + header.mmax = bio.read_led64(fp) return header end @@ -198,8 +153,8 @@ function SF:read_shx() self.header = read_sf_header(fp) self.index = {} for i = 1, self.nrecs do - local offset = read_bei32(fp) - local length = read_bei32(fp) + local offset = bio.read_bei32(fp) + local length = bio.read_bei32(fp) table.insert(self.index, {offset=offset, length=length}) end io.close(fp) @@ -208,12 +163,27 @@ end function SF:read_record_header(index) local fp = self.fp fp:seek("set", self.index[index].offset * 2) - local num = read_bei32(fp) - local len = read_bei32(fp) - local shape = shape_name[read_lei32(fp)] + local num = bio.read_bei32(fp) + local len = bio.read_bei32(fp) + local shape = shape_name[bio.read_lei32(fp)] return num, len, shape end +function SF:read_bbox(index) + local shape_name = self.header.shape + assert( + startswith(shape_name, "polyline") or + startswith(shape_name, "polygon") or + startswith(shape_name, "multipoint"), + ("%s shape doesn't have bbox"):format(shape_name) + ) + local num, len, shape = self:read_record_header(index) + if shape == "null" then + return nil + end + return read_bbox(self.fp) +end + function SF:read_polygons(index) assert(startswith(self.header.shape, "polygon"), "type mismatch") local fp = self.fp @@ -222,13 +192,13 @@ function SF:read_polygons(index) return nil end local bbox = read_bbox(fp) - local nparts = read_lei32(fp) + local nparts = bio.read_lei32(fp) --~ print(nparts) - local npoints = read_lei32(fp) - local total = read_lei32(fp) -- first index is always zero? + local npoints = bio.read_lei32(fp) + local total = bio.read_lei32(fp) -- first index is always zero? local lengths = {} for i = 2, nparts do - local index = read_lei32(fp) + local index = bio.read_lei32(fp) table.insert(lengths, index - total) total = index end @@ -242,8 +212,8 @@ function SF:read_polygons(index) return function() j = j + 1 if j <= length then - local x = read_led64(fp) - local y = read_led64(fp) + local x = bio.read_led64(fp) + local y = bio.read_led64(fp) return {x, y} end end