# Save scipy interpolations as a raster (GeoTIFF) format using Python

I am using `scipy.interpolate` (Python 3.9) to interpolate some data. However, after using a nearest neighbor interpolator, I need to save it as a raster (GeoTIFF) image, preserving the resolution of 0.5m x 0.5m. How can I do this?

Here is the code I am using to generate my interpolations:

``````# Import modules
import geopandas as gpd
import matplotlib.pyplot as plt
plt.rcParams['axes.facecolor'] = 'black'
import pandas as pd
import numpy as np
from shapely.geometry import box
from scipy.interpolate import griddata

# Reproject shapefile boundary
bnd  = bnd.to_crs("epsg:4326")

# Convert DataFrame to GeoDataFrame

# Define coordinate systems
wgs84 = 4326 # WGS84

# Convert latitude and longitude to point data
crs = wgs84)

# Plotting of map boundary and sampled points
fig, ax = plt.subplots(1,1, figsize=(16,16))
ax.ticklabel_format(useOffset=False)
bnd.plot(ax=ax, facecolor='w', edgecolor='k')
plt.show()

# Create grid
epsg_utm = 32614 # UTM Zone 14
xmin, ymin, xmax, ymax = bnd.to_crs(epsg=epsg_utm).total_bounds #Represents field in terms of meters

# Define cell size
xdelta = 0.5 # meters
ydelta = 0.5 # meters

# Create an empty array to save the grid
grid = np.array([])

for x in np.arange(xmin, xmax, xdelta): #min, max step
for y in np.arange(ymin, ymax, ydelta): #min, max step
cell = box(x,y, x+xdelta, y+ydelta)
grid = np.append(grid, cell)

gdf_grid = gpd.GeoDataFrame(grid, columns=['geometry'], crs=epsg_utm)
gdf_grid['centroids'] = gdf_grid['geometry'].centroid

# Convert CRS back to Lat/Long
gdf_grid['geometry'] = gdf_grid['geometry'].to_crs(crs=wgs84)
gdf_grid['centroids'] = gdf_grid['centroids'].to_crs(crs=wgs84)

# Clip cells
gdf_grid = gpd.clip(gdf_grid, bnd['geometry'].iloc[0])
gdf_grid.reset_index(inplace=True, drop=True)

# Interpolate the values of observed weed density values to each centroid

xq = gdf_grid['centroids'].x
yq = gdf_grid['centroids'].y

# Arrange variables in griddata input format
points = (x,y)
values = z
xi = (xq, yq)

values,
xi,
method='nearest')

#Save interpoalted data values into the geodataframe

# Plot
fig, ax = plt.subplots(1,1, figsize=(10,11))
ax.ticklabel_format(useOffset=False)
gdf_grid.plot(ax=ax,
edgecolor='none',
cmap='RdYlGn',
antialiased=False)
plt.show()

# Convert to UTM (EPGS:3721)

print(gdf_grid.crs)
``````

After running this code, I can export the data in WKT with `gdf_grid.to_csv('data_export.csv')`, but I wish to export the data as a raster layer. How can I do this?

TLDR

``````#Your code
...
plt.show()

#Skip the last 6 lines of your code.

import rasterio
from rasterio import features
import affine

width = int((xmax-xmin)/xdelta)
height = int((ymax-ymin)/ydelta)

transform = rasterio.transform.from_bounds(xmin,
ymin,
xmax,
ymax,
width=width,
height=height)

gdf_grid['centroids'] = gdf_grid['centroids'].to_crs(crs=epsg_utm)
gdf2 = gdf_grid.drop('geometry', axis=1).values.tolist()

result = features.rasterize(
gdf2,
out_shape=(height, width),
transform=transform)

with rasterio.open(
'rasterized.tif', 'w',
driver='GTiff',
dtype=rasterio.uint8,
count=1,
width=width,
height=height,
crs=gdf_grid.crs,
transform = transform) as dst:
dst.write(result, indexes=1)
``````

The "export the data as a raster layer" part is not possible with a pure `geopandas` solution. I propose a solution using the library `rasterio` (another solution with `GDAL` is also possible).

In the code provided, the frequent conversions from one projection system to another using `geopandas.to_crs` function does not get reflected in the `dataframe.crs` property (see pull request). Therefore, the last projection change can be avoided (see my code comment "Skip the last 6 lines of your code"). The remaining code is just preparing data for `rasterio's - rasterize and write` functions.

Tip: Try using `GDAL's GRID` function to convert from point layer to raster layer.