Getting the coordinates of the points inside of a shapefile/Getting the coordinates of the closest point if there are none

I’m working with an ERA5 (0.25×0.25 degree) temperature dataset for Europe in order to get average temperature in small statistical regions called NUTS 3. So far I’ve used a shapefile of the NUTS 3 regions to visualize the temperature on an specific day and time.I also plotted a 0.25×0.25 degree grid. The code can be seen below:

from netCDF4 import Dataset
import numpy as np
import matplotlib.pyplot as plt
import xarray as xr
from mpl_toolkits.basemap import Basemap

ds=xr.load_dataset('Europe_temp_dataset-test.nc')
ds_celsius=ds-273
ds_celsius.t2m.attrs['units']='Celsius'
df=ds_celsius.to_dataframe()
df.reset_index(inplace=True)
df=df.set_index('time')
df.to_csv('Europe_temperature_dataset.csv')


ds_celsius.to_netcdf('Europe_temp_in_celsius_dataset-test4.nc')  #To solve: If file already exists, it gives [Errno 13] error


data=Dataset('C:/Users/sand_jv/Desktop/NetCDF4/Europe_temp_in_celsius_dataset-test4.nc')
lats=data.variables['latitude'][:]
lons=data.variables['longitude'][:]
time=data.variables['time'][:]
temp=data.variables['t2m'][:]
    
map=Basemap(projection='merc',
            llcrnrlon=-11.0,
            llcrnrlat=28.0,
            urcrnrlon=32.0,
            urcrnrlat=72.0,
            resolution='i'
            )

lon,lat=np.meshgrid(lons,lats)
x,y=map(lon,lat)
    
c_scheme=map.pcolor(x,y,np.squeeze(temp[0,:,:]),cmap='jet',shading='auto')
    #map.drawcoastlines()
    #map.drawcountries(linewidth=1,color='k')
map.drawmeridians(np.arange(-11,32,0.25))
map.drawparallels(np.arange(28,72,0.25))
    
map.readshapefile('NUTS_BN_01M_2016_4326','NUTS 3',color='r')

cbar=map.colorbar(c_scheme,location='right', pad='10%')
plt.title('Temperature 2m, 01.01.2019, 00:00')
plt.show()

The results can be seen in the next picture:

As an example, zooming in, the plot looks like this:

As it can be observed in the second picture, there are two cases:

1. Regions with one or more grid points inside
2. Regions with no grid points inside (circled in green in the example)

For case 1, I would like to get the coordinates of every grid point contained in each region so I can get an average of the temperature. For case 2, if there are no grid points inside of the region, I would like to find the coordinates of the closest grid point.

Any ideas on how to achieve this?