1

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

1 Answer 1

2

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)
2
  • Thanks! I tried this with my real example. For some reason, the output raster is being "stretched" now to the left. I tried multiple transformations and each returned the same. Any idea what could be the cause? Commented Sep 28, 2023 at 13:25
  • 1
    The transform is not the issue. You say your extent is 4km x 4km but it's not, Check them, you'll see xmax-xmin is 10.6km and ymax-ymin is 5.8km.
    – user2856
    Commented Sep 28, 2023 at 21:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.