Open ctoney opened 5 months ago
This already covers a lot what vapour does with vector layers.
mvfile <- system.file(file.path("extdata/tab", file), package="vapour")
x <- new(gdalraster:::GDALVector, mvfile)
x$getDriverLongName() ## vapour_driver
## vapour_layer_info, vapour_report_fields
x$bbox()
x$getFeatureCount()
x$getSpatialRef()
x$getLayerDefn()
x$getName()
## we would want some helpers to avoid materializing geometry and attributes, so we can skip/limit, and get just bbox etc
## vapour_read_geometry(limit_n = 2, skip_n = 4)
## vapour_geom_summary()/vapour_read_extent()
x$getNextFeature(); x$getNextFeature();
## now iterate for the next 4 features to match the skip_n/limit_n of vapour
x$getGeometryColumn() ##vapour_geom_name()
I can't spend more time this week but I feel like this could easily replace the innards of vapour. I'll have a closer look and see if there's anything missing across the raster stuff as well (I'm pretty well established in what and how I use vapour now so it won't take me much to investigate the core pieces I really need). I would look to replacing my dependence on vapour with this, and probably just let vapour sit as-is rather than work more on it (though it might be fun to modify it to dep on this, and I will need to think about the one non-me revdep vapour has {ursa}).
I think ideally we would want to be able to spit out vectors of geos or ogr pointers, or wkb/wkt with support by {wk} (or even the arrow stream from GDAL itself) . There's certainly significant interest in arrow and geoarrow, so it might be easy to get contributions to help streamline this with the wk/geoarrow approach with this as the primary GDAL-API package for R (but I'll have to think about that some more at a later time.)
@mdsumner, thank you for the helpful comments. Look forward to working through more of the details.
I added “vectors of geos or ogr pointers, or wkb/wkt with support by {wk}” in the section "Further consideration / TBD".
Added GDALVector::setIgnoredFields()
in the class definition (not in prototype yet).
Added arguments fields
and geom_column
in GDALVector::getFeatureSet()
(not in prototype yet).
I had the OGR Arrow interface under "TBD" for GDALVector::getArrowStream()
(GDAL >= 3.6) and GDALVector::writeArrowBatch()
(GDAL >= 3.8), but I don't have that expertise currently. A contributor would be awesome if that were possible. The basic implementation of getFeatureSet()
would populate a data frame by iterating feature rows and appending across the column vectors. The ability to access the data by column in a memory efficient way seems like a major improvement for large datasets.
This is really cool to experiment with so far!
I am currently trying to clip a source layer by another one but hitting Error: could not find valid method.
(this would be a very common use case inside the mapme package).
Below is my code, does anyone have an idea what is going wrong?
library(sf)
library(gdalraster)
nc_org <- system.file("shape/nc.shp", package="sf")
nc_out <- "/vsimem/nc.gpkg"
ogr2ogr(nc_org, nc_out, cl_arg = c("-t_srs", "EPSG:4326", "-nlt", "PROMOTE_TO_MULTI"))
url <- "/vsicurl/https://data.source.coop/vida/google-microsoft-open-buildings/flageobuf/by_country/country_iso=USA/USA.fgb"
buildings_usa <- new(GDALVector, url)
nc <- new(GDALVector, nc_out)
srs <- buildings_usa$getSpatialRef()
dsn_out <- "/vsimem/nc_buildings.gpkg"
gdalraster:::.create_ogr("GPKG", dsn_out, 0, 0, 0, "Unknown", "nc_buildings", srs)
nc_buildings <- new(GDALVector, dsn_out, "nc_buildings", read_only = FALSE)
buildings_usa$layerClip(nc, nc_buildings)
buildings_usa$close()
nc$close()
nc_buildings$close()
the layerClip call is missing two required arguments:
atm you'll have to specify all args explicitly, usually we would do that in the R wrapper with default values for most args
(untested, just I'm accustomed to seeing that error message) 🙏
Ah, right! That seems to make it work. Thanks.
The layer geoprocessing functions were done as class methods in GDALVector
mainly because it mimics the underlying API (class methods of OGRLayer
in the C++ API). It was quick and easy to add them in GDALVector
with RCPP_EXPOSED_CLASS
, for prototype purposes. We might instead make them stand-alone functions in the R interface. That would look something like
layerIntersection(input_layer, method_layer, result_layer, quiet = FALSE, options = NULL)
in which case we need to supply the input_layer
too. That works fine since we can pass around objects of class GDALVector
which may have an attribute and/or spatial filter already defined. Also, we could optionally create the output layer. So result_layer
could be specified as currently a GDALVector
object, or alternatively with additional arguments for DSN/layer potentially to be created.
Added a note to that effect in the draft document for TBD.
Also added in the draft document, a link to the header file src/gdalvector.h
. That file gives a quick reference to the exposed class methods that have been implemented so far in the prototype (a subset of the draft class definition).
@goergen95 , I experimented with your example above. I noticed that the USA buildings .fgb file is pretty large at 36 GB with 145 million features. There are close to 9 million features that intersect the NC bbox. When I try to run the clip for all NC counties it processes for a long time without getting very far. I'm curious if you would expect that to work in reasonable time, i.e., is that a representative data size/processing task?
The FlatGeobuf file spec also says:
If you're accessing a FlatGeobuf file over HTTP, consider using a CDN to minimize latency. In particular, when using the spatial filter to get a subset of features, multiple requests will be made. Often round-trip latency, rather than throughput, is the limiting factor. A caching CDN can be especially helpful here.
So latency could be an issue in this case, or maybe Source Cooperative does the caching by default?
The code below filters the NC layer to one county. It also sets a spatial filter on the buildings layer first. The clip for Ashe county then runs pretty fast:
Hi @ctoney,
processing this amount of data is not representative of what we are dealing with (though there are some requesters interested in the building footprints data set). The most complex vector data source we are currently working with most probably are the IUCN species ranges. These have to be downloaded locally anyway, so latency is not the limiting factor here. I would not expect gdalraster
to handle the building footprint data set while being located on a remote server with ease. If one was really going to analyze it globally, copying the data to a server with lower latency would be the way to go, I guess.
With that in mind, I wanted to say that your vignette on GDAL's block cache is really helpful to optimize read operations on rasters. I wonder if it was worthwhile to work on something along the same lines for vector data? We could consider e.g. number of features, complexity of geometries, GDAL options, and file formats with regard to I/O performance?
@goergen95 , thanks for the additional information. I'm still experimenting with the building footprints. It's likely just throughput that's limiting for the ~9 million in NC. It sounds like Source Cooperative provides what the FlatGeobuf spec recommends: https://github.com/flatgeobuf/flatgeobuf#optimizing-remotely-hosted-flatgeobufs. It's good to have a couple of examples like this for testing. I'll take a look at the IUCN species ranges too.
Right, there are several format-specific options and other configuration with regard to I/O performance for vector processing. That's a good idea to start collecting those into one place for reference. I could try to do that and eventually put a draft somewhere that we can add to.
Hi,
I am currently experimenting a bit more with the vector API. I hit a wall when I was trying to write a vector data set from scratch. Below code assumes that I obtained a number of bounding boxes somehow (a real use-case might be the footprints from items in a STAC catalog). I seem not to be able to find any methods to write features to a newly created GDALVector
dataset. The examples in the vignette seem to operate on copies of datasets. Is it the case that we are currently not able to write features or am I missing something obvious?
Hi @goergen95,
That is correct, we are currently not able to write features, at least with the prototype of GDALVector
. The class definition in the specification document is intended to reflect proposed/planned functionality. Currently, the only way to know what has actually been implemented in the prototype is to scan through the header file in the gdalvector
branch: https://github.com/USDAForestService/gdalraster/blob/gdalvector/src/gdalvector.h. Under class GDALVector
, then under the public:
specifier at line 33, the signatures for the public constructors are given, and then signatures for the public methods that have been implemented so far. The prototype is basically read-only at this point, except that the OGR geoprocessing methods are available to write new layers (layerIntersection()
etc.)
Note that ogrinfo
could be used to edit data with SQL statements. The documentation for gdalraster::ogrinfo()
has a couple of examples. But in general this would be a little clunky for writing data compared with what we're envisioning for the class-based interface, having methods for (at least) createFeature()
, setFeature()
, upsertFeature()
, and deleteFeature()
, potentially used inside startTransaction()
/ commitTransaction()
.
Thanks for the nudge though. I expect to refocus on the vector API sometime soon and prioritize that for development over the coming weeks. It's a good idea to get some basic write capability into the prototype sooner than later, so that it can be included in early testing/feedback as you're attempting to do now. I'll post updates here once anything meaningful is added. In the meantime, please do send feedback or suggestions if anything else comes up.
@goergen95, Correction to the above
Currently, the only way to know what has actually been implemented in the prototype is to scan through the header file in the
gdalvector
branch
It may be more convenient to view that information in R. Calling GDALVector
will return the built-in documentation provided by Rcpp
. Once an object has been created from a RCPP_EXPOSED_CLASS
, str()
will return the list of methods that can be called on the object via the $
operator. Example output below:
Another minor correction to the above.
Using ogrinfo()
to edit data might be a little clunky only in the cases when we want to edit by operating on individual feature records (e.g., accessing the dataset via cursor, consistent with the use case that @goergen95 described above). In those cases, it is likely better to construct a list of field name=value for a feature, including a geom field containing WKT, and pass to an insert/update/upsert method on the layer, with FID of the affected feature as return value. (Of course, that functionality is one of the motivations to develop this interface in the first place, and I hope to have some write capability available for testing in the near future).
We will probably also have an execute_sql()
that operates on OGR data sources (wrapper of GDALDatasetExecuteSQL()
). We can already create a layer object from a SELECT statement in the constructor for GDALVector
, but execute_sql()
would be used for ALTER TABLE, CREATE INDEX, INSERT, UPDATE etc. that do not return a result set. This would be redundant with functionality that is now provided by gdalraster::ogrinfo()
, but could be more specialized. Also, note that ogrinfo
the command-line utility, was exposed via the GDAL utils API somewhat recently at GDAL 3.7, so has that version requirement for using via the wrapper in gdalraster
.
The clarification here is mainly that ogrinfo()
quite helpfully provides functionality to edit data with SQL statements, and it is not especially clunky for that usage in general. It would even be preferred for some cases (the example in https://usdaforestservice.github.io/gdalraster/reference/ogrinfo.html). It can be used for editing any writeable OGR vector data source, not just RDBMS-based. Also worth noting that we can use both the OGR SQL and SQLite dialects, the latter providing INSERT/UPDATE/DELETE:
https://gdal.org/user/ogr_sql_dialect.html https://gdal.org/user/sql_sqlite_dialect.html
(apologies to the extent that any of the above may be well-known already, but hopefully worth clarifying)
Utility functions for managing vector data sources have been added in main (#343), to be included in a 1.11 release in the next week or two hopefully. Documented at: ogr_manage - https://usdaforestservice.github.io/gdalraster/reference/ogr_manage.html and ogr_define - https://usdaforestservice.github.io/gdalraster/reference/ogr_define.html
These changes have also been merged into the gdalvector branch, and the examples in the draft document for the vector API have been updated: https://usdaforestservice.github.io/gdalraster/articles/gdalvector-draft.html
In case you install from the updated gdalvector branch, another notable change affecting dependencies was:
initial int64 support; now linking to RcppInt64, and importing bit64; FID and OFTInteger64 fields are now returned in R as integer64; updated the examples (2024-04-06)
I just noticed that the gdalvector branch source and documentation in the draft are not in-sync, i.e. getFeatureSet
is listed in the doc but not in the source:
https://usdaforestservice.github.io/gdalraster/articles/gdalvector-draft.html#class-gdalvector
Not unexpected for in-dev work but perhaps there's commits at your end not pushed for the branch? (Also fine of course just wondering, I've been unable to look at this in recent weeks)
btw I love how the class is self-documenting with print(GDALVector)
- I had no idea!!
just a note to think about, sometimes geometry doesn't have a name - in OGRSQL you can use "_ogrgeometry\" but I'm not sure that's really desirable.
tf <- tempfile(fileext = ".fgb")
ogr2ogr(src, tf, cl_arg = c("-f", "FlatGeobuf"))
ds <- new(GDALVector, tf)
ds$resetReading()
names(ds$getNextFeature())
## [1] "FID" "event_id" "incid_name" "incid_type" "map_id" "burn_bnd_ac"
## [7] "burn_bnd_lat" "burn_bnd_lon" "ig_date" "ig_year" "" ## << geometry name is empty
ds$close()
It just means that if we tibble::as_tibble(ds$getNextFeature())
we get an error "elements must be named" - but as.data.frame is ok it crafts a long name from the WKT (yukky). I don't know exactly what a good answer here would be but there might be some convention in GDAL itself. Also it's a development tool so not having a name also might be fine as a stance.
I just updated branch gdalvector from main. I generally try to keep it up to date, and I usually install from there for my local working instance of gdalraster. But it was a little behind after recent activity.
Also, note that the doc states:
The draft class definition below has been partially implemented in: https://github.com/USDAForestService/gdalraster/blob/gdalvector/src/gdalvector.cpp
The header file can be referenced for the public class methods that have been implemented so far in the prototype (a subset of the draft class definition below): https://github.com/USDAForestService/gdalraster/blob/gdalvector/src/gdalvector.h
So not everything in the draft spec has been implemented yet.
Status update: I'm working toward a CRAN release of gdalraster 1.11 in the next day or two. After that, I plan to refocus on gdalvector and work toward a release version for it in coming weeks if possible.
Excellent, thank you. I'll keep collecting feedback on GDALVector for when you get back to it.
Loving this package it's really excellent!!
The document "Draft Bindings to OGR" has been updated for recent changes and ongoing development in branch gdalvector.
Mainly, sub-headings were added under the heading Description of the interface. New specifications were added in the sub-heading "Feature retrieval" (adds GDALVector::fetch()
, analog of DBI::dbFetch()
, and support for all OGR field types mapped to native R types), and sub-heading "Geometry" (adds support for WKB, and a GDALVector
per-object setting returnGeomAs
).
The draft class definition has been updated under heading class GDALVector.
Code has been updated under Example: usage for class GDALVector to include the new additions, mainly fetch()
and the new setting for geometry return format. Note that the new field returnGeomAs
currently defaults to NONE
, and can be set back and forth as needed to one of WKT
, WKT_ISO
, WKB
, WKB_ISO
or NONE
. There are also two other new read/write fields, defaultGeomFldName
(applies when the geometry column name in the source layer is empty, like with shapefiles etc.) and wkbByteOrder
(LSB
the default, or MSB
).
Branch gdalvector
has been merged into main
and further development will be done there.
Documentation for class GDALVector
is now available in the development version of the package, and online at:
https://usdaforestservice.github.io/gdalraster/articles/gdalvector-draft.html
Note that the class methods for layer geoprocessing are not documented (layerIntersection()
etc.). Those are being re-implemented as a stand-alone ogr_proc()
interface that will accept objects of class GDALVector
as input, and create the output layer if needed.
awesome!
ogr_proc()
was added in #451, with online documentation at:
https://usdaforestservice.github.io/gdalraster/reference/ogr_proc.html
Feedback is welcome if we need improvements to the interface or other changes.
This issue is for potential discussion of bindings to the GDAL/OGR Vector API. A draft vector interface is described at: https://usdaforestservice.github.io/gdalraster/articles/gdalvector-draft.html
A prototype with partial implementation is in the
gdalvector
branch at: https://github.com/USDAForestService/gdalraster/tree/gdalvectorComments and suggestions are appreciated.