# Python - How to write interpolated raster to netcdf file

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!

• I have modified my answer to account for the change in coordinate reference system. Probably it's good to stop altering your question multiple times, and to simply ask a new question. For example, of an example of how to create a netCDF file that contains data on a "completely irregular array". That should be able to accomodate your data without having to go through the griddata etc. steps, to only supply the data you have, and nothing more.
– Erik
Oct 28 at 15:43
• Hey @Erik is there a name or literature for this algorithm? like buffer/proximity analysis, euklid distance, euklid allocation? because it seems like a mix of all these methods.. ? Nov 23 at 0:07

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

• thanks for your answer Erik! I commented Ingvar´s answer as well, can you refer to it? Its probably just a transformation issue.. Oct 28 at 10:12
• As you use scipy.interpolate.griddata, you get a linear interpolation of the data onto the grid; that's where the problem resides... not on the netCDF side. I don't have knowledge of whether netCDF and QGIS support unstructured (i.e., non-regular lat/long) grid data. but it sounds like that might be what you really want. Or you can, e.g., simply make a .CSV file with column1=lon, colum2=lat, column3=value, and load it as a set of points into QGIS? Wouldn't that be a better solution?
– Erik
Oct 28 at 11:01
• Well, within the constraints of a raster file, then I'd suggest to set zi[:,:]=0 and loop over your data coordinates and fill zi only for those xi, yi coordinates that are minimal at the corresponding lon and lat coordinates. Assuming rectangular grids preserve Pythagorean distance, idx=np.argmin(np.sqrt( (xi-lon[i])**2 + (yi-lat[i])**2)) should give you the location of where to set zi[idx]=temp[i], where you should just loop over i...
– Erik
Oct 28 at 11:16
• I modified the answer to do what I described above.
– Erik
Oct 28 at 11:32
• Hi @Erik thank you so much, i start to understand how to work in this environment. I adjusted your source code with one little edit: I transformed the lat/lon coordinates into epsg:3857 format. The matplotlib output looks good, but when i place the file in QGIS there is no result. Can you check my updated source code? Oct 28 at 12:23

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'

• Hey Ingvar, thanks for your answer, i edited the source code and there is no more error. When i drag and drop the netcdf file into QGIS however is draws a polygon... but i just wanted to display six colored squares on the coordinates e.g. [50.1, 8.1, white(1)], [50.2, 8.2, lightgray(2)], ... , [52, 12, black(6)], is there any way to achieve this? Oct 28 at 10:08
• @Creativecrypter Your comment is another question and should be posted as such, but maybe the GIS Stack Exchange would be more appropriate. Oct 28 at 11:35
• @Creativecrypter Link to the gis stack exchange gis.stackexchange.com Oct 28 at 15:53
• Alright i will work on the codebase and if anything fails i will go there Oct 28 at 18:03