从给定的 [lat, lon, density] 数据集生成 KML 热图

Question

我希望构建一个静态 KML（Google 地球标记）文件，该文件以 [lat, lon, density] 元组的形式显示一些给定数据集的热图样式渲染。

我拥有的一个非常简单的数据集是人口密度。

我的要求是：

• 必须能够输入给定纬度、经度的数据
• 必须能够指定该纬度、经度的数据密度
• 必须导出到 KML

此项目的要求与语言无关，因为我将离线生成这些文件，以便构建在其他地方使用的 KML。

我查看了一些项目，最值得注意的是heatmap.py，它是带有 KML 导出功能的 Python 中 gheat 的一个端口。从某种意义上说，我发现迄今为止的项目都依赖于从输入算法的 [lat, lon] 点的密度构建热图，因此我遇到了障碍。

如果我缺少一种明显的方法来调整我的数据集以仅输入 [lat, lon] 元组，但使用我拥有的密度值调整我输入它们的方式，我很想知道！

Answer 1

嘿，威尔，heatmap.py 是我。您的请求很常见，并且在我的待处理事项清单上。我还不太确定如何以一般方式这样做。在 heatmap.py 的说法中，使用每点 dotsize 而不是像现在这样的全局点大小会很简单，但我不确定这是否能满足真正的需求。我的目标是在 2010 年夏季发布，但你可以自己制作这个模组。

您可以尝试搜索Kernel Density Estimator工具；这就是统计学家所说的热图。 R 有一些不错的内置工具可供您使用，它们可能会更快地满足您的需求。

祝你好运！

Answer 2

我更新了heatmap.py 脚本，以便您可以为每个点指定密度。我uploaded my changes to my blog。不确定它是否会完全符合您的要求！

干杯， 亚历克斯

更新 [2020 年 11 月 13 日] 不久前我存档了我的博客，因此链接不再有效，因此以下是更改以供参考：

差异文件：

--- __init__.py 2010-09-14 08:40:35.829079482 +0100
+++ __init__.py.mynew   2010-09-06 14:50:10.394447647 +0100
@@ -1,5 +1,5 @@
 #heatmap.py v1.0 20091004
-from PIL import Image,ImageChops
+from PIL import Image,ImageChops,ImageDraw
 import os
 import random
 import math
@@ -43,10 +43,13 @@
     Most of the magic starts in heatmap(), see below for description of that function.
     """
     def __init__(self):
+        self.minIntensity = 0
+        self.maxIntensity = 0
         self.minXY = ()
         self.maxXY = ()
+       
 
-    def heatmap(self, points, fout, dotsize=150, opacity=128, size=(1024,1024), scheme="classic"):
+    def heatmap(self, points, fout, dotsize=150, opacity=128, size=(4048,1024), scheme="classic", area=(-180,180,-90,90)):
         """
         points  -> an iterable list of tuples, where the contents are the 
                    x,y coordinates to plot. e.g., [(1, 1), (2, 2), (3, 3)]
@@ -59,33 +62,41 @@
         size    -> tuple with the width, height in pixels of the output PNG 
         scheme  -> Name of color scheme to use to color the output image.
                    Use schemes() to get list.  (images are in source distro)
+        area    -> specify the coordinates covered by the resulting image 
+                   (could create an image to cover area larger than the max/
+                   min values given in the points list) 
         """
-        
+        print("Starting heatmap")
         self.dotsize = dotsize
         self.opacity = opacity
         self.size = size
         self.imageFile = fout
- 
+
         if scheme not in self.schemes():
             tmp = "Unknown color scheme: %s.  Available schemes: %s"  % (scheme, self.schemes())                           
             raise Exception(tmp)
 
-        self.minXY, self.maxXY = self._ranges(points)
-        dot = self._buildDot(self.dotsize)
+        self.minXY = (area[0],area[2])
+        self.maxXY = (area[1],area[3])
 
+        self.minIntensity, self.maxIntensity = self._intensityRange(points)
+        
         img = Image.new('RGBA', self.size, 'white')
-        for x,y in points:
+        for x,y,z in points:
+            dot = self._buildDot(self.dotsize,z)
             tmp = Image.new('RGBA', self.size, 'white')
             tmp.paste( dot, self._translate([x,y]) )
             img = ImageChops.multiply(img, tmp)
 
-
+        print("All dots built")
         colors = colorschemes.schemes[scheme]
         img.save("bw.png", "PNG")
+        print("Saved temp b/w image")
+        print("Colourising")
         self._colorize(img, colors)
 
         img.save(fout, "PNG")
-
+        print("Completed colourising and saved final image %s" % fout)
     def saveKML(self, kmlFile):
         """ 
         Saves a KML template to use with google earth.  Assumes x/y coordinates 
@@ -110,17 +121,19 @@
         """
         return colorschemes.schemes.keys() 
 
-    def _buildDot(self, size):
+    def _buildDot(self, size,intensity):
         """ builds a temporary image that is plotted for 
             each point in the dataset"""
+        
+        intsty = self._calcIntensity(intensity)
+        print("building dot... %d: %f" % (intensity,intsty))
+
         img = Image.new("RGB", (size,size), 'white')
-        md = 0.5*math.sqrt( (size/2.0)**2 + (size/2.0)**2 )
-        for x in range(size):
-            for y in range(size):
-                d = math.sqrt( (x - size/2.0)**2 + (y - size/2.0)**2 )
-                rgbVal = int(200*d/md + 50)
-                rgb = (rgbVal, rgbVal, rgbVal)
-                img.putpixel((x,y), rgb)
+        draw  =  ImageDraw.Draw(img)
+        shade = 256/(size/2)
+        for x in range (int(size/2)):
+            colour = int(256-(x*shade*intsty))
+            draw.ellipse((x,x,size-x,size-x),(colour,colour,colour))
         return img
 
     def _colorize(self, img, colors):
@@ -139,7 +152,7 @@
                 rgba.append(alpha) 
 
                 img.putpixel((x,y), tuple(rgba))
-            
+     
     def _ranges(self, points):
         """ walks the list of points and finds the 
         max/min x & y values in the set """
@@ -153,6 +166,23 @@
             
         return ((minX, minY), (maxX, maxY))
 
+    def _calcIntensity(self,z):
+        return (z/self.maxIntensity)        
+               
+    def _intensityRange(self, points):
+        """ walks the list of points and finds the 
+        max/min points of intensity
+        """   
+        minZ = points[0][2]
+        maxZ = minZ
+        
+        for x,y,z in points:
+            minZ = min(z, minZ)
+            maxZ = max(z, maxZ)
+        
+        print("(minZ, maxZ):(%d, %d)" % (minZ,maxZ))
+        return (minZ, maxZ)
+        
     def _translate(self, point):
         """ translates x,y coordinates from data set into 
         pixel offsets."""

和一个演示脚本：

import heatmap
import random
import MySQLdb
import math

print "starting script..."

db = MySQLdb.connect(host="localhost", # your host, usually localhost
                     user="username", # your username
                      passwd="password", # your password
                      db="database") # name of the data base
cur = db.cursor() 

minLng = -180
maxLng = 180
minLat = -90
maxLat = 90

# create and execute the query
query = "SELECT lat, lng, intensity FROM mytable \
            WHERE %f<=tllat AND tllat<=%f \
            AND %f<=tllng AND tllng<=%f" % (minLat,maxLat,minLng,maxLng)

cur.execute(query)

pts = []
# print all the first cell of all the rows
for row in cur.fetchall() :
    print (row[1],row[0],row[2])
    # work out the mercator projection for latitute x = asinh(tan(x1))
    proj = math.degrees(math.asinh(math.tan(math.radians(row[0]))))
    print (row[1],proj,row[2])
    print "-"*15
    if (minLat < proj and proj < maxLat):
        pts.append((row[1],proj,row[2]))

print "Processing %d points..." % len(pts)

hm = heatmap.Heatmap()
hm.heatmap(pts, "bandwidth2.png",30,155,(1024,512),'fire',(minLng,maxLng,minLat,maxLat))

Answer 3

我认为一种方法是创建一个（更大的）元组列表，每个点根据该点的密度重复。密度高的点由许多点相互叠加，而密度低的点则很少。因此，您可以使用[(120.7, 82.5), (120.7, 82.5), (130.6, 81.5)]（一个相当枯燥的数据集）来代替：[(120.7, 82.5, 2), (130.6, 81.5, 1)]。

一个可能的问题是您的密度很可能是浮点数，而不是整数，因此您应该对数据进行规范化和四舍五入。进行转换的一种方法是这样的：

def dens2points (dens_tups):
    min_dens = dens_tups[0][2]
    for tup in dens_tups:
        if (min_dens > tup[2]):
           min_dens = tup[2]
    print min_dens

    result = []
    for tup in dens_tups:
        for i in range(int(tup[2]/min_dens)):
            result.append((tup[0],tup[1]))
    return result

if __name__ == "__main__":
    input = [(10, 10, 20.0),(5, 5, 10.0),(10,10,0.9)]
    output = dens2points(input)
    print input
    print output

（这不是很pythonic，但似乎适用于简单的测试用例）。此子例程应将您的数据转换为 heatmap.py 接受的形式。稍加努力，我认为子程序可以减少到两行。