python - 边界内的等高线图数据(纬度、经度、值)并导出 GeoJSON

Question

我正在尝试根据来自 csv 文件的纬度、经度和值数据在边界内插入数据并绘制等高线(多边形)。

我想将结果打印为 geojson。

我被困在 csv 的基本等高线图上。非常感谢您的帮助。

这就是我此刻得到的，但无法让它发挥作用。

import numpy as np
import matplotlib.pyplot as plt


data = np.genfromtxt('temp.csv', delimiter=',')

x = data[:1]
y = data[:2]
[x,y] = meshgrid(x,y);
z = data[:3];

plt.contour(x,y,z)
plt.show()

CSV 看起来像这样

2015-10-28T09:25:12.800Z,51.56325,17.529043333,4.6484375,19.8046875
2015-10-28T09:25:12.900Z,51.56325,17.529041667,4.4921875,19.6875
2015-10-28T09:25:13.000Z,51.563248333,17.529041667,4.453125,19.9609375
2015-10-28T09:25:13.200Z,51.563245,17.529041667,4.140625,19.765625
2015-10-28T09:25:13.300Z,51.563243333,17.529041667,4.609375,19.6484375
2015-10-28T09:25:13.500Z,51.563241667,17.529041667,4.609375,19.53125
2015-10-28T09:25:13.600Z,51.56324,17.529041667,4.4921875,19.375
2015-10-28T09:25:13.700Z,51.563238333,17.529041667,4.4140625,19.765625
2015-10-28T09:25:13.800Z,51.563236667,17.529041667,4.453125,20.234375
2015-10-28T09:25:13.900Z,51.563235,17.529041667,4.3359375,19.84375
2015-10-28T09:25:14.000Z,51.563233333,17.529041667,4.53125,19.453125
2015-10-28T09:25:14.100Z,51.563231667,17.529041667,4.53125,19.8046875
2015-10-28T09:25:14.200Z,51.563228333,17.529041667,4.1796875,19.4921875
2015-10-28T09:25:14.300Z,51.563226667,17.529041667,4.2578125,19.453125
2015-10-28T09:25:14.400Z,51.563225,17.529041667,4.4921875,19.4921875
2015-10-28T09:25:14.500Z,51.563223333,17.529041667,4.375,19.453125
2015-10-28T09:25:14.600Z,51.563221667,17.529041667,4.609375,18.90625
2015-10-28T09:25:14.700Z,51.563218333,17.529041667,4.53125,19.6875
2015-10-28T09:25:14.900Z,51.563215,17.529041667,4.140625,18.75
2015-10-28T09:25:15.000Z,51.563213333,17.529041667,4.453125,18.828125

第 1 列 - 纬度第 2 列 - 经度第 3 列 - 值

对于轮廓线，我还需要限制 - 例如 (4.1,4.3,4.6)

Answer 1

我猜你的代码有一些错误(根据你的数据你不应该做 x = data[:1] 但更多 x = data[..., 1 ]).

根据您的数据，我将按照基本步骤插入 z 值并获取输出作为 geojson 至少需要 shapely模块(这里使用 geopandas 是为了方便起见)。

import numpy as np
import matplotlib.pyplot as plt
from geopandas import GeoDataFrame
from matplotlib.mlab import griddata
from shapely.geometry import Polygon, MultiPolygon

def collec_to_gdf(collec_poly):
    """Transform a `matplotlib.contour.QuadContourSet` to a GeoDataFrame"""
    polygons, colors = [], []
    for i, polygon in enumerate(collec_poly.collections):
        mpoly = []
        for path in polygon.get_paths():
            try:
                path.should_simplify = False
                poly = path.to_polygons()
                # Each polygon should contain an exterior ring + maybe hole(s):
                exterior, holes = [], []
                if len(poly) > 0 and len(poly[0]) > 3:
                    # The first of the list is the exterior ring :
                    exterior = poly[0]
                    # Other(s) are hole(s):
                    if len(poly) > 1:
                        holes = [h for h in poly[1:] if len(h) > 3]
                mpoly.append(Polygon(exterior, holes))
            except:
                print('Warning: Geometry error when making polygon #{}'
                      .format(i))
        if len(mpoly) > 1:
            mpoly = MultiPolygon(mpoly)
            polygons.append(mpoly)
            colors.append(polygon.get_facecolor().tolist()[0])
        elif len(mpoly) == 1:
            polygons.append(mpoly[0])
            colors.append(polygon.get_facecolor().tolist()[0])
    return GeoDataFrame(
        geometry=polygons,
        data={'RGBA': colors},
        crs={'init': 'epsg:4326'})


data = np.genfromtxt('temp.csv', delimiter=',')
x = data[..., 1]
y = data[..., 2]
z = data[..., 3]
xi = np.linspace(x.min(), x.max(), 200)
yi = np.linspace(y.min(), y.max(), 200)
zi = griddata(x, y, z, xi, yi, interp='linear') # You could also take a look to scipy.interpolate.griddata

nb_class = 15 # Set the number of class for contour creation
# The class can also be defined by their limits like [0, 122, 333]
collec_poly = plt.contourf(
    xi, yi, zi, nb_class, vmax=abs(zi).max(), vmin=-abs(zi).max())

gdf = collec_to_gdf(collec_poly)
gdf.to_json()
# Output your collection of features as a GeoJSON:
# '{"type": "FeatureCollection", "features": [{"type": "Feature", "geometry": {"type": "Polygon", "coordinates": [[[51.563214073747474,
# (...)

编辑:您可以通过使用 get_facecolor() 方法(和然后使用它们填充 GeoDataFrame 的一个(或 4 个)列:

colors = [p.get_facecolor().tolist()[0] for p in collec_poly.collections]
gdf['RGBA'] = colors

一旦你有了 GeoDataFrame，你就可以很容易地让它与你的边界相交。使用这些边界来制作一个具有形状的多边形并计算与你的结果的交集:

mask = Polygon([(51,17), (51,18), (52,18), (52,17), (51,17)])
gdf.geometry = gdf.geometry.intersection(mask)

或者将您的 geojson 读取为 GeoDataFrame:

from shapely.ops import unary_union, polygonize
boundary = GeoDataFrame.from_file('your_geojson')
# Assuming you have a boundary as linestring, transform it to a Polygon:
mask_geom = unary_union([i for i in polygonize(boundary.geometry)])
# Then compute the intersection:
gdf.geometry = gdf.geometry.intersection(mask_geom)