轮廓曲线数据（lat，lon，值）在边界内并导出GeoJSON [英] Contour plot data (lat,lon,value) within boundaries and export GeoJSON
问题描述
我尝试在边界内插入数据，并根据纬度，经度，csv文件中的值数据绘制轮廓线（多边形）。
as 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
z = data [：3];
plt.contour（x，y，z）
plt.show（）
b $ b
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）
我想你的代码有一些错误（根据你的数据，你不应该做 x = data [：1] 但更多 x = data [...，1] ）。
使用你的数据，我将遵循插入 z 值和获取一个输出作为geojson的基本步骤将需要至少< a href =https://pypi.python.org/pypi/Shapely =nofollow> shapely 模块（此处 geopandas 用于方便）。
import numpy as np
import matplotlib.pyplot as plt
从geopandas import GeoDataFrame
从matplotlib.mlab import griddata
from shapely.geometry import Polygon，MultiPolygon
def collec_to_gdf（collec_poly）：
将一个`matplotlib.contour.QuadContourSet`转换为一个GeoDataFrame
多边形，colors = []，[]
对于i，多边形在枚举（collec_poly.collections）：
mpoly = []
对于polygon.get_paths
try：
path.should_simplify = False
poly = path.to_polygons（）
＃每个多边形都应该包含外环+可能的孔：
，holes = []，[]
if len（poly）> 0和len（poly [0]）> 3：
＃列表的第一个是外部环：
exterior = poly [0]
＃其他是孔：
如果len ）> 1：
holes = [h for poly in poly [1：] if len（h）> 3
mpoly.append（polygon（exterior，holes））
except：
print（'警告：制作多边形时出现几何错误＃{}'
.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 b geometry = polygons，
data = {'RGBA'：colors}，
crs = {'init'：'epsg：4326'}）
= 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 ），200）
zi = griddata（x，y，z，xi，yi，interp ='linear'）＃你也可以看看scipy.interpolate.griddata
nb_class = 15＃设置轮廓创建类的数量
＃类也可以通过它们的限制来定义，如[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 ）
＃以GeoJSON形式输出特征集合：
＃'{type：FeatureCollection，features：[{type：Feature，geometry：{ type：Polygon，coordinates：[[[51.563214073747474，
＃（...）
Edit：
您可以通过使用
来获取matplotplib使用的颜色值（作为范围0-1中的4个值的数组，表示RGBA值） c $ c> get_facecolor（）方法（然后使用它们填充您的GeoDataFrame的一个（或4）列：
colors = [p.get_facecolor（）。tolist（）[0] for p in collec_poly.collections]
gdf ['RGBA'] = colors
一旦你拥有GeoDataFrame，你就可以很容易地得到它与你的边界相交。使用这些边界来形成多边形并计算与您的结果的交集：
mask = Polygon ），（51,18），（52,18），（52,17），（51,17）]）
gdf.geometry = gdf.geometry.intersection（mask）
或以GeoDataFrame的形式读取geojson：
from shapely.ops import unary_union，polygonize
boundary = GeoDataFrame.from_file（'your_geojson'）
＃假设你有一个边界作为linestring，将它转换为多边形：
mask_geom = unary_union（[i for polygon in polygonize（boundary.geometry）]）
＃然后计算交集：
gdf.geometry = gdf.geometry.intersection（mask_geom）
I'm trying to interpolate data within boundaries and plot contour lines (polygons) based on Latitude, longitude, value data from csv files.
Results I want print as geojson.
I'm stuck on the basic contour plot from csv. I really appreciate help here.
This is what I got in this moment but can't make it work.
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()
The CSV looks like this
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
Column 1 - Latitude Column 2 - Longitude Column 3 - Value
For contour lines I need also limits - for example (4.1,4.3,4.6)
I guess there is some mistake in your code (according to your data you shouldn't do x = data[:1] but more x = data[..., 1]).
With your of data, the basic steps I will follow to interpolate the z value and fetch an output as a geojson would require at least the shapely module (and here geopandas is used for the convenience).
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,
# (...)
Edit:
You can grab the colors values (as an array of 4 values in range 0-1, representing RGBA values) used by matplotplib by fetching them on each item of the collection with the get_facecolor() method (and then use them to populate one (or 4) column of your GeoDataFrame :
colors = [p.get_facecolor().tolist()[0] for p in collec_poly.collections]
gdf['RGBA'] = colors
Once you have you GeoDataFrame you can easily get it intersected with your boundaries. Use these boundaries to make a Polygon with shapely and compute the intersection with your result:
mask = Polygon([(51,17), (51,18), (52,18), (52,17), (51,17)])
gdf.geometry = gdf.geometry.intersection(mask)
Or read your geojson as a 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)
这篇关于轮廓曲线数据（lat，lon，值）在边界内并导出GeoJSON的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
