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

Question

I'm working with an ERA5 (0.25x0.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.25x0.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?