nmarchio
commented
1 year ago Hi @martinfleis -- not sure if this helps (or if this is duplicative with functionality in this repo) but I developed these two functions to convert GeoTiff rasters to square vector geometries using xarray (so the flow is, GeoTiff > xarray > pd.DataFrame > gpd.GeoDataFrame > GeoParquet). Below is a reproducible example working with a population raster dataset. I'm sure the process could be improved and generalized further, but figured I'd share in case its useful.
Nico
# %%
import pygeos
import numpy as np
import pandas as pd
import geopandas as gpd
import rasterio as rio
from rasterio import features
import xarray as xr
import rioxarray as rxr
# %%
def vectorize_pixels(xr_array: xr.DataArray) -> xr.Dataset:
"""Vectorize pixel coordinates in a xarray.DataArray raster
Args:
xr_array: xarray.DataArray of target raster containing dimensions 'y' and 'x'
Returns:
xarray.Dataset with corner coordinates of each pixel as data variables
"""
bounds = xr_array.rio.bounds()
transform_constructor = rio.transform.from_bounds(west = bounds[0], south = bounds[1], east = bounds[2], north = bounds[3], width = len(xr_array['x']), height = len(xr_array['y']))
height = len(xr_array['y'])
width = len(xr_array['x'])
cols, rows = np.meshgrid(np.arange(width), np.arange(height))
xs, ys = rio.transform.xy(transform = transform_constructor, rows = rows, cols = cols, offset = 'center')
xs_ul, ys_ul = rio.transform.xy(transform = transform_constructor, rows = rows, cols = cols, offset = 'ul')
xs_ur, ys_ur = rio.transform.xy(transform = transform_constructor, rows = rows, cols = cols, offset = 'ur')
xs_ll, ys_ll = rio.transform.xy(transform = transform_constructor, rows = rows, cols = cols, offset = 'll')
xs_lr, ys_lr = rio.transform.xy(transform = transform_constructor, rows = rows, cols = cols, offset = 'lr')
lons = np.array(xs)
lats = np.array(ys)
lons_ul = np.array(xs_ul)
lats_ul = np.array(ys_ul)
lats_ll = np.array(ys_ll)
lons_ur = np.array(xs_ur)
dims_remove = [x for x in list(xr_array.dims) if x not in set(['x','y'])]
xr_data = xr_array.to_dataset(name='data').drop_dims(drop_dims = dims_remove)
xr_data['ymax'] = xr.DataArray(data=lats_ul, coords = xr_data.coords, dims=('y', 'x'))
xr_data['ymin'] = xr.DataArray(data=lats_ll, coords = xr_data.coords, dims=('y', 'x'))
xr_data['xmax'] = xr.DataArray(data=lons_ur, coords = xr_data.coords, dims=('y', 'x'))
xr_data['xmin'] = xr.DataArray(data=lons_ul, coords = xr_data.coords, dims=('y', 'x'))
return xr_data
def shapelify_pixels(pd_data: pd.DataFrame) -> gpd.GeoDataFrame:
"""Convert pandas.DataFrame containing bounding box columns to a geopandas.GeoDataFrame
Args:
pd_data: pandas.DataFrame returned from merge_pixels() function
Returns:
pandas.DataFrame converted to geopandas.GeoDataFrame with pixels formatted as geometries
"""
coords = [pygeos.box(xmin, ymin, xmax, ymax) for xmin, ymin, xmax, ymax in zip(pd_data['xmin'], pd_data['ymin'], pd_data['xmax'], pd_data['ymax'])]
data = gpd.GeoDataFrame(pd_data, geometry = pygeos.to_shapely(coords), crs ="EPSG:4326")
return data
# %%
# Download data
# wget -O /Users/nm/Downloads/worldpop/DJI.tif https://data.worldpop.org/GIS/Population/Global_2000_2020_Constrained/2020/maxar_v1/DJI/dji_ppp_2020_UNadj_constrained.tif
# Read in geotif data
xr_array = rxr.open_rasterio(filename = '/Users/nm/Downloads/worldpop/DJI.tif')
# %%
# Vectorize geotif
raster_geoms = vectorize_pixels(xr_array = xr_array)
xr_array.name = 'djipop'
raster = xr.merge(objects = [xr_array, raster_geoms], join= 'left')
# %%
# Convert xarray to df
df = raster.to_dataframe()
# Convert df to gdf
gdf = shapelify_pixels(pd_data = df)
# %%
# Plot as a gdf
gdf[gdf['djipop'] > 0].plot(column = 'djipop', cmap='viridis').reset_index()
rabernat
dcherian
e-kotov
edzer
martinfleis
benbovy
thomcom
paleolimbot
dblodgett-usgs
kylebarron
At the moment, we have a decent set of functionality to work with in-memory vector data cubes. What we do not have is a way of saving them to a disk and reading back again, apart from using
picklewhich is not very interoperable.We should find a way of saving to a file format that allows us interoperability with @edzer's
starsin RSpatial. @edzer, you were mentioning that NetCDF is able to handle VDC. Doesstarsimplements that? Can we create a data cube indexed by geometry and save it via GDAL to NetCDF? I tried only very briefly to play with thewrite_starsfunction to no avail[^1].The spec for NetCDF way of storing geometries is here but as someone who never worked with the file format, it is not super clear to me how it all works together.
Over in Python, there is some work @dcherian posted on conversion of shapely geometries to CF geometries and back (https://cf-xarray.readthedocs.io/en/latest/generated/cf_xarray.shapely_to_cf.html) [^2]. There is also https://github.com/twhiteaker/cfgeom which may potentially be used (although it is only a reference implementation that is not maintained).
So it seems that NetCDF should be able to store vector data cubes and we may just need to figure out some wrappers for a convenient IO from xvec-backed xarray objects. But I'd welcome some help or at least a guidance from someone more familiar with CF conventions.
The other option seems to be Zarr as discussed here but that is at this stage only an idea and (Geo)Zarr spec is not ready as of now.
The last option is to convert the cube to a long-form dataframe and save it as a GeoParquet but that kind of breaks the point of having a cube in the first place.
[^1]: Getting
Error in !all.equal(match(xydims, names(d)), 1:2) : invalid argument typewhen trying to usewrite_stars(st, 'test.nc')on the example cube from https://r-spatial.github.io/stars/articles/stars1.html#vector-data-cube-example. [^2]: only points are implemented so far but that may change