用Python在地图上模拟疫情扩散

Dear 丶 2022-06-09 09:23 312阅读 0赞

用Python在地图上模拟疫情扩散

受杰森的《Almost Looks Like Work》启发，我来展示一些病毒传播模型。需要注意的是这个模型并不反映现实情况，因此不要误以为是西非可怕的传染病。相反，它更应该被看做是某种虚构的僵尸爆发现象。那么，让我们进入主题。

image

这就是SIR模型，其中字母S、I和R反映的是在僵尸疫情中，个体可能处于的不同状态。

  • S 代表易感群体，即健康个体中潜在的可能转变的数量。
  • I 代表染病群体，即僵尸数量。
  • R 代表移除量，即因死亡而退出游戏的僵尸数量，或者感染后又转回人类的数量。但对与僵尸不存在治愈者，所以我们就不要自我愚弄了（如果要把SIR模型应用到流感传染中，还是有治愈者的）。

至于β(beta)和γ(gamma):

  • β(beta)表示疾病的传染性程度，只要被咬就会感染。
  • γ(gamma)表示从僵尸走向死亡的速率，取决于僵尸猎人的平均工作速率，当然，这不是一个完美的模型，请对我保持耐心。

S′=−βIS告诉我们健康者变成僵尸的速率，S′是对时间的导数。

I′=βIS−γI告诉我们感染者是如何增加的，以及行尸进入移除态速率（双关语）。

R′=γI只是加上(gamma I)，这一项在前面的等式中是负的。

上面的模型没有考虑S/I/R的空间分布，下面来修正一下！

一种方法是把瑞典和北欧国家分割成网格，每个单元可以感染邻近单元，描述如下：

image

实验完整代码如下：

Main.py

  1. # -*- coding: utf-8 -*-
  2. import numpy as np
  3. import math
  4. import matplotlib.pyplot as plt    
  5. from matplotlib import rcParams
  6. import matplotlib.image as mpimg
  7. from PIL import Image
  9. rcParams['font.family'] = 'serif'
  10. rcParams['font.size'] = 16
  11. rcParams['figure.figsize'] = 12, 8
  14. beta = 0.010
  15. gamma = 1
  17. def euler_step(u, f, dt):
  18.     return u + dt * f(u)
  20. def f(u):
  21.     S = u[0]
  22.     I = u[1]
  23.     R = u[2]
  25.     new = np.array([-beta*(S[1:-1, 1:-1]*I[1:-1, 1:-1] + \
  26.                             S[0:-2, 1:-1]*I[0:-2, 1:-1] + \
  27.                             S[2:, 1:-1]*I[2:, 1:-1] + \
  28.                             S[1:-1, 0:-2]*I[1:-1, 0:-2] + \
  29.                             S[1:-1, 2:]*I[1:-1, 2:]),
  30.                      beta*(S[1:-1, 1:-1]*I[1:-1, 1:-1] + \
  31.                             S[0:-2, 1:-1]*I[0:-2, 1:-1] + \
  32.                             S[2:, 1:-1]*I[2:, 1:-1] + \
  33.                             S[1:-1, 0:-2]*I[1:-1, 0:-2] + \
  34.                             S[1:-1, 2:]*I[1:-1, 2:]) - gamma*I[1:-1, 1:-1],
  35.                      gamma*I[1:-1, 1:-1]
  36.                     ])
  38.     padding = np.zeros_like(u)
  39.     padding[:,1:-1,1:-1] = new
  40.     padding[0][padding[0] < 0] = 0
  41.     padding[0][padding[0] > 255] = 255
  42.     padding[1][padding[1] < 0] = 0
  43.     padding[1][padding[1] > 255] = 255
  44.     padding[2][padding[2] < 0] = 0
  45.     padding[2][padding[2] > 255] = 255
  47.     return padding
  49. img = Image.open('popdens2.png')
  50. img = img.resize((img.size[0]/2,img.size[1]/2)) 
  51. img = 255 - np.asarray(img)
  52. imgplot = plt.imshow(img)
  53. imgplot.set_interpolation('nearest')
  55. S_0 = img[:,:,1]
  56. I_0 = np.zeros_like(S_0)
  57. I_0[309,170] = 1 # patient zero
  59. R_0 = np.zeros_like(S_0)
  61. T = 900                         # final time
  62. dt = 1                          # time increment
  63. N = int(T/dt) + 1               # number of time-steps
  64. t = np.linspace(0.0, T, N)      # time discretization
  66. # initialize the array containing the solution for each time-step
  67. u = np.empty((N, 3, S_0.shape[0], S_0.shape[1]))
  68. u[0][0] = S_0
  69. u[0][1] = I_0
  70. u[0][2] = R_0
  72. import matplotlib.cm as cm
  73. theCM = cm.get_cmap("Reds")
  74. theCM._init()
  75. alphas = np.abs(np.linspace(0, 1, theCM.N))
  76. theCM._lut[:-3,-1] = alphas
  78. for n in range(N-1):
  79.     u[n+1] = euler_step(u[n], f, dt)
  81. from images2gif import writeGif
  83. keyFrames = []
  84. frames = 60.0
  86. for i in range(0, N-1, int(N/frames)):
  87.     imgplot = plt.imshow(img, vmin=0, vmax=255)
  88.     imgplot.set_interpolation("nearest")
  89.     imgplot = plt.imshow(u[i][1], vmin=0, cmap=theCM)
  90.     imgplot.set_interpolation("nearest")
  91.     filename = "outbreak" + str(i) + ".png"
  92.     plt.savefig(filename)
  93.     keyFrames.append(filename)
  95. images = [Image.open(fn) for fn in keyFrames]
  96. gifFilename = "outbreak.gif"
  97. writeGif(gifFilename, images, duration=0.3)
  98. plt.clf()

image2gif.py

  1. """ MODULE images2gif
  3. Provides a function (writeGif) to write animated gif from a series
  4. of PIL images or numpy arrays.
  6. This code is provided as is, and is free to use for all.
  8. Almar Klein (June 2009)
  10. - based on gifmaker (in the scripts folder of the source distribution of PIL)
  11. - based on gif file structure as provided by wikipedia
  13. """
  15. try:
  16.     import PIL
  17.     from PIL import Image, ImageChops
  18.     from PIL.GifImagePlugin import getheader, getdata
  19. except ImportError:
  20.     PIL = None
  22. try:
  23.     import numpy as np
  24. except ImportError:
  25.     np = None    
  27. # getheader gives a 87a header and a color palette (two elements in a list).
  28. # getdata()[0] gives the Image Descriptor up to (including) "LZW min code size".
  29. # getdatas()[1:] is the image data itself in chuncks of 256 bytes (well
  30. # technically the first byte says how many bytes follow, after which that
  31. # amount (max 255) follows).
  34. def intToBin(i):
  35.     """ Integer to two bytes """
  36.     # devide in two parts (bytes)
  37.     i1 = i % 256
  38.     i2 = int( i/256)
  39.     # make string (little endian)
  40.     return chr(i1) + chr(i2)
  43. def getheaderAnim(im):
  44.     """ Animation header. To replace the getheader()[0] """
  45.     bb = "GIF89a"
  46.     bb += intToBin(im.size[0])
  47.     bb += intToBin(im.size[1])
  48.     bb += "\x87\x00\x00"
  49.     return bb
  52. def getAppExt(loops=0):
  53.     """ Application extention. Part that secifies amount of loops. 
  54.     if loops is 0, if goes on infinitely.
  55.     """
  56.     bb = "\x21\xFF\x0B"  # application extension
  57.     bb += "NETSCAPE2.0"
  58.     bb += "\x03\x01"
  59.     if loops == 0:
  60.         loops = 2**16-1
  61.     bb += intToBin(loops)
  62.     bb += '\x00'  # end
  63.     return bb
  66. def getGraphicsControlExt(duration=0.1):
  67.     """ Graphics Control Extension. A sort of header at the start of
  68.     each image. Specifies transparancy and duration. """
  69.     bb = '\x21\xF9\x04'
  70.     bb += '\x08'  # no transparancy
  71.     bb += intToBin( int(duration*100) ) # in 100th of seconds
  72.     bb += '\x00'  # no transparant color
  73.     bb += '\x00'  # end
  74.     return bb
  77. def _writeGifToFile(fp, images, durations, loops):
  78.     """ Given a set of images writes the bytes to the specified stream.
  79.     """
  81.     # init
  82.     frames = 0
  83.     previous = None
  85.     for im in images:
  87.         if not previous:
  88.             # first image
  90.             # gather data
  91.             palette = getheader(im)[1]
  92.             data = getdata(im)
  93.             imdes, data = data[0], data[1:]            
  94.             header = getheaderAnim(im)
  95.             appext = getAppExt(loops)
  96.             graphext = getGraphicsControlExt(durations[0])
  98.             # write global header
  99.             fp.write(header)
  100.             fp.write(palette)
  101.             fp.write(appext)
  103.             # write image
  104.             fp.write(graphext)
  105.             fp.write(imdes)
  106.             for d in data:
  107.                 fp.write(d)
  109.         else:
  110.             # gather info (compress difference)              
  111.             data = getdata(im) 
  112.             imdes, data = data[0], data[1:]       
  113.             graphext = getGraphicsControlExt(durations[frames])
  115.             # write image
  116.             fp.write(graphext)
  117.             fp.write(imdes)
  118.             for d in data:
  119.                 fp.write(d)
  121. #             # delta frame - does not seem to work
  122. #             delta = ImageChops.subtract_modulo(im, previous)            
  123. #             bbox = delta.getbbox()
  124. #             
  125. #             if bbox:
  126. #                 
  127. #                 # gather info (compress difference)              
  128. #                 data = getdata(im.crop(bbox), offset = bbox[:2]) 
  129. #                 imdes, data = data[0], data[1:]       
  130. #                 graphext = getGraphicsControlExt(durations[frames])
  131. #                 
  132. #                 # write image
  133. #                 fp.write(graphext)
  134. #                 fp.write(imdes)
  135. #                 for d in data:
  136. #                     fp.write(d)
  137. #                 
  138. #             else:
  139. #                 # FIXME: what should we do in this case?
  140. #                 pass
  142.         # prepare for next round
  143.         previous = im.copy()        
  144.         frames = frames + 1
  146.     fp.write(";")  # end gif
  147.     return frames
  150. def writeGif(filename, images, duration=0.1, loops=0, dither=1):
  151.     """ writeGif(filename, images, duration=0.1, loops=0, dither=1)
  152.     Write an animated gif from the specified images. 
  153.     images should be a list of numpy arrays of PIL images.
  154.     Numpy images of type float should have pixels between 0 and 1.
  155.     Numpy images of other types are expected to have values between 0 and 255.
  156.     """
  158.     if PIL is None:
  159.         raise RuntimeError("Need PIL to write animated gif files.")
  161.     images2 = []
  163.     # convert to PIL
  164.     for im in images:
  166.         if isinstance(im,Image.Image):
  167.             images2.append( im.convert('P',dither=dither) )
  169.         elif np and isinstance(im, np.ndarray):
  170.             if im.dtype == np.uint8:
  171.                 pass
  172.             elif im.dtype in [np.float32, np.float64]:
  173.                 im = (im*255).astype(np.uint8)
  174.             else:
  175.                 im = im.astype(np.uint8)
  176.             # convert
  177.             if len(im.shape)==3 and im.shape[2]==3:
  178.                 im = Image.fromarray(im,'RGB').convert('P',dither=dither)
  179.             elif len(im.shape)==2:
  180.                 im = Image.fromarray(im,'L').convert('P',dither=dither)
  181.             else:
  182.                 raise ValueError("Array has invalid shape to be an image.")
  183.             images2.append(im)
  185.         else:
  186.             raise ValueError("Unknown image type.")
  188.     # check duration
  189.     if hasattr(duration, '__len__'):
  190.         if len(duration) == len(images2):
  191.             durations = [d for d in duration]
  192.         else:
  193.             raise ValueError("len(duration) doesn't match amount of images.")
  194.     else:
  195.         durations = [duration for im in images2]
  198.     # open file
  199.     fp = open(filename, 'wb')
  201.     # write
  202.     try:
  203.         n = _writeGifToFile(fp, images2, durations, loops)
  204.         print n, 'frames written'
  205.     finally:
  206.         fp.close()
  209. if __name__ == '__main__':
  210.     im = np.zeros((200,200), dtype=np.uint8)
  211.     im[10:30,:] = 100
  212.     im[:,80:120] = 255
  213.     im[-50:-40,:] = 50
  215.     images = [im*1.0, im*0.8, im*0.6, im*0.4, im*0]
  216.     writeGif('lala3.gif',images, duration=0.5, dither=0)

实验原始图像与实验后图像如下：

image

image

发表评论

表情：
评论列表 （有 0 条评论，312人围观）

还没有评论，来说两句吧...

相关阅读