Create multiband raster in Python from X, Y

Question

How do you create a multiband raster in Python and save it as a .tif file from a long dataframe format?

import pandas as pd
import numpy as np
from rasterio.crs import CRS
import rasterio
from rasterio.transform import from_origin

# Define the extent of the raster (4 km on each side)
xmin, ymin, xmax, ymax = (325631.8, 3051603, 336281.8, 3057483)  # in meters

# Define the resolution (30 meters)
res_x = 30
res_y = -30

# Create a grid of points within the specified extent and resolution
x_values = np.arange(xmin, xmax, 30)
y_values = np.arange(ymin, ymax, 30)

# Create an empty list to store the data
data = []

# Iterate over the grid and populate the list with random values
for x in x_values:
    for y in y_values:
        # Random value
        value1 = np.random.rand()
        value2 = np.random.rand()
        
        # Append the data as a dictionary
        data.append({'x': x, 'y': y, 'value1': value1, 'value2': value2})

# Create the GeoDataFrame from the list of dictionaries
df = pd.DataFrame(data)

What I tried so far (with rasterio library) was to pivot the table first, then do some transformation according to the documentation (although I can't exactly tell which one to use and when). But I suppose from_origin() is appropriate to use here:

df = df.pivot(index = 'y', columns = 'x', values=columns[2:]) # pivot table
transform = from_origin(xmin, ymin, res_x, res_y) #  create a tranformation
crs = CRS.from_epsg(32617)
height = df['x'].nunique()
width = df['y'].nunique()
columns = ['x', 'y', 'value1', 'value2']

And then something like:

# Create a multiband raster from the DataFrame
with rasterio.open('myrast.tif', 'w', driver='GTiff', height=height, width=width, count=len(columns) - 2, dtype='float32', crs=crs, transform=transform) as dst:
    for i, column in enumerate(columns[2:], start=1):  # Skip the first two columns (x and y)
        dst.write(allDataPivot[column].values.astype('float32'), i)
# Close the raster file
dst.close()

This is working but I found that in some cases, I get a "shifted" raster to the East (like with several kilometers) which raises me some questions about the method I used.

But it's that such complicated indeed? Worst thing is that I need to wrap all this in a PyInstaller where I get other issues like bundling the rasterio library with all its modules inside the .exe (so fewer imports, better the chances).

As a comparison to R, that would be such a simple job:

library(terra)
myrast = rast(df) # that's all you need to do.
writeRaster(myrast, 'myraster.tif', filetype = 'GTiff') # I think it was just like that, doing it from memories

Answer 1

You're doing a lot of unnecessary looping. You can create your random values in a single array and you don't need a dataframe.

import numpy as np
import rasterio as rio
from rasterio.transform import from_bounds

dst_filename = 'raster.tif'

# Define the extent of the raster
xmin, ymin, xmax, ymax = (325631.8, 3051603, 336281.8, 3057483)  # in meters

# Define the resolution (30 meters)
res_x = 30
res_y = -30

#Raster size
height = int((ymax - ymin) / res_y)
width = int((xmax - xmin) / -res_x)
count = 2

dtype = np.float32

data = np.random.random((count, height, width)).astype(dtype)
crs = "EPSG:32617"

# You can create a transform in a couple of ways, either
transform = rio.Affine(res_x, 0, xmin, 0, res_y, ymax)
# or
transform = from_bounds(xmin, ymin, xmax, ymax, width, height)

# Now write out the data in one step    
with rio.open(
        dst_filename, 'w',
        driver='GTIFF',
        height=height,
        width=width,
        count=count,
        dtype=dtype,
        crs=crs,
        transform=transform
    ) as ds:
    ds.write(data)