英文:
From numpy to geotif having a georeferenced box
问题
I created a numpy array by calculating the density of dwellings within an area through the following code:
这是我通过以下代码计算区域内住宅密度而创建的NumPy数组：
def myplot(x, y, z, s, bins=10000):
heatmap, xedges, yedges = np.histogram2d(x, y, bins=bins, weights=z)
heatmap = gaussian_filter(heatmap, sigma=s)
extent = [xedges[0], xedges[-1], yedges[0], yedges[-1]]
return heatmap.T, extent
fig, axs = plt.subplots(2, 2)
# Generate some test data
x = buildings["x"]
y = buildings["y"]
weights = buildings["Area"]
sigmas = [0, 16, 32, 64]
for ax, s in zip(axs.flatten(), sigmas):
if s == 0:
ax.plot(x, y, weights, 'k.', markersize=5)
ax.set title("Scatter plot")
else:
img, extent = myplot(x, y, weights, s)
ax.imshow(img, extent=extent, origin='lower', cmap=cm.jet)
ax.set_title("Smoothing with $\sigma$ = %d" % s)
plt.savefig('export_'+str(s)+'.png', dpi=150, bbox_inches='tight')
plt.show()
This is the result and works fine:
这是结果，运行良好：
Now I need to save it as a geotif and I know the extreme coordinates of the box angles. I tried to do that using the following code:
现在我需要将其保存为geotif，并且我知道盒子角的极端坐标。我尝试使用以下代码来执行此操作：
# create a georeferenced box
transform = from_bounds(extent[0], extent[1], extent[2], extent[3], 10000, 10000)
# save the georeferenced tif
with rio.open('data.tif', 'w', driver='GTiff', height=10000, width=10000, count=1, dtype='float64', nodata=0, crs=32632, transform=transform) as dst:
dst.write(img, 1)
The problem is that the result is transpose and not in the right position. Could you help me to find the solution?
问题是结果是转置的，位置不正确。您能帮我找到解决方案吗？
英文:
I created a numpy array by calculating the density of dwellings within an area through the following code:
def myplot(x, y, z, s, bins=10000):
heatmap, xedges, yedges = np.histogram2d(x, y, bins=bins, weights=z)
heatmap = gaussian_filter(heatmap, sigma=s)
extent = [xedges[0], xedges[-1], yedges[0], yedges[-1]]
return heatmap.T, extent
fig, axs = plt.subplots(2, 2)
# Generate some test data
x = buildings["x"]
y = buildings["y"]
weights = buildings["Area"]
sigmas = [0, 16, 32, 64]
for ax, s in zip(axs.flatten(), sigmas):
if s == 0:
ax.plot(x, y, weights, 'k.', markersize=5)
ax.set_title("Scatter plot")
else:
img, extent = myplot(x, y, weights, s)
ax.imshow(img, extent=extent, origin='lower', cmap=cm.jet)
ax.set_title("Smoothing with $\sigma$ = %d" % s)
plt.savefig('export_'+str(s)+'.png', dpi=150, bbox_inches='tight')
plt.show()
This is the result and works fine:
enter image description here
Now I need to save it as a geotif and I know the extreme coordinates of the box angles. I tried to do that using the following code:
# create a georeferenced box
transform = from_bounds(extent[0], extent[1],extent[2], extent[3], 10000, 10000)
# save the georeferenced tif
with rio.open('data.tif', 'w', driver='GTiff', height=10000, width=10000, count=1, dtype='float64', nodata=0, crs=32632, transform=transform) as dst:
dst.write(img, 1)
The problem is that the result is transpose and not on the right position. Could you help me to find the solution?
I tried to develop the code but did not work
答案1
得分: 0
你只需在数组上简单使用 numpy.transpose - 这是一个非常快速的操作，不会复制数组。
GDAL使用传统的C样式栅格坐标。在numpy中，形状为（x，y）的数组表示x行y列像素，而在GDAL中则相反。
# 保存地理参考的tif
with rio.open('data.tif', 'w', driver='GTiff', height=10000, width=10000, count=1, dtype='float64', nodata=0, crs=32632, transform=transform) as dst:
dst.write(img.transpose(), 1)
英文:
You should simply use numpy.transpose on your array - it is a very fast operation that does not copy the array.
GDAL uses traditional C style raster coordinates. In numpy an array with shape (x, y) is x lines of y pixels, while in GDAL it is the other way around.
# save the georeferenced tif
with rio.open('data.tif', 'w', driver='GTiff', height=10000, width=10000, count=1, dtype='float64', nodata=0, crs=32632, transform=transform) as dst:
dst.write(img.tranpose(), 1)
通过集体智慧和协作来改善编程学习和解决问题的方式。致力于成为全球开发者共同参与的知识库，让每个人都能够通过互相帮助和分享经验来进步。
评论