clean up grid rounding error

This commit is contained in:
admin 2026-02-02 06:20:32 -08:00
parent 9adc57e7a0
commit 8904c7d966

View File

@ -31,6 +31,8 @@ class Topology:
self.iso_country = iso_country
self.subregion = subregion
self.linemap_size = linemap_size
self.linemap_xscale = 1.0 * self.linemap_size[0] / linemap_cols
self.linemap_yscale = 1.0 * self.linemap_size[1] / linemap_rows
self.linemap_rows = linemap_rows
self.linemap_cols = linemap_cols
self.linemap_width = linemap_width
@ -55,6 +57,7 @@ class Topology:
self.target_region_name=self.subregion
self.elevation_filename = self.storage_directory + "/" + self.target_region_name+".png"
print("Elevation file name is [{}]".format(self.elevation_filename))
print("requested linemap size is {}, xscale is [{}], yscale is [{}]".format(self.linemap_size, self.linemap_xscale, self.linemap_yscale))
def list_subunits(self):
@ -118,6 +121,7 @@ class Topology:
frame_image.paste(elevation_data,
((1224 - elevation_data.size[0]) // 2,
(1224 - elevation_data.size[1]) // 2))
frame_image.save(self.elevation_filename+".frame.png")
if rotation is None:
elevation_data = frame_image
@ -125,21 +129,33 @@ class Topology:
elevation_data = frame_image.rotate(rotation, expand=False)
width, height = elevation_data.size
print(width)
print(height)
grid_w, grid_h = self.linemap_rows, self.linemap_cols
# Calculate cell dimensions
cell_w = width // grid_w
cell_h = height // grid_h
cell_w = 1.0 * width / grid_w
cell_h = 1.0 * height / grid_h
print("Cell width is w,h [{},{}]".format(cell_w, cell_h))
# Calculate mean values for each cell
# future improve -- only average non-zero values to get coastlines and edges right
# need to use xscale yscale to get cell window, adjust for rounding/floor error
values = np.zeros((grid_h, grid_w))
for r in range(grid_h):
for c in range(grid_w):
# Define box (left, upper, right, lower)
box = (c * cell_w, r * cell_h, (c + 1) * cell_w, (r + 1) * cell_h)
box = (int(c * cell_w),
int(r * cell_h),
int((c + 1) * cell_w),
int((r + 1) * cell_h)
)
cell = elevation_data.crop(box)
values[r, c] = np.mean(np.array(cell))
this_array = np.array(cell)
if np.max(this_array) == 0:
values[r,c] = 0
else:
values[r, c] = np.mean(this_array[this_array != 0])
return values
def create_svg_file(self, rotation=None,file_index=None):
@ -177,33 +193,34 @@ class Topology:
cr = cairo.Context(surface)
# Background
cr.set_source_rgb(1, 1, 1)
#cr.set_source_rgb(1, 1, 1)
cr.set_source_rgb(0, 0, 0)
cr.paint()
# Simple Orthographic Projection Matrix (pseudo-3D)
# Projects (x,y,z) -> (x+z/2, y+z/2)
# max height should be no more than 3 times the distance between rows, this could be a setting
def project(x, y, z, xscale=20, yscale=10, offset=(0, 0)):
shelf = 0 if z == 0 else z + 2
shelf = z
return (offset[0] + x * xscale, offset[1] + y * yscale + shelf)
def project(x, y, z):
return (x * self.linemap_xscale, y * self.linemap_yscale - z)
# return (offset[0] + x * scale + z * 2, offset[1] + y * scale + z * 2)
rows, cols = values.shape
max_val = values.max()
# Normalize data for visualization height (0-100)
norm_data = (values / max_val) * 100
# Normalize data for visualization height (0-100) -- max is no more than three row heights
norm_data = (values / max_val) * self.linemap_xscale * 3
# Draw Lines along rows
cr.set_source_rgb(0, 0, 0) # Black lines
cr.set_line_width(1.0)
cr.set_line_width(2.0)
cr.set_line_cap(cairo.LINE_CAP_ROUND)
cr.set_line_join(cairo.LINE_JOIN_ROUND)
x = 0
y = 0
for r in range(rows):
for c in range(cols):
z = norm_data[r, c]
x, y = project(c, r, -z, xscale=800 / cols, yscale=800 / rows) # -z to make higher values go "up" in 2D
x, y = project(c, r, z) # -z to make higher values go "up" in 2D
if c == 0:
cr.move_to(x, y + 50)
cr.line_to(x, y)
@ -211,27 +228,29 @@ class Topology:
cr.line_to(x, y)
cr.line_to(x, y + 50)
cr.close_path()
cr.set_source_rgb(1, 1, 1)
cr.fill_preserve()
cr.set_source_rgb(0, 0, 0)
cr.set_line_width(1.0)
cr.fill_preserve()
cr.set_source_rgb(1, 1, 1)
cr.set_line_width(2.0)
cr.stroke()
surface.write_to_png(svg_filename.replace(".svg", ".png"))
surface.write_to_png(svg_filename.replace(".svg", ".topo.png"))
surface.finish()
def main():
topology = Topology(country="Spain",
topology = Topology(country="Switzerland",
storage_directory="/Users/he/PycharmProjects/toposhirt/outputfiles",
subregion_mode=False,
linemap_cols=170,
linemap_rows=270)
linemap_cols=150,
linemap_rows=150,
linemap_size=(1224,1224))
topology.list_subunits()
print(topology.get_boundingbox())
topology.create_svg_file()
sys.exit(0)
for r in range(4*360):
print("Rotating by [{}] degrees".format(r/4.0))
topology.create_svg_file(rotation=(r / 4.0), file_index=r)
sys.exit(0)
t2=Topology(country="United Kingdom",
subregion_mode=True,