1

I'm trying to understand why my multiband raster output created with GDAL is distorted.

Short summary: I got several hundreds of US counties data in csv format. For each one, I need to convert the CSV into a multiband raster. Each csv file have 24 columns, where the first two columns are X and Y UTM coordinates. For most of the counties, the script is running without issues. However, some counties have some barrier inlands which I suspects distorts margins of the raster.

I will provide the reproducible example along with the code I am using. I know is pretty a long code and I apologies for the size of the shapefile (that was the only I could find that had the same boundary with my raster), but is already been like several days since I trying to debug this. I tried rasterio as well and I got the same issue. Only geocube library was working properly but I would prefer not to use it.

import numpy as np
import geopandas as gpd
import random
import pandas as pd
from osgeo import gdal, osr
from urllib.request import urlretrieve
from zipfile import ZipFile 

######## TO REPRODUCE THE ISSUE, I WILL CREATE A GRID WITH POINTS AT 90 M EACH ############
# Download file from ESRI
url = (
       'https://www.arcgis.com/sharing/rest/content/items/715db3ed501b42fe9581caaa5c56caf9/data'
       )
filename = 'countyshapefiles.zip'
urlretrieve(url, filename) # the zip has a size of ~80 mb

# Unzip
with ZipFile('countyshapefiles.zip', 'r') as z:
    z.extractall(
        path = 'countyshapefiles'
    )

# Read the file
shapefile = r"countyshapefiles\USA_Counties.shp"
epsg = 32617

df = gpd.read_file(shapefile)
df = df.to_crs(f'EPSG:{epsg}')

df = df[(df['NAME'] == 'Hillsborough') & (df['STATE_NAME'] == 'Florida')] # this is the tricky county
df.plot()

# Extract the boundaries
xmin, ymin, xmax, ymax = df.total_bounds

# Set distance between points
spacing = 90 # in reality my data is at 30 m but for time/space saving purpose, I'll set it to 90.

# Get the coordinates
xcoords = [i for i in np.arange(xmin, xmax, spacing)]
ycoords = [i for i in np.arange(ymin, ymax, spacing)]

pointcoords = np.array(np.meshgrid(xcoords, ycoords)).T.reshape(-1, 2) #A 2D array like [[x1,y1], [x1,y2], ...
points = gpd.points_from_xy(x=pointcoords[:,0], y=pointcoords[:,1])
grid = gpd.GeoSeries(points, crs=df.crs)
grid.name = 'geometry'

# find points falling inside the poly
gridinside = gpd.sjoin(gpd.GeoDataFrame(grid), df[['geometry']], how="inner")

# Generate some random data
gridinside['variable1'] = gridinside.apply(lambda _: random.randint(0, len(gridinside.index)), axis=1)
gridinside['variable2'] = gridinside.apply(lambda _: random.randint(0, len(gridinside.index)), axis=1)


# Extract the X, Y coords
gdf = gridinside.get_coordinates(ignore_index=True)
gdf['var1'] = gdf.apply(lambda _: random.randint(0, len(gdf.index)), axis=1)
gdf['var2'] = gdf.apply(lambda _: random.randint(0, len(gdf.index)), axis=1)


######## LETS'S CREATE THE RASTER NOW ############
# Transform to dataframe
df = pd.DataFrame(gdf)

# Pivot the table 
dfP = df.pivot(index = 'y', columns='x', values=df.iloc[:, 2:4].columns)

# Get the width/height of the raster
height = df['y'].nunique()
width = df['x'].nunique()

# Create the bounds of the raster 
xmin, ymin, xmax, ymax = min(df['x']), min(df['y']), max(df['x']), max(df['y'])

# Specify the cell size
res_x = spacing
res_y = -spacing

# Output raster name
output_tif = '12057_fake_90m.tif'

# Create the dataset
driver = gdal.GetDriverByName('GTiff')
dataset = driver.Create(output_tif, width, height, 2, gdal.GDT_Int16) # 2 bands
dataset.SetGeoTransform((xmin, res_x, 0, ymax, 0, res_y))

srs = osr.SpatialReference()
srs.ImportFromEPSG(epsg)  # Set the projection
dataset.SetProjection(srs.ExportToWkt())

columns = df.iloc[:, 2: 4].columns
columns

# Write the data to the UTM raster bands
for i, column in enumerate(columns, start=1):
    band = dataset.GetRasterBand(i)
    column_data = np.array(dfP[column][::-1])
    column_data[np.isnan(column_data)] = -9999
    band.WriteArray(column_data)
    band.SetNoDataValue(-9999)
    band.SetDescription(column)
dataset = None
2
  • 1
    Rather than reinventing the wheel, why not just use basic gdal functionality? If the CSV XY values form a regular grid use gdal.Translate or if irregular, with gdal.Grid. You may need to extract out just the X, Y and specific "Z" column from the CSV for each band. Then clip to shapefile.
    – user2856
    Commented Nov 2, 2023 at 0:53
  • so gdal.Translate would avoid this whole part where I need to pivot the table? But thanks, I'll look into that. Commented Nov 2, 2023 at 9:06

1 Answer 1

2

Assuming the CSV XY values form a regular grid, you can treat the CSV as an XYZ raster, though you will need to extract out just the X, Y and specific "Z" column from the CSV for each band that will form the multiband output. You can then clip to your shapefile with gdal.Warp and stack the clipped results.

import pandas as pd
from osgeo import gdal
from urllib.request import urlretrieve
from zipfile import ZipFile

# # Download file from ESRI
# url = (
#        'https://www.arcgis.com/sharing/rest/content/items/715db3ed501b42fe9581caaa5c56caf9/data'
#        )
# filename = 'countyshapefiles.zip'
# urlretrieve(url, filename) # the zip has a size of ~80 mb
#
# # Unzip
# with ZipFile('countyshapefiles.zip', 'r') as z:
#     z.extractall(
#         path = 'countyshapefiles'
#     )

shapefile = r"countyshapefiles/USA_Counties.shp"
epsg = 32617
csvfile = "hill_multicolumn.csv"  # Assumes X, Y, col1,...,colN

df = pd.read_csv(csvfile)
df = df.sort_values(["Y", "X"], ascending=True)
vrts = []
for col in df.columns[2:]:
    df_col = df.filter(["X", "Y", col])
    df_col.to_csv(f"hill_{col}.csv", index=False)
    rds = gdal.OpenEx(f"hill_{col}.csv", gdal.OF_RASTER)
    vrt = gdal.Warp(
        f"/vsimem/hill_{col}.vrt",
        rds,
        cutlineDSName=shapefile,
        cutlineWhere="NAME = 'Hillsborough' AND STATE_NAME = 'Florida'",
        cropToCutline=True,
        srcSRS=f"EPSG:{epsg}"
    )
    vrts.append(vrt)

vrt = gdal.BuildVRT(f"/vsimem/hill.vrt", vrts, separate=True)
multiband = gdal.Translate("hill_multiband.tif", vrt)

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.