vect() incorrectly returns zero geometries when using filter option

[bengannon-fc]

First, thank you to the development team for advancing the terra package. It is very useful.

I have encountered an issue where the vect() function incorrectly returns zero geometries when reading a point feature class using the filter option. I can read the same layer in without a filter and produce the correct number of features with the intersect() function. Unfortunately, this is not a practical fix for my workflow because of how much faster I can read in data from large feature classes using the filter option. I suspect that vect() function may be incorrectly testing for overlapping extents as a precursor to the spatial filter because I can get the process to work fine when I clip my points layer to a different extent in a GIS.

Example:

PODs is polygon SpatVector, i is my iterator for a looped process

PODs[i,]
# class       : SpatVector 
# geometry    : polygons 
# dimensions  : 1, 5  (geometries, attributes)
# extent      : -2181524, -2154610, 2442862, 2477299  (xmin, xmax, ymin, ymax)
# coord. ref. : USA_Contiguous_Albers_Equal_Area_Conic_USGS_version (ESRI:102039) 
# names       :    Id POD_ID Region     Acres   Area_ac
# type        : <int>  <int>  <int>     <num>     <num>
# values      :     0    525      6 9.252e+04 9.252e+04

# Reading in with filter produces object with no features
inVect <- vect('./INPUT/VECTOR_INPUT.gdb',layer='cr_bldg_stru_pt',filter=PODs[i,])
inVect
# class       : SpatVector 
# geometry    : none 
# dimensions  : 0, 0  (geometries, attributes)
# extent      : 0, 0, 0, 0  (xmin, xmax, ymin, ymax)
# source      : VECTOR_INPUT.gdb (cr_bldg_stru_pt)
# coord. ref. : USA_Contiguous_Albers_Equal_Area_Conic_USGS_version (ESRI:102039) 

Reading in without filter works fine

inVect <- vect('./INPUT/VECTOR_INPUT.gdb',layer='cr_bldg_stru_pt')
inVect
# class       : SpatVector 
# geometry    : points 
# dimensions  : 63037, 34  (geometries, attributes)
# extent      : -11316921, 6920808, -72902.24, 7531605  (xmin, xmax, ymin, ymax)
# source      : VECTOR_INPUT.gdb (cr_bldg_stru_pt)
# coord. ref. : USA_Contiguous_Albers_Equal_Area_Conic_USGS_version (ESRI:102039) 
# names       :           CR_ID SURVEY_ID         GEOM_ID         RESNAME        BND_TYPE       BND_OTHER IS_EXTANT
# type        :           <chr>     <chr>           <chr>           <chr>           <chr>           <chr>     <chr>
# values      : {91BD962F-96E8~           {EA172C5C-F070~ Vista Avenue V~ Arbitrary point Point is unspe~   Unknown
#               {6AC088F7-6C47~           {B1CC00E2-B08B~  Astoria Column Arbitrary point Point is unspe~   Unknown
#               {4B723942-5706~           {4F43D9C4-049B~ Rock O' the Ra~ Arbitrary point Point is unspe~   Unknown

intRes <- intersect(inVect,PODs[i,])
intRes
# class       : SpatVector 
# geometry    : points 
# dimensions  : 49, 39  (geometries, attributes)
# extent      : -2180693, -2160744, 2445613, 2471349  (xmin, xmax, ymin, ymax)
# coord. ref. : USA_Contiguous_Albers_Equal_Area_Conic_USGS_version (ESRI:102039) 
# names       :           CR_ID SURVEY_ID         GEOM_ID         RESNAME        BND_TYPE       BND_OTHER IS_EXTANT
# type        :           <chr>     <chr>           <chr>           <chr>           <chr>           <chr>     <chr>
# values      : {D405996C-6E1D~           {D27621AF-F864~ Clark, Frank C~ Arbitrary point Point is unspe~   Unknown
#               {2B707EFF-650A~           {A31B2B82-FBBC~ Shone-Charley ~ Arbitrary point Point is unspe~   Unknown
#               {042D9A78-3FF3~           {8C307DEA-B876~    Barnum Hotel Arbitrary point Point is unspe~   Unknown

[rhijmans]

This seems to work in the example below.

library(terra)
f <- system.file("ex/lux.shp", package="terra")
v <- vect(f)
s <- spatSample(v, 1000)
fout <- paste0(tempfile(), ".shp")
writeVector(s, fout, overwrite=TRUE)

x <- vect(fout, filter = v[1,])

plot(v); points(s);  points(x, col="red")

Does it fail if you store "fout" in a GDB? (unfortunately, GDAL can not write these, so trying that would take some manual intervention, and an ArcGIS license)

[rhijmans]

This also works for me with a GDB. So this is a bit mysterious. Perhaps an old (version) GDB file?

[bengannon-fc]

Thank you Robert.

Your example works fine for me, including when I save the shapefile contents to a geodatabase and read it back into R. I'm not sure what counts as an "old" file GDB. I created this file GDB yesterday with ArcGIS 10.8.

I should have mentioned that I'm using terra version 1.5.34.

I have successfully used the filter method with vect() to read in geodatabase feature classes in several other workflows without obvious issues, so I don't question that the filter generally works. The reason I suspect an issue with a test for overlapping extent is that the failure occurred with data from United States-wide coverage (https://irma.nps.gov/DataStore/Reference/Profile/2210280), but everything works fine after I clipped my inputs to CONUS coverage. I don't think the issue is data size related because filter works in this same script on bigger feature classes. I also don't detect any corrupt geometries or other data anomalies in layers that don't load.

Surprisingly, I also avoid the error by exporting the original point feature class to a shapefile and reading it in from there.

[rhijmans]

I'm not sure what counts as an "old" file GDB

All I know is that the format has gradually changed over time. So I wondered if this was perhaps a GDB that was created many years ago; or with an old version of ArcMap. This does not appear to be the issue at hand.

Can you point me to the file you are using? I downloaded the gdb from the hyperlink above but that is something else (the crs is lonlat and the closest layer name is crbldg_pt).

[bengannon-fc]

Thank you Robert.

This behavior is very confusing. My idea about the feature class extent was a red herring. Clipping the data to a smaller extent changed some other property of the data. The attached example contrasts the results I get for two identical feature classes except one was created within a feature dataset and the other was created outside a feature dataset. The feature class created in the feature dataset does not properly load with the filter option, but the one created outside the feature dataset works fine. I can't see any differences in the feature class properties that would cause them to behave differently.

I have also observed that multipart polygon feature classes do not properly load with the filter option. The st_read() function in sf has the same behavior. Converting from multipart to singlepart in a GIS is an acceptable solution to me, but it might be worth expanding on the limitations of the filter method in the documentation so users are warned of this.

Failure_case.zip

[rhijmans]

This suggests that GDAL::SetSpatialFilter does not work well for layers in a "feature dataset" inside of a GDB. If this is consistent, it should be reported at GDAL.

Can you expand on your additional point that polygon classes do not properly load? Do you get too many features? Too few?

[bengannon-fc]

Thank you Robert. For now, I'll use filter cautiously. I have many feature classes properly loading from within feature datasets, so it is not a universal problem with feature datasets. It is disconcerting that a rather simple point feature class will not properly load with the filter method. When loading fails, it does so across all summary units, indicating it is a problem with the input layer instead of the spatial filter(s).

The polygon failure is similar to the point failure. The vect() function incorrectly returns an empty SpatVector with the filter option, but the data loads fine without the filter. I tried to create a simple example to demonstrate the multipart polygon issue, but the filter option behaves properly with my test case. Perhaps it is a geometry issue commonly associated with complex multipart polygon features and not a universal incompatibility with multipart features? I'll report this issue separately if I narrow in more on what the trigger is. The feature class I have issue with contains > 3,000 very complex multipart polygons. Common ArcGIS geometry clean up operations like repair geometry or applying a small buffer did not fix the issue. It only started loading after converting from multipart to singlepart.

[rhijmans]

Closing this as, for now, the problem is not sufficiently defined to take action.