SQL使用PostgreSQL
使用场景为生成马赛克风格的图片，即把某个区域分割成若干个小的区域，计算每个小的区域中包含的某些体量，用不同颜色显示。
1. 矩形（rectangle）
原理：每一行/列小矩形的中心点连线都是平行的，根据横纵切割的数量可以算出所有的中心点，再根据每个中心点计算出四个顶点的坐标。（如果想展现正方形可以1.手动控制横纵切割的比;2.简单修改代码逻辑，只传入横向切割数量，用边长计算纵向数量。）
def cal_rasterized_points_from_boundary_v2(obj, boundary_li, lng_count, lat_count):
"""
计算边界内栅格化后的中心点，返回若干个矩形的坐标数组 -- 马赛克凹凸有致版
:param obj: 多边形的object
:param boundary_li: multipolygon.boundary.coords转化后的list
:param lng_count: 横向切割数量
:param lat_count: 纵向切割数量
:return:
"""
# 目前只考虑MultiPolygon包含一个多边形的情况。
boundary_list = boundary_li[0]
lngs = []
lats = []
for bl in boundary_list:
lngs.append(bl[0])
lats.append(bl[1])
max_lng = float(max(lngs))
min_lng = float(min(lngs))
max_lat = float(max(lats))
min_lat = float(min(lats))
max_lng_amap, max_lat_amap = wgs84togcj02(max_lng, max_lat)
min_lng_amap, min_lat_amap = wgs84togcj02(min_lng, min_lat)
lng_gap = (max_lng_amap - min_lng_amap) / lng_count
lat_gap = (max_lat_amap - min_lat_amap) / lat_count
data_lngs = []
data_lats = []
cps = []
result = []
for i in range(lng_count):
lng = min_lng_amap + lng_gap * 0.5 + lng_gap * i
if lng < max_lng_amap:
data_lngs.append(lng)
for i in range(lat_count):
lat = min_lat_amap + lat_gap * 0.5 + lat_gap * i
if lat < max_lat_amap:
data_lats.append(lat)
for lat in data_lats:
for lng in data_lngs:
cps.append([lng, lat])
for cp in cps:
# TODO: To judge cp whether in boundary or not
lng, lat = gcj02towgs84(cp[0], cp[1])
point = Point((lng, lat))
li = obj.filter(boundary__contains=point)
if not li:
continue
result.append([[cp[0] - lng_gap * 0.5, cp[1] - lat_gap * 0.5],
[cp[0] + lng_gap * 0.5, cp[1] - lat_gap * 0.5],
[cp[0] + lng_gap * 0.5, cp[1] + lat_gap * 0.5],
[cp[0] - lng_gap * 0.5, cp[1] + lat_gap * 0.5]])
return result
2. 六边形（hexagonal）
原理：
20170124_01_pic_004.png
六边形可以拆分成六个等边三角形。
20170124_01_pic_010.png
每个中心点的位置每两行一个循环。同行两个中心点间距为三倍的边长，行间距为√3/2倍边长（经纬度间距离转换存在差异，现采用纬度间*0.88）
def cal_rasterized_points_from_boundary(obj, boundary_li, lng_count):
"""
计算边界内栅格化后的中心点，返回若干个六边形的坐标数组 -- 马赛克凹凸有致版
蜂巢六边形切割 -- 粗略用lng的差来表示边长，因为分割后数值很小，可以忽略坐标和距离转换的误差
:param obj: 多边形的object
:param boundary_li: multipolygon.boundary.coords转化后的list
:param lng_count: 横向切割数量
:return:
"""
# 目前只考虑MultiPolygon包含一个多边形的情况。
boundary_list = boundary_li[0]
lngs = []
lats = []
for bl in boundary_list:
lngs.append(bl[0])
lats.append(bl[1])
max_lng = float(max(lngs))
min_lng = float(min(lngs))
max_lat = float(max(lats))
min_lat = float(min(lats))
max_lng_amap, max_lat_amap = wgs84togcj02(max_lng, max_lat)
min_lng_amap, min_lat_amap = wgs84togcj02(min_lng, min_lat)
lng_gap = (max_lng_amap - min_lng_amap) / lng_count
side_length = lng_gap / 3
lat_gap = side_length * math.sqrt(3) / 2 * 0.88
lat_count = int((max_lat_amap - min_lat_amap) / lat_gap) + 1
data_lngs_first = []
data_lngs_second = []
data_lats = []
cps = []
result = []
for i in range(lng_count):
lng = min_lng_amap + lng_gap * 0.5 + lng_gap * i
if lng < max_lng_amap:
data_lngs_first.append(lng)
for lng in data_lngs_first:
sec_lng = lng - side_length * 3 / 2
if sec_lng < max_lng_amap:
data_lngs_second.append(sec_lng)
for i in range(lat_count):
lat = min_lat_amap + lat_gap * 0.5 + lat_gap * i
if lat < max_lat_amap:
data_lats.append(lat)
for i, lat in enumerate(data_lats):
# 六边形单双行中心点不对齐
if i % 2 == 0:
for lng in data_lngs_first:
cps.append([lng, lat])
else:
for lng in data_lngs_second:
cps.append([lng, lat])
for cp in cps:
lng, lat = gcj02towgs84(cp[0], cp[1])
point = Point((lng, lat))
li = obj.filter(boundary__contains=point)
if not li:
continue
# 左上角为第一个点
result.append([[cp[0] - side_length * 0.5, cp[1] + lat_gap],
[cp[0] + side_length * 0.5, cp[1] + lat_gap],
[cp[0] + side_length, cp[1]],
[cp[0] + side_length * 0.5, cp[1] - lat_gap],
[cp[0] - side_length * 0.5, cp[1] - lat_gap],
[cp[0] - side_length, cp[1]]])
return result
def generate_comb_mosaic_item(bc_id, lng_count):
area = BusinessCircle.objects.filter(business_circle_id=bc_id)
if not area:
return ""
mp = area[0].boundary
if len(mp) > 1:
return "返回的是多个多边形，这样不ok"
boundary_li = [[list(list(location)) for location in list(polygon)] for polygon in list(mp.boundary.coords)]
result = cal_rasterized_points_from_boundary(area, boundary_li, lng_count)
# data = [[j, i] for i in range(lat_count) for j in range(lng_count)]
data = []
cursor = connection.cursor()
SQL_string_list = []
dp_result = copy.deepcopy(result)
for index, points in enumerate(dp_result):
MULTIPOLYGON_txt = ""
points.append(points[0])
for point in points:
lng, lat = gcj02towgs84(point[0], point[1])
MULTIPOLYGON_txt += ", {} {}".format(lng, lat)
MULTIPOLYGON_txt = MULTIPOLYGON_txt[2:]
SQL_RAW = """SELECT avg(forecast_price) FROM economic_dqchina_loupaninfo
WHERE st_contains(st_geomfromtext('MULTIPOLYGON((({})))', 4326) :: geometry, center_point :: geometry) UNION ALL """.format(MULTIPOLYGON_txt)
SQL_string_list.append(SQL_RAW)
SQL_RAWS = "".join(SQL_string_list)
cursor.execute(SQL_RAWS[:-10])
rows = cursor.fetchall()
for index, row in enumerate(rows):
data.append(row[0])
item = {
"data": [float(d) if d else 0 for d in data],
"result": result,
}
# print(item["data"])
# print(item["result"])
return item
