Store raster output in GeoPackage, along with vector data
Closed this issue · 3 comments
In GitLab by @Huite on Sep 27, 2022, 16:51
rioxarray should be able to write a dataset, via rasterio, as it supports subdatasets.
with rasterio.open(path, "r+", raster_table="new_table", append_subdataset=True) as destination:
destination.write(array)[I think rioxarray will do this, provided we add the additional raster_table and append_subdataset kwargs in the to_raster call.
This might be a better solution that using the netCDF as currently, if we can make a stack of rasters in the geopackage into a Temporal raster entity.
In GitLab by @Huite on Sep 27, 2022, 21:39
This is indeed easy to setup:
import numpy as np
import rioxarray
import xarray as xr
import imod
import geopandas as gpd
# %%
like = imod.util.empty_2d(
dx=100.0, xmin=0.0, xmax=20_000.0, dy=-100.0, ymin=300_000.0, ymax=320_000.0
)
da = xr.full_like(like, 5.0)
# %%
xy = np.array([
[500.0, 300_500.0],
[1500.0, 301_500.0],
[2500.0, 302_500.0],
])
gdf = gpd.GeoDataFrame(geometry=gpd.points_from_xy(xy[:, 0], xy[:, 1]))
vector_outpath = "results0.gpkg"
gdf = gdf.set_crs(epsg=28992)
gdf.to_file(vector_outpath, driver="GPKG", layer="observations")
# %%
da.rio.write_crs(28992, inplace=True)
da.astype(np.float32).rio.to_raster(
vector_outpath,
driver="GPKG",
raster_table="head_layer0",
append_subdataset=True,
)In GitLab by @Huite on Sep 27, 2022, 21:42
Having time dependent rasters still requires loading all of the rasters into the Layers panel, then setting the temporal properties on them. This is not ideal, though obviously manageable. netCDF raster support will not do the trick either, and meshes have too little features available.
For steady-state, it seems like storing the rasters in a geopackage is almost definitely a win.