<p>I changed your model solution a little bit; but it's like Ingvar Lukas wrote in his answer: you redefined <code>xi</code> and <code>yi</code> in the process, so when you later on define the netCDF values for <code>lat</code> and <code>lon</code> you try to fill small 1D arrays with a 2D array. That is the source of your error.</p>
<pre><code>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... =&gt; 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
</code></pre>
<p>That yields a georeferenced picture when loaded into QGIS (I reduced the step size of your <code>xi</code> and <code>yi</code> values to make the <code>pcolor</code> step a bit more performant; but you can change it back in your code).</p>
<p><a href="https://i.sstatic.net/iWJAB.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/iWJAB.png" alt="enter image description here" /></a></p>
<p><a href="https://i.sstatic.net/1MyVt.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/1MyVt.png" alt="enter image description here" /></a></p>

CC BY-SA 4.0
2021-10-28T10:12:02.407Z
thanks for your answer Erik! I commented Ingvar´s answer as well, can you refer to it? Its probably just a transformation issue..
Creative crypter
CC BY-SA 4.0
2021-10-28T11:01:13.477Z
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
CC BY-SA 4.0
2021-10-28T11:03:28.977Z
Hey, yea but i need to do it via a netcdf file generated in python.. is there no way to prepare the data in python so i can write it to netcdf? I thought i should first create a raster from the scattered data(latitude/longitude) and then write it to netcdf... When you look at the plot being generated, it actually looks correct..
Creative crypter
CC BY-SA 4.0
2021-10-28T11:16:45.183Z
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
CC BY-SA 4.0
2021-10-28T11:32:37.917Z
I modified the answer to do what I described above.
Erik
CC BY-SA 4.0
2021-10-28T12:23:22.887Z
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?
Creative crypter
CC BY-SA 4.0
2021-10-28T15:04:03.920Z
Let us [continue this discussion in chat](https://chat.stackexchange.com/rooms/130941/discussion-between-creative-crypter-and-erik).
Creative crypter
CC BY-SA 4.0
2021-10-28T09:37:53.643Z
2021-10-28T15:40:46.260Z
2021-10-28T15:40:46.260Z
Erik
<p>You have a shape mismatch, as you are overwriting<code>xi</code> and <code>yi</code> using <code>np.meshgrid</code>, then assigning wrong dimensions with <code>dim_lat</code> and <code>dim_lon</code> and eventually trying to fill <code>latitude</code> and <code>longitude</code> with your initial values of length 6.</p>
<p>Try modifying the grid preparation and interpolation</p>
<pre><code># 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')
</code></pre>
<p>.. and writing the values</p>
<pre><code># fill with values
latitude[:] = yi
longitude[:] = xi
value[0,:,:] = zi
value.grid_mapping = 'WGS84'
</code></pre>

CC BY-SA 4.0
2021-10-28T10:08:06.750Z
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?
Creative crypter
CC BY-SA 4.0
2021-10-28T11:35:51.050Z
@Creativecrypter Your comment is another question and should be posted as such, but maybe the GIS Stack Exchange would be more appropriate.
Jean-Marie Prival
CC BY-SA 4.0
2021-10-28T15:53:29.150Z
@Creativecrypter Link to the gis stack exchange https://gis.stackexchange.com/
Sam Dean
CC BY-SA 4.0
2021-10-28T18:03:37.637Z
Alright i will work on the codebase and if anything fails i will go there
Creative crypter
CC BY-SA 4.0
2021-10-28T08:32:46.410Z
2021-10-28T08:32:46.410Z
Ingvar Lukas
<p>Currently i am dealing with a bunch of coordinates (epsg:4326 =&gt; 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)</p>
<p>Currently the data is in scattered format, thats why its first interpolated.</p>
<p>After interpolation i tried to write the data (raster format) into a netcdf file, but thats failing:</p>
<pre><code>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 = &quot;&quot;&quot;GEOGCS[&quot;WGS 84&quot;,DATUM[&quot;WGS_1984&quot;,SPHEROID[&quot;WGS 84&quot;,6378137,298.257223563,AUTHORITY[&quot;EPSG&quot;,&quot;7030&quot;]],AUTHORITY[&quot;EPSG&quot;,&quot;6326&quot;]],PRIMEM[&quot;Greenwich&quot;,0,AUTHORITY[&quot;EPSG&quot;,&quot;8901&quot;]],UNIT[&quot;degree&quot;,0.01745329251994328,AUTHORITY[&quot;EPSG&quot;,&quot;9122&quot;]],AUTHORITY[&quot;EPSG&quot;,&quot;4326&quot;]]&quot;&quot;&quot;

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 # =&gt; ERROR: ValueError: shape mismatch: objects cannot be broadcast to a single shape
longitude[:] = lon
value[0,:,:] = zi
value.grid_mapping = 'WGS84'

ds.close()
</code></pre>
<p>When i try to write to e.g. latitude variable it throws:</p>
<pre><code>ValueError: shape mismatch: objects cannot be broadcast to a single shape
</code></pre>
<p>I debugged throw the code but wasnt able to find the error, any experts around?</p>
<p>Cheers!</p>
<p><strong>UPDATE #1:</strong></p>
<p>I adjusted your source code to work with coordinates in epsg:3857 format (it is plotted correctly in matplotlib but invisible in QGIS):</p>
<pre><code>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... =&gt; 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 = &quot;&quot;&quot;GEOGCS[&quot;WGS 84&quot;,DATUM[&quot;WGS_1984&quot;,SPHEROID[&quot;WGS 84&quot;,6378137,298.257223563,AUTHORITY[&quot;EPSG&quot;,&quot;7030&quot;]],AUTHORITY[&quot;EPSG&quot;,&quot;6326&quot;]],PRIMEM[&quot;Greenwich&quot;,0,AUTHORITY[&quot;EPSG&quot;,&quot;8901&quot;]],UNIT[&quot;degree&quot;,0.01745329251994328,AUTHORITY[&quot;EPSG&quot;,&quot;9122&quot;]],AUTHORITY[&quot;EPSG&quot;,&quot;4326&quot;]]&quot;&quot;&quot;    

    # --- 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
</code></pre>
<p>Do i have to make any changes to netcdf variables in netcdf? (like units or axis?)</p>
<p>Can i keep crs.spatial_ref / value.grid_mapping / value.grid_mapping_name ?</p>
<p>Cheers!</p>

CC BY-SA 4.0
2021-10-28T15:43:47.627Z
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
CC BY-SA 4.0
2021-11-23T00:07:03.683Z
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.. ?
Creative crypter
<geology><temperature><python><netcdf><interpolation>
Python - How to write interpolated raster to netcdf file
Creative crypter
# Python - How to write interpolated raster to netcdf file

(Asked by: username_0 on 2021-10-27)

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!

**Tags:** geology, temperature, python, netcdf, interpolation

> **Comments:**
> **username_1**: 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.  
> **username_0**: Hey @username_1 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.. ?  

---

## 2 answers

### Answer by username_1 on 2021-10-28. Score: 7

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

[![enter image description here](https://i.sstatic.net/iWJAB.png)](https://i.sstatic.net/iWJAB.png)

[![enter image description here](https://i.sstatic.net/1MyVt.png)](https://i.sstatic.net/1MyVt.png)

> **Comments:**
> **username_0**: thanks for your answer username_1! I commented username_2´s answer as well, can you refer to it? Its probably just a transformation issue..  
> **username_1**: 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?  
> **username_0**: Hey, yea but i need to do it via a netcdf file generated in python.. is there no way to prepare the data in python so i can write it to netcdf? I thought i should first create a raster from the scattered data(latitude/longitude) and then write it to netcdf... When you look at the plot being generated, it actually looks correct..  
> **username_1**: 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`...  
> **username_1**: I modified the answer to do what I described above.  
> **username_0**: Hi @username_1 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?  
> **username_0**: Let us [continue this discussion in chat](https://chat.stackexchange.com/rooms/130941/discussion-between-username_0-username_0-and-username_1).  

---

### Answer by username_2 on 2021-10-28. Score: 6

You have a shape mismatch, as you are overwriting`xi` 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'

```

> **Comments:**
> **username_0**: Hey username_2, 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?  
> **username_3**: @username_0 Your comment is another question and should be posted as such, but maybe the GIS Stack Exchange would be more appropriate.  
> **username_4**: @username_0 Link to the gis stack exchange https://gis.stackexchange.com/  
> **username_0**: Alright i will work on the codebase and if anything fails i will go there  

---

