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
gdal.Translate
or if irregular, withgdal.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.gdal.Translate
would avoid this whole part where I need to pivot the table? But thanks, I'll look into that.