2016-01-19 62 views
0

我試圖在csv文件中基於緯度,經度,值數據的邊界內插入數據並繪製輪廓線(多邊形)。邊界內的輪廓圖數據(lat,lon,value)並導出GeoJSON

結果我想打印爲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)

+0

我想你的代碼和geojson之間有幾步作爲輸出。運行代碼時是否遇到過一些錯誤消息? – mgc

回答

3

我想你的代碼有一些錯誤(根據你的數據你不應該這樣做)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) 
+0

謝謝,這會幫助我很多。我試圖首先繪製結果,但如果我使用所有關於40k點的數據,則會出現錯誤。 griddata中的錯誤(RuntimeError:初始化:Triangulation無效)。也許我需要先排序這些數據?但是如何? – seek

+0

我在geojson文件中將邊界作爲多邊形(有時甚至是帶孔的多邊形)。 geometry.intersection給了我很奇怪的結果。裏面只顯示最小的多邊形。 關於RGBA顏色我有錯誤:「ValueError:值的長度與索引長度不匹配」在gdf ['RGBA'] =顏色 – seek

+0

對不起。我沒有在我的邊界多邊形上使用這個unary_union函數。現在它看起來像它的工作..如果你可以檢查爲什麼這種顏色不適合我的工作會很好。 – seek

0

我幾乎得到了我除外。根據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() 
+0

我編輯了我關於顏色和邊界的答案。 – mgc