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

Question

I am new to using python for processing rasters. I have an interpolated map that I have divided into a grid system. I need to overlay a field boundary (shapefile) and clip it to be within the boundary. I would also like to display it around the field after clipping. However, I am struggling with how to do this. Both the input data (.csv file) and shapefile can be found here. This code is run on Jupyter Lab, so I have blocked it into corresponding code chunks.

Here is the code to reproduce the project:

#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

Set graph parameters

#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)

Cubic interpolation

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

Final product

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

After printing out the plot, I need to overlay it with the shapefile and clip it so that the contents of the graph are only inside the boundary.

Answer 1

if you can save your map_cubic as geotiff you can use fiona and rasterio.mask to do the clipping:

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)

for writing ndarray to geotiff see for example Simplest way to save array into raster file in Python