我有大量的多边形(~100000),并尝试找到一种智能的方法来计算它们与常规网格单元的交叉区域.
目前,我正在使用形状(基于它们的角坐标)创建多边形和网格单元.然后,使用一个简单的for循环我遍历每个多边形并将其与附近的网格单元进行比较.
只是一个小例子来说明多边形/网格单元格.
from shapely.geometry import box, Polygon
# Example polygon
xy = [[130.21001, 27.200001], [129.52, 27.34], [129.45, 27.1], [130.13, 26.950001]]
polygon_shape = Polygon(xy)
# Example grid cell
gridcell_shape = box(129.5, -27.0, 129.75, 27.25)
# The intersection
polygon_shape.intersection(gridcell_shape).area
Run Code Online (Sandbox Code Playgroud)
(顺便说一句:网格单元的尺寸为0.25x0.25,多边形的最大尺寸为1x1)
实际上,对于具有大约0.003秒的单个多边形/网格单元组合,这是非常快的.但是,在大量多边形上运行此代码(每个多边形可以与几十个网格单元相交)在我的机器上花费大约15分钟以上(最多30分钟,具体取决于交叉网格单元的数量),这是不可接受的.不幸的是,我不知道如何编写多边形交集的代码来获得重叠区域.你有什么建议吗?有一种形状的替代品吗?
我想知道你是否可以用厘米来指定matplotlib中数字的大小.目前我写道:
def cm2inch(value):
return value/2.54
fig = plt.figure(figsize=(cm2inch(12.8), cm2inch(9.6)))
Run Code Online (Sandbox Code Playgroud)
但是有原生方法吗?
我有一个大约60000个形状的日期集(每个角的纬度/经度坐标),我想用matplotlib和底图在地图上绘制.
这就是我现在这样做的方式:
for ii in range(len(data)):
lons = np.array([data['lon1'][ii],data['lon3'][ii],data['lon4'][ii],data['lon2'][ii]],'f2')
lats = np.array([data['lat1'][ii],data['lat3'][ii],data['lat4'][ii],data['lat2'][ii]],'f2')
x,y = m(lons,lats)
poly = Polygon(zip(x,y),facecolor=colorval[ii],edgecolor='none')
plt.gca().add_patch(poly)
Run Code Online (Sandbox Code Playgroud)
但是,我的机器需要大约1.5分钟,我在想是否有可能加快速度.有没有更有效的方法来绘制多边形并将它们添加到地图中?
我正在使用以下代码在地图上绘制数据:
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from mpl_toolkits.basemap import Basemap
from matplotlib.patches import Polygon
from scipy.io import netcdf
ncfile = netcdf.netcdf_file(myfile.nc,'r')
lon = ncfile.variables['longitude'][:]
lat = ncfile.variables['latitude'][:]
data = ncfile.variables['mydata'][:]
ncfile.close()
m = Basemap(projection='nplaea', boundinglat=40, lon_0=270)
m.drawcoastlines(linewidth=.6, zorder=2)
m.drawparallels(np.arange(-80.,81.,20.), zorder=1)
m.drawmeridians(np.arange(-180.,181.,20.), zorder=1)
cNorm = mpl.colors.Normalize(vmin=0, vmax=np.nanmax(data))
cmap = plt.get_cmap('jet')
lons, lats = np.meshgrid(lon, lat)
x, y = m(lons, lats)
datamap = m.pcolor(x, y, data, zorder=0)
datamap.set_norm(cNorm)
plt.colorbar(datamap, cmap=cmap, norm=cNorm, …Run Code Online (Sandbox Code Playgroud) 我使用matplotlib.patches.Polygon类在地图上绘制多边形.实际上,给出了关于多边形角的坐标和每个"多边形"的浮点数据值的信息.现在,我想将这些数据值(范围从0到3e15)转换为颜色信息,以便将其可视化.在Python中执行此操作的最佳做法是什么?
我的代码片段:
poly = Polygon( xy, facecolor=data, edgecolor='none')
plt.gca().add_patch(poly)
Run Code Online (Sandbox Code Playgroud) 我有两个数组,比如varx和vary.两者都包含不同位置的NAN值.但是,我想对两者进行线性回归,以显示两个数组的相关程度.到目前为止,这非常有用:http://glowingpython.blogspot.de/2012/03/linear-regression-with-numpy.html
但是,使用这个:
slope, intercept, r_value, p_value, std_err = stats.linregress(varx, vary)
Run Code Online (Sandbox Code Playgroud)
导致每个输出变量的nans.将两个数组中的有效值作为线性回归的输入的最方便的方法是什么?我听说过屏蔽数组,但我不确定它是如何工作的.
执行一个python脚本(这里包括很长的方式)我写了一个警告消息.我不知道我的代码中的哪一行会被引发.我怎样才能获得这些信息?
此外,这究竟意味着什么?事实上,我不知道我使用的是某种蒙面数组?
/usr/lib/pymodules/python2.7/numpy/ma/core.py:3785: UserWarning: Warning: converting a masked element to nan.
warnings.warn("Warning: converting a masked element to nan.")
Run Code Online (Sandbox Code Playgroud) 我有数千个以表格形式存储的多边形(给定它们的4个角坐标),它们代表地球的小区域.另外,每个多边形都有一个数据值.该文件看起来像这样的例子:
lat1, lat2, lat3, lat4, lon1, lon2, lon3, lon4, data
57.27, 57.72, 57.68, 58.1, 151.58, 152.06, 150.27, 150.72, 13.45
56.96, 57.41, 57.36, 57.79, 151.24, 151.72, 149.95, 150.39, 56.24
57.33, 57.75, 57.69, 58.1, 150.06, 150.51, 148.82, 149.23, 24.52
56.65, 57.09, 57.05, 57.47, 150.91, 151.38, 149.63, 150.06, 38.24
57.01, 57.44, 57.38, 57.78, 149.74, 150.18, 148.5, 148.91, 84.25
...
Run Code Online (Sandbox Code Playgroud)
许多多边形相交或重叠.现在我想创建一个*m矩阵,范围从-90°到90°纬度和-180°到180°经度,例如0.25°x0.25°,以存储(面积加权)平均数据落在每个像素内的所有多边形的值.
因此,常规网格中的一个像素将获得一个或多个多边形的平均值(如果没有多边形与像素重叠,则为无).每个多边形应该根据其在该像素内的面积分数贡献该平均值.
基本上常规网格和多边形看起来像这样:

如果查看像素2,您会看到这个像素中有两个多边形.因此,我必须考虑它们的面积分数来取两个多边形的平均数据值.然后应将结果存储在常规网格像素中.
我环顾网络,到目前为止找不到令人满意的方法.由于我使用Python/Numpy进行日常工作,我想坚持下去.这可能吗?该软件包匀称看起来很有希望,但我不知道从哪里开始......一切移植到PostGIS的数据库的努力一个可怕的量,我想还会有我的方式相当多的障碍.
我有包含表示ISO格式日期的字符串列表的数据文件.目前,我正在阅读他们使用:
mydates = [ datetime.datetime.strptime(timdata[x], "%Y-%m-%dT%H:%M:%S") for x in range(len(timedata)) ]
Run Code Online (Sandbox Code Playgroud)
这看起来非常简单,但是在大约25000个日期的大型列表上运行时速度非常慢 - >每个转换列表大约0.34秒.由于我有数千个此类列表,我正在寻找更快的方法.但是,我还没找到.dateutil解析器的性能更差......
在具有4GB内存的计算机上,这种简单的插值会导致内存错误:
(基于:http://docs.scipy.org/doc/scipy/reference/tutorial/interpolate.html)
import numpy as np
from scipy.interpolate import interp1d
x = np.linspace(0, 10, 80000)
y = np.cos(-x**2/8.0)
f2 = interp1d(x, y, kind='cubic')
Run Code Online (Sandbox Code Playgroud)
我想过将数据切割成块,但有没有办法可以执行这种三次样条插值而不需要太多内存?为什么它甚至遇到麻烦?