Python

Question

我不熟悉使用 python 处理栅格。我有一个已划分为网格系统的插值地图。我需要覆盖一个字段边界（shapefile）并将其剪辑到边界内。我还想在裁剪后在场地周围显示它。但是，我正在为如何做到这一点而苦苦挣扎。输入数据（.csv 文件）和 shapefile can be found here。此代码在 Jupyter Lab 上运行，因此我已将其分块到相应的代码块中。

这里是重现项目的代码：

#Import modules
import pandas as pd
import numpy as np
import os
import shapely
import geopandas as geo
import glob
import holoviews as hv
hv.extension('bokeh')
from scipy.interpolate import griddata, interp2d
import fiona
import gdal
import ogr

#Read files
nut = pd.read_csv('nutrition.csv') #Data to be interpolated
boundary = fiona.open('field_shapefile.shp') #Shapefile, border around field

#Visualize map points to understand overlay of interpolated points
map_points = hv.Scatter(nut, 'longitude', ['latitude', 'n_ppm']).opts(size = 8,
                                                                      line_color = 'black',
                                                                      colorbar = True,
                                                                      cmap = 'Reds',
                                                                      aspect = 'equal',
                                                                      title = 'Title')

map_points

设置图形参数

#Minimum and maximum longtude values
lon_min = nut['longitude'].min()
lon_max = nut['longitude'].max()

#Create range of nitrogen values
lon_vec = np.linspace(lon_min, lon_max, 30) #Set grid size

#Find min and max latitude values
lat_min = nut['latitude'].min()
lat_max = nut['latitude'].max()

# Create a range of N values spanning the range from min to max latitude
# Inverse the order of lat_min and lat_max to get the grid correctly
lat_vec = np.linspace(lat_max,lat_min,30,)

# Generate the grid
lon_grid, lat_grid = np.meshgrid(lon_vec,lat_vec)

三次插值

points = (nut['longitude'], nut['latitude'])
targets = (lon_grid, lat_grid)
grid_cubic = griddata(points, nut['n_ppm'], targets, method='cubic', fill_value=np.nan)

最终产品

map_bounds=(lon_min,lat_min,lon_max,lat_max)


map_cubic = hv.Image(grid_cubic, bounds=map_bounds).opts(aspect='equal',
                                                         colorbar=True,
                                                         title='Cubic',
                                                         xlabel='Longitude',
                                                         ylabel='Latitude',
                                                         cmap='Reds')

map_cubic

打印出图后，我需要用 shapefile 覆盖它并对其进行裁剪，使图形的内容仅在边界内。

Answer 1

如果您可以将 map_cubic 保存为 geotiff，您可以使用 fiona 和 rasterio.mask 进行剪辑：

with fiona.open('field_shapefile.shp', "r") as shapefile:
    shapes = [feature["geometry"] for feature in shapefile]

with rasterio.open('map_cubic.tif') as src:
    mc, mc_transform = rasterio.mask.mask(src, shapes, nodata=np.nan, crop=False)

用于将 ndarray 写入 geotiff，例如参见 [=15=]

Python - 如何使用插值 python 地图叠加和裁剪 shapefile？

Python - How to overlay and clip a shapefile with an interpolated python map?

clip