边界内的轮廓图数据（lat，lon，value）并导出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, 
# (...) 

编辑： 可以抓取的颜色值（如4个值的范围在0-1数组，表示RGBA值）通过获取他们在与get_facecolor()方法集合的每个项目由matplotplib使用（然后用它们来填充一个（或4）您GeoDataFrame列：

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) 

Answer 2

我几乎得到了我除外。根据mgc答案。我建议您将scipy.interpolate.griddata和插值方法更改为最接近的griddata。现在它的作品就像我想要的。

只需要最后一件事情就是限制插值到geoJson的多边形（边界）。

另一个问题是将颜色导出为geojson多边形作为属性。

import numpy as np 
import matplotlib.pyplot as plt 
#from matplotlib.mlab import griddata 
from scipy.interpolate import griddata 

from geopandas import GeoDataFrame 
from shapely.geometry import Polygon, MultiPolygon 

def collec_to_gdf(collec_poly): 
    """Transform a `matplotlib.contour.QuadContourSet` to a GeoDataFrame""" 
    polygons = [] 
    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) 
     elif len(mpoly) == 1: 
      polygons.append(mpoly[0]) 
    return GeoDataFrame(geometry=polygons, crs={'init': 'epsg:4326'}) 

data = np.genfromtxt('temp2.csv', delimiter=',') 
x = data[..., 1] 
y = data[..., 2] 
z = data[..., 4] 
xi = np.linspace(x.min(), x.max(), num=100) 
yi = np.linspace(y.min(), y.max(), num=100) 

#zi = griddata(x, y, z, xi, yi, interp='nn') # You could also take a look to scipy.interpolate.griddata 
zi = griddata((x, y), z, (xi[None,:], yi[:,None]), method='nearest') 

nb_class = [5,10,15,20,25,30,35,40,45,50] # 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() 
print gdf.to_json() 

plt.plot(x,y) 

plt.show()