I changed your model solution a little bit; but it's like <PERSON> wrote in his answer: you redefined xi and yi in the process, so when you later on define the netCDF values for lat and lon you try to fill small 1D arrays with a 2D array. That is the source of your error.
```python import numpy as np from scipy.interpolate import griddata import xarray as xr import pandas as pd import rioxarray import netCDF4 as nc4 import matplotlib.pyplot as plt from shapely.geometry import Point import geopandas as gpd from pyproj import CRS
CRS_3857 = CRS.from_epsg(3857) CRS_wkt = CRS_3857.to_wkt(pretty=False) print(CRS_wkt)
--- Input data....
lat = [50.1, 50.2, 50.3, 50.4, 50.5, 62] lon = [ 8.1, 8.2, 8.3, 8.4, 8.5, 12] temp = [1 , 2 , 3 , 4 , 5 , 6]
Put into pandas Dataframe
df = pd.DataFrame( { 'latitude': lat, 'longitude': lon, 'temp': temp } )
Prepare geometry
pointShp = [Point(x, y) for x, y in zip(df.longitude, df.latitude)] pointGpd = gpd.GeoDataFrame(df, geometry=pointShp, crs='EPSG:4326')
Reproject
point3857 = pointGpd.to_crs('EPSG:3857') point3857['x'] = point3857.apply(lambda x: x.geometry.centroid.x, axis=1) point3857['y'] = point3857.apply(lambda x: x.geometry.centroid.y, axis=1) df = point3857[['x','y','temp']] lon = list(df['x']) lat = list(df['y']) temp = list(df['temp'])
Proceed with your source code
--- Project input data on a regular grid
xi = np.linspace(min(lon), max(lon), 500) yi = np.linspace(min(lat), max(lat), 500) xi, yi = np.meshgrid(xi, yi) zi = np.ones_like(xi,dtype=np.float32) * np.NaN
zi = griddata((lon, lat), temp, (xi, yi), method='linear')
for i in range(len(temp)): idx = np.argmin( np.sqrt( (xi-lon[i])2 + (yi-lat[i])2) ) zi[np.unravel_index(idx, xi.shape)] = temp[i]
--- Check... => LOOKS GOOD
plt.figure(figsize=(15,7)) plt.subplot(1,2,1) plt.scatter(lon, lat, temp, temp) plt.subplot(1,2,2) plt.pcolor(xi, yi, np.where(np.isnan(zi),0,zi)) plt.show()
--- Open NetCDF file to write on
with nc4.Dataset('test.nc', 'w' , format='NETCDF3_CLASSIC') as ds: # --- Initialize the dimensions of the dataset dim_time = ds.createDimension('time', 0) dim_lat = ds.createDimension('lat', yi.shape[0]) dim_lon = ds.createDimension('lon', xi.shape1)
# --- Create the corresponding variables for the dimensions
time = ds.createVariable('time', np.float32, 'time')
latitude = ds.createVariable('lat', np.float32, 'lat')
latitude.units = ['degrees north']
latitude.axis = ['Y']
latitude.standard_name = ['latitude']
longitude = ds.createVariable('lon', np.float32, 'lon')
longitude.units = ['degrees east']
longitude.axis = ['X']
longitude.standard_name = ['longitude']
# --- Fill with 1D (! 509