scipy Map上的叠加插值

k75qkfdt  于 8个月前  发布在  其他
关注(0)|答案(1)|浏览(59)

我终于要寻求帮助了。在尝试了所有的选择之后,我被卡住了。
现在我有数据是从一个至少有105个村庄的地区的10个村庄抽样的。有了这10个村庄的样本数据,我做了预测,效果也很好,我的最终预测值表看起来像这样(不幸的是,我无法将此表转换为可以在这里分享的东西):

现在我的问题是插值。我想插入这些数据,以覆盖其他未采样的村庄,这是我如何做到的:

from scipy.interpolate import griddata

# Extract the longitude, latitude, and prediction columns from the decoded dataframe
interpolation_data = decoded_df[['longitude', 'latitude', 'prediction']]

# Remove any rows with missing values
interpolation_data = interpolation_data.dropna()

# Convert the data to numpy arrays
points = interpolation_data[['longitude', 'latitude']].values
values = interpolation_data['prediction'].values

# Define the grid points for interpolation
grid_points = np.vstack((grid_lon.flatten(), grid_lat.flatten())).T

# Perform IDW interpolation
interpolated_values = griddata(points, values, grid_points, method='linear')

interpolated_values = interpolated_values.reshape(grid_lon.shape)

# Create a contour plot of the interpolated predictions
plt.contourf(grid_lon, grid_lat, interpolated_values)
plt.colorbar()
plt.scatter(decoded_df['longitude'], decoded_df['latitude'], c=decoded_df['prediction'], cmap='viridis', edgecolors='black')
plt.xlabel('Longitude')
plt.ylabel('Latitude')
plt.title('Interpolated Predictions')
plt.show()

这给了我这个

现在,下一步是将插值结果叠加到该区域的Map上。我是这样做的:

import geopandas as gpd
from mpl_toolkits.axes_grid1 import make_axes_locatable

import geopandas as gpd
import matplotlib.pyplot as plt

# Read shapefile data of Babati
shapefile_path = "Babati Villages/Babati_villages.shp"  # Replace with the actual path to the shapefile
gdf_babati = gpd.read_file(shapefile_path)

gdf_bti= gdf_babati[gdf_babati["District_N"] == "Babati"]
gdf_bti.head()

# Define the grid points for interpolation
grid_points = np.vstack((grid_lon.flatten(), grid_lat.flatten())).T

# Perform IDW interpolation
interpolated_values = griddata(points, values, grid_points, method='linear')

# Reshape the interpolated values to match the grid shape
interpolated_values = interpolated_values.reshape(grid_lon.shape)

from shapely.geometry import box
# Create a bounding box geometry of the Babati region
bbox = box(gdf_bti.total_bounds[0], gdf_bti.total_bounds[1],
           gdf_bti.total_bounds[2], gdf_bti.total_bounds[3])

# Clip the interpolated predictions to the extent of the Babati region
interpolated_predictions = gpd.clip(interpolated_predictions, bbox)

# Create subplots
fig, ax = plt.subplots(figsize=(10, 10))

# Plot the shapefile of the Babati region
gdf_bti.plot(ax=ax, facecolor='none', edgecolor='black')

# Plot the interpolated predictions
interpolated_predictions.plot(ax=ax, column='prediction', cmap='viridis', markersize=30, legend=True)

# Add colorbar
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.1)
interpolated_predictions.plot(ax=cax, column='prediction', cmap='viridis', legend=True, cax=cax)

# Set plot title and labels
ax.set_title('Interpolated Predictions in Babati Region')
ax.set_xlabel('Longitude')
ax.set_ylabel('Latitude')

# Show the plot
plt.show()

这里是现在的问题是,因为插值的覆盖是完全关闭的。我希望它能覆盖所有插入的村庄,但它没有。这就是我得到的

我做错了什么,有什么办法解决这个问题吗?

f0ofjuux

f0ofjuux1#

  • 免责声明:我自己不是地理Maven,只习惯于用Python绘制Map,下面可能包括一些(希望如此!)小错误我很乐意纠正。*

我高度怀疑你的问题来自不管理CRS
实际上,原始数据集中的坐标是用经度和纬度表示的,这意味着这里的CRS是一个地理坐标系。简而言之,地球被近似为球体或椭球体,坐标给予其表面上的位置。
但是,您的shapefile是在投影坐标系中定义的,因此每个点的位置都用二维平面中的(x,y)元组表示。正如其名称所暗示的,这种类型的CRS定义了将球体/椭球体表面上的每个点Map到平面的过程(例如,投影),地理学家们想出了很多方法来做到这一点。

如何处理CRS?

您必须首先找出数据集和shapefile都在哪个CRS中表示。它通常表示为EPSG代码,您应该在找到数据集/shapefile的位置找到此信息。它也可能已经包含在您使用的文件中,在这种情况下,它可能已经存在于您的GeoDataFrame中,您可以使用以下命令找到它:

>>> gdf.crs
<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

如果它还没有,你需要设置它:

# set it in place...
gdf.set_crs(my_crs, inplace=True)

# ... or not
gdf = gdf.set_crs(my_crs)

正确定义CRS后,您只需决定要在哪个目标CRS中绘制Map,并将数据集和shapefile的CRS更改为该目标CRS。幸运的是,geopandas知道如何做到这一点:

# set it in place...
gdf.to_crs(target_crs, inplace=True)

# ... or not
gdf = gdf.to_crs(target_crs)

既然您的数据集和shapefile共享相同的CRS,那么您的图应该没问题!
最后一点注意:有一个合理的机会,你的数据集使用EPSG:4326和你的shapefile使用EPSG:3857

相关问题