Matplotlib 线段宽度的右对齐

Matplotlib 中Text有horizontalAlignment 和 verticalAlignment属性对应文字的水平对齐和垂直对齐方式,但是对于线段(plot),找遍全网却没有找到令线段宽度根据中法线左右对齐的解决方法。

好在通过查阅Matplotlib的Transformations教程,找到了使用ScaledTranslation另线段偏移(offset)的方法(https://matplotlib.org/3.2.2/tutorials/advanced/transforms_tutorial.html#using-offset-transforms-to-create-a-shadow-effect)。

通过简单的三角计算可以实现目标(图1)。

图1

直接上代码:

def segmentOffset(fig, x, y, linewidth, lr="left"):
    import numpy as np
    import matplotlib.transforms as transforms

    r = {"left": (-1, 1), "right": (1, -1)}[lr]
    delta = np.squeeze(np.diff([x, y]))
    del_p = np.hypot(*delta)
    dpi_scale = np.diag(fig.dpi_scale_trans.get_matrix())
    dx, dy = delta[::-1] * linewidth / dpi_scale[:2] / del_p / 2.0 * r
    offset = transforms.ScaledTranslation(dx, dy, fig.dpi_scale_trans)
    return offset

示例和图

fig, ax = plt.subplots(figsize=(5, 5))
x = [1, 4]
y = [1, 5]
plt.axis("equal")
lw = 10
ax.plot(x, y, color="k", ls=":", lw=2)
ax.plot(x, y, color="g", lw=lw, alpha=0.5, label="origin")
ax.plot(
    x,
    y,
    color="r",
    lw=lw,
    alpha=0.5,
    label="offset_left",
    transform=ax.transData + segmentOffset(fig, x, y, lw, lr="left"),
)
ax.legend()
图2

为什么要弄这个,因为之前在描最大风暴潮沿岸分布的时候,内部线的线宽在拐角地形处会盖过外边的线段。图3,图4

图3 线宽10,中对齐
图4 线宽5,右对齐(left offset),空白是因为网格边界和实际大陆边界有点差异。

还是有瑕疵但这也算性价比较高的解决方案了。

Python脚本下台风数据

从浙江水利厅的台风路径网下的,但是数据和中国气象网是一样的。

import requests
import json
import os
from urllib import parse as URL_PARSE

url_search = r"http://typhoon.zjwater.gov.cn/Api/TyphoonSearch/"
url_info   = r"http://typhoon.zjwater.gov.cn/Api/TyphoonInfo/"
output_path= r"./"

def getTYinfo(str):
    r = requests.get(url_info + str[:6])
    if r.status_code == 200:
        data = r.json()
        savePath = os.path.join(output_path + str + ".json")
        print("Typhoon Info ", str, "get, save in " + savePath)
    else:
        print("Typhoon Info ", str, "download failed.")
    
def main():
    s = input("Input search information: ")
    if not s:
        print("Incorrect input. Quiting")
        return
    r = requests.get(url_search + URL_PARSE.quote(s, encoding="GBK"))
    info = r.json()
    for i, sg  in enumerate(info):
        print(i, sg['name'])
    s = input("Select ids(seperate with [Space]):").split()
    spl = [int(v) for v in s]
    for i in spl:
        try:
            getTYinfo(info[i]['name'])
        except Exception as e:
            pass

if __name__ == "__main__":
    main()

shapely&非结构化网格:随机生成点产生三角形并融合

简单弄了稀疏化非结构网格的算法,再准备产生结构化网格点进行插值的时候碰要得到网格边界的问题,用shapely.union解决,下面用随机生成点进行示例。

import numpy as np
import matplotlib.pyplot as plt
from shapely.geometry import MultiPoint
from shapely.ops import triangulate
from descartes.patch import PolygonPatch
from functools import reduce

points_num = 15
x = np.random.rand(points_num)
y = np.random.rand(points_num)
points = MultiPoint(list(zip(x, y)))
triangles = triangulate(points)
plt.figure()
for triangle in triangles:
    patch = PolygonPatch(triangle, alpha=0.5, zorder=2)
    plt.gca().add_patch(patch)
for point in points:
    plt.plot(point.x, point.y, 'o', color='grey')
plt.xlim((0,1))
plt.ylim((0,1))
    
poly = reduce(lambda x, y: x.union(y), triangles)
plt.figure()
plt.gca().add_patch(PolygonPatch(poly, alpha=0.2))
outline = poly.boundary.coords.xy
plt.plot(*outline, 'r')
plt.xlim((0,1))
plt.ylim((0,1))

善用python第三方库:如何精确统计海南岛台风登陆点热点以及强度

海南岛由于其独特的地理位置,每年遭受大量来自南海或热带西太平洋的台风影响,素有“南海台风走廊”之称。建立海南岛的风暴潮模型对预防台风带来的风暴潮灾害有重要意义,而模型由台风风场驱动,如何根据不同的台风登陆点、来向、强度等因素设计敏感性试验很重要,因而我们需要根据历史资料对登陆海南岛的台风进行统计和分类。假设我们已整理好历年的台风6h数据,包括路径的经纬、强度、风速等,如何统计海南岛的台风登陆点分布呢?

首先,最容易想到的,也是前人使用较多的方法是,在研究区域及周边建立一个粗网格(如图1a所示),然后将台风路径插值到这些网格点中,对每个网格点进行统计计数,得到所有台风经过点的热力图,再根据这个图分析海南岛沿岸的登陆热点情况。这种方法虽然简单,但是缺点也很明显:过于粗糙、二次登陆和台风登出点均会对结果造成很大的误差。第二种方法是用几个矩形框将海南岛沿岸分为几个区域(如图1b),将台风路径进行加密插值,然后对每个台风路径按时间顺序遍历,判断路径点是否落在几个矩形区域内,若在区域内,进行记录,然后跳入下一个台风路径的遍历。这种方法可以避免对一次台风的多次统计,同时也对区域进行了划分,可以统计到台风的登陆热点分布情况。然而,这种方法依然不是我们想要的,不够精确,也不够优雅。

第二种方法里有“判断点在区域内”的思路,不过是判断矩形区域的,只需判断经纬坐标的范围就行,方法比较简单。我们可以在这基础上扩展一下思路:如何判断点在多边形内(Polygon Inside)。判断点在多边形内有几种方法,面积判别、夹角判别或射线法,但是自己再写这个算法显然太麻烦,这么常见的功能应该有python第三方包。上网查阅资料后发现Shapely(https://github.com/Toblerity/Shapely)符合我们的要求,它是用于处理和分析平面几何对象的一个python第三方库,基于广泛运用的GEOS(PostGIS的引擎)和JTS库,具有强大的平面几何图形(点、线、多边形等)处理和运算功能。有了Shapely,我们还需要海南岛的轮廓线。GADM(http://www.gadm.org/)提供全球的行政区划数据库,包括了几乎所有国家和地区的国界、子级行政区域划分数据,其提供的shapefile格式数据可以直接用于python的Basemap包(https://matplotlib.org/basemap/)的绘图中,也包含了所有区域轮廓的经纬坐标点。GADM提供的中国行政区域划分中分为0,1,2,3级数据,对应国、省、县(市)、乡的行政区划精度。在获得以上的数据和扩展包之后,我们便可以对台风进行登陆点统计。

首先,从GADM中国1级数据中获取海南岛的轮廓线,其包含35749个坐标点,数据量过于庞大,我们需要先用Shapely的Polygon生成海南岛轮廓线的多边形,剔除掉面积过小的区划(海南周边海岛),然后用.simplify(0.05)将多边形进行简化为49个坐标点的多边形(如图1c)。对每个台风,将台风路径用Shapely的LineString生成折线段,用Polygon.intersection与LineString进行相交得到交点。通过试验得知,Shapely中Polygon与LineString类的Intersection是会考虑LineString的坐标顺序的,因而只要LineString中的坐标顺序是按台风路径的时间顺序,我们只需要取intersection的第一个交点即可。得到所有台风的登陆点位置后,我们用GADM的3级数据的乡级行政区划来进行统计(图1d),用每个区划的轮廓经纬生成Polygon,用Polygon.contains(Point)来判断登陆点是否在区域内。但是,由于我们上面生成的海南岛的轮廓线是简化过的,登陆点可能落在区域外,因而,我们需要换个思路:用Point.buffer(0.04)将台风登陆点稍微地放大成圆形区域,然后再用Polygon.intersects(..)判断区域与圆形是否相交,若相交,则进行统计,并在接下来的区划遍历中移除该台风(防止buffered后的点与多个区划相交造成重复统计)。最后,使用Basemap将台风登陆点的分布画出,使用windrose的python包(https://github.com/python-windrose/windrose)将统计得到的数据进行展示(图2展示台风强度的风玫瑰图)。

总结:Python是一门优雅且灵活的脚本编程语言,其数量众多的丰富成熟的第三方库可以大大提高我们的工作效率。在解决实际工作问题的时候,要善于在网上寻找是否有第三方库可提供成熟的算法或解决方案,阅读这些库的Api说明文档,从中找到可行的解决路线。同样,互联网上也有许多公开的数据库或数据集,有效的利用这些数据也可以使我们的工作更加精确和高效。

Image copy from Clipboard

The win32clipboard module provides method of getting content from clipboard. But, when it comes to image content, bytes data are returned and are hard to convert to image type(involving Bitmap Header and ‘PIL.Image.frombytes() ‘ requires size though is unknown).

However, PIL module offers function can grab image from clipboard.

# from PIL/ImageGrab.py
...
def grabclipboard():
    if sys.platform == "darwin":
        fh, filepath = tempfile.mkstemp('.jpg')
        os.close(fh)
        commands = [
            "set theFile to (open for access POSIX file \""+filepath+"\" with write permission)",
            "try",
                "write (the clipboard as JPEG picture) to theFile",
            "end try",
            "close access theFile"
        ]
        script = ["osascript"]
        for command in commands:
            script += ["-e", command]
        subprocess.call(script)

        im = None
        if os.stat(filepath).st_size != 0:
            im = Image.open(filepath)
            im.load()
        os.unlink(filepath)
        return im
    else:
        data = Image.core.grabclipboard()
        if isinstance(data, bytes):
            from . import BmpImagePlugin
            import io
            return BmpImagePlugin.DibImageFile(io.BytesIO(data))
        return data

So, ‘BmImagePlugin.DibImageFile()’ is exactly what we need.

import win32clipboard as w
from PIL import Image, BmpImagePlugin
import io
#
w.OpenClipboard()
if w.IsClipboardFormatAvailable(w.CF_DIB):
    c = w.GetClipboardData(w.CF_DIB)
    im = BmpImagePlugin.DibImageFile(io.BytesIO(c))
    im.show() # or other operations
w.CloseClipboard()

多线程批量下载ASCAT风场数据

ftp进不去,只能用http接口这样子。
主页:http://www.remss.com/missions/ascat/
说明:http://data.remss.com/ascat/metopa/readme_ASCAT_metopA.txt
bytemap格式数据:http://data.remss.com/ascat/metopa/bmaps_v02.1/

from bs4 import BeautifulSoup
import requests
import threading
import os
import time

threadslimit = 50
global count
count = 0

domain = r"http://data.remss.com"
urlprefix = r"/ascat/metopa/bmaps_v02.1/y2018/m"

def download(url):
    global count
    count += 1
    id = count
    filename = url.split('/')[-1]
    print("【%03d】Downloading: %s" % (id, filename))
    r = requests.get(url, stream=True)
    f = open(os.path.join(os.getcwd(),'download', filename), 'wb')
    for chunk in r.iter_content(chunk_size=1024*1024):
        if chunk:
            f.write(chunk)
    f.close()
    print("【%03d】download completed." % id)
  
threads = []
for m in range(1, 13):
    url = domain + urlprefix + "%02d" % m
    r = requests.get(url)
    r.encoding = "utf-8"
    soup = BeautifulSoup(r.text)
    print(m)
    for item in soup.select("a"):
        if item['href'].endswith('1.gz'):
            downlink = domain + item['href']
            threads.append(threading.Thread(target=download, args=(downlink, )))
    r.close()

print('Threads initialization complited, downloading...')
for t in threads:
    t.start()
    while True:
        if(len(threading.enumerate()) <= threadslimit):
            break
        time.sleep(0.1)

记笔记——Python、Matplotlib、Basemap

Basemap:

 

basemap在初始化地图的时候实在是太慢了,如果需要画多个图,速度感人。解决思路是把一个画完图的对象存进内存(做个备份),然后再需要新的一张图的时候再读取这个图,需要注意的是deepcopy并不管用,需要更深层次的pickle.dumps和pickle.loads。 代码大致是:

def plotbasemap(self, ...):
    if not self.m:
    self.fig = plt.figure(figsize=(xx,xx))
    map = Basemap(...)
    ...
    self.m = pickle.dumps(map)
    return map
return pickle.loads(self.m)

值得注意的是,这种方法也可以用来/快速/保存读取变量。虽然占用空间大,但好在比较便捷,省空间的存法推荐用numpy存成npy。

import pickle
Dataset = ...
pickle.dump(Dataset, open('save.txt','wb')
data = pickle.load(open('save.txt','wb')

数据转化, DictInList -> ListInDict:

如果你有一个数据格式类似
[
    {
        'a' : 1,
        'b' : -1
    },
    {
        'a' : 2,
        'b' : -2
    },
    {
        'a' : 3,
        'b' : -3
    }
]

想转化为

{
    'a' : [1, 2, 3], 
    'b' : [-1, -2, -3]
}

 

def dataFormatConvert(data):
    return { key : [ p[key] for p in data ] for key in data[0].keys() } if type(data) == list else ...