2

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?

The data needed to reproduce this can be downloaded here. Total download size = 10 KB.

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

# Read in data
toad = pd.read_csv('stack_question_csv.csv')
bnd = gpd.read_file("stack_boundary_final.shp")

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

# Convert DataFrame to GeoDataFrame
toad = gpd.GeoDataFrame(toad)
toad.head()

# Define coordinate systems
wgs84 = 4326 # WGS84

# Convert latitude and longitude to point data
toad['geometry'] = gpd.points_from_xy(toad['LON'], 
                                     toad['LAT'],
                                    crs = wgs84)

toad.head()

# 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')
toad.plot(ax=ax, marker='x', facecolor='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
gdf_grid.head()

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

gdf_grid.head()

# 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
x = toad['geometry'].x
y = toad['geometry'].y
z = toad['toads']

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

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

toad_nn = griddata(points,
                    values,
                    xi,
                    method='nearest')

#Save interpoalted data values into the geodataframe
gdf_grid['toads'] = toad_nn.round(1)

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

# Convert to UTM (EPGS:3721)
nad83 = 3721
gdf_grid['geometry'] = gdf_grid['geometry'].to_crs(crs=nad83)
gdf_grid['centroids'] = gdf_grid['centroids'].to_crs(crs=nad83)

gdf_grid.head()

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?

0

1 Answer 1

2
+100

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.

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.