Python - How to write interpolated raster to netcdf file

Question

Currently i am dealing with a bunch of coordinates (epsg:4326 => lat/lon) and their corresponding value (temperature). The goal is to write these coordinates and their values into a simple netcdf file and display it in e.g. QGIS (So you will have a colored square for each pixel/coordinate)

Currently the data is in scattered format, thats why its first interpolated.

After interpolation i tried to write the data (raster format) into a netcdf file, but thats failing:

import numpy as np
from scipy.interpolate import griddata
import xarray as xr
import pandas as pd
import rioxarray
import netCDF4 as nc4


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]

# prepare a grid for interpolation
xi = np.arange(6.0, 14.0, 0.001)
yi = np.arange(48.0, 64.0, 0.001)
xi, yi = np.meshgrid(xi, yi)
# As you can see the grid is slightly bigger then the used coordinates

# interpolate
zi = griddata((lon, lat), temp, (xi, yi), method='linear')

# time to write this into a netcdf file
ds = nc4.Dataset('test.nc', 'w', format='NETCDF3_CLASSIC')

dim_time = ds.createDimension('time', 0)
dim_lat = ds.createDimension('lat', len(yi))
dim_lon = ds.createDimension('lon', len(xi))

# is this correct or how should i set CRS to epsg:4326 ?
crs = ds.createVariable('WGS84', 'c')
crs.spatial_ref = """GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]"""

time = ds.createVariable('time', 'f4', 'time')

latitude = ds.createVariable('lat', 'f4', 'lat')
latitude.units = 'degrees north'

longitude = ds.createVariable('lon', 'f4', 'lon')
longitude.units = 'degrees east'

value = ds.createVariable('temperature', 'f4', ('time', 'lat', 'lon',))

# fill with values
latitude[:] = lat # => ERROR: ValueError: shape mismatch: objects cannot be broadcast to a single shape
longitude[:] = lon
value[0,:,:] = zi
value.grid_mapping = 'WGS84'

ds.close()

When i try to write to e.g. latitude variable it throws:

ValueError: shape mismatch: objects cannot be broadcast to a single shape

I debugged throw the code but wasnt able to find the error, any experts around?

Cheers!

UPDATE #1:

I adjusted your source code to work with coordinates in epsg:3857 format (it is plotted correctly in matplotlib but invisible in QGIS):

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

# --- 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.arange(min(lon), max(lon), 1000)
yi = np.arange(min(lat), max(lat), 1000)
xi, yi = np.meshgrid(xi, yi)
zi = np.zeros_like(xi,dtype=np.float32) * -999
# 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]

# Replace -0 with nan value (so i have invisible pixels instead of black background)
np.place(zi, zi == -0, None)

# --- 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.shape[1])

    # --- 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 (!) arrays of xi/yi, as the meshgrid returns 2D arrays...
    time[:] = 0
    latitude[:] = yi[:,0]
    longitude[:] = xi[0,:]
    
    # --- Create a coordinate reference system
    crs = ds.createVariable('WGS84', 'c')
    crs.spatial_ref = """GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]"""    

    # --- Ready the Temperature data field
    value = ds.createVariable('temperature', np.float32, ('time','lat','lon'))
    
    value.grid_mapping = 'WGS84' # the crs variable name
    value.grid_mapping_name = 'latitude_longitude'
    
    # --- Fill with values
    value[0,:,:] = zi

Do i have to make any changes to netcdf variables in netcdf? (like units or axis?)

Can i keep crs.spatial_ref / value.grid_mapping / value.grid_mapping_name ?

Cheers!

Answer 1

I changed your model solution a little bit; but it's like Ingvar Lukas 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.

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.shape[1])

    # --- 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 (!) arrays of xi/yi, as the meshgrid returns 2D arrays...
    time[:] = 0
    latitude[:] = yi[:,0]
    longitude[:] = xi[0,:]
    
    # --- Create a coordinate reference system
    crs = ds.createVariable('CRS', 'c')
    crs.spatial_ref = CRS_wkt

    # --- Ready the Temperature data field
    value = ds.createVariable('temperature', np.float32, ('time','lat','lon'))
    
    value.grid_mapping = 'CRS' # the crs variable name
    value.grid_mapping_name = 'latitude_longitude'
    
    # --- Fill with values
    value[0,:,:] = zi

That yields a georeferenced picture when loaded into QGIS (I reduced the step size of your xi and yi values to make the pcolor step a bit more performant; but you can change it back in your code).

Answer 2

You have a shape mismatch, as you are overwritingxi and yi using np.meshgrid, then assigning wrong dimensions with dim_lat and dim_lon and eventually trying to fill latitude and longitude with your initial values of length 6.

Try modifying the grid preparation and interpolation

# prepare a grid for interpolation
xi = np.arange(6.0, 14.0, 0.001)
yi = np.arange(48.0, 64.0, 0.001)
xi_mesh, yi_mesh = np.meshgrid(xi, yi)

# interpolate
zi = griddata((lon, lat), temp, (xi_mesh, yi_mesh), method='linear')

.. and writing the values

# fill with values
latitude[:] = yi
longitude[:] = xi
value[0,:,:] = zi
value.grid_mapping = 'WGS84'