Getting started examples

sfarrow is designed to help read/write spatial vector data in “simple feature” format from/to Parquet files while maintaining coordinate reference system information. Essentially, this tool is attempting to connect R objects in sf and in arrow and it relies on these packages for its internal work.

A key goal is to support interoperability of spatial data in Parquet files. R objects (including sf) can be written to files with arrow; however, these do not necessarily maintain the spatial information or can be read in by Python. sfarrow implements a metadata format also used by Python GeoPandas, described here: https://github.com/geopandas/geo-arrow-spec. Note that these metadata are not stable yet, and sfarrow will warn you that it may change.

Reading and writing single files

A Parquet file (with .parquet extension) can be read using st_read_parquet() and pointing to the file system. This will create an sf spatial data object in memory which can then be used as normal using functions from sf.

# read an example dataset created from Python using geopandas
world <- st_read_parquet(system.file("extdata", "world.parquet", package = "sfarrow"))

class(world)
#> [1] "sf"         "data.frame"
world
#> Simple feature collection with 177 features and 5 fields
#> Geometry type: GEOMETRY
#> Dimension:     XY
#> Bounding box:  xmin: -180 ymin: -90 xmax: 180 ymax: 83.64513
#> Geodetic CRS:  WGS 84
#> First 10 features:
#>      pop_est     continent                     name iso_a3 gdp_md_est
#> 1     920938       Oceania                     Fiji    FJI  8.374e+03
#> 2   53950935        Africa                 Tanzania    TZA  1.506e+05
#> 3     603253        Africa                W. Sahara    ESH  9.065e+02
#> 4   35623680 North America                   Canada    CAN  1.674e+06
#> 5  326625791 North America United States of America    USA  1.856e+07
#> 6   18556698          Asia               Kazakhstan    KAZ  4.607e+05
#> 7   29748859          Asia               Uzbekistan    UZB  2.023e+05
#> 8    6909701       Oceania         Papua New Guinea    PNG  2.802e+04
#> 9  260580739          Asia                Indonesia    IDN  3.028e+06
#> 10  44293293 South America                Argentina    ARG  8.794e+05
#>                          geometry
#> 1  MULTIPOLYGON (((180 -16.067...
#> 2  POLYGON ((33.90371 -0.95, 3...
#> 3  POLYGON ((-8.66559 27.65643...
#> 4  MULTIPOLYGON (((-122.84 49,...
#> 5  MULTIPOLYGON (((-122.84 49,...
#> 6  POLYGON ((87.35997 49.21498...
#> 7  POLYGON ((55.96819 41.30864...
#> 8  MULTIPOLYGON (((141.0002 -2...
#> 9  MULTIPOLYGON (((141.0002 -2...
#> 10 MULTIPOLYGON (((-68.63401 -...
plot(sf::st_geometry(world))

Similarly, a Parquet file can be written from an sf object using st_write_parquet() and specifying a path to the new file. Non-spatial objects cannot be written with sfarrow, and users should instead use arrow.

# output the file to a new location
# note the warning about possible future changes in metadata.
st_write_parquet(world, dsn = file.path(tempdir(), "new_world.parquet"))
#> Warning: This is an initial implementation of Parquet/Feather file support and
#> geo metadata. This is tracking version 0.1.0 of the metadata
#> (https://github.com/geopandas/geo-arrow-spec). This metadata
#> specification may change and does not yet make stability promises.  We
#> do not yet recommend using this in a production setting unless you are
#> able to rewrite your Parquet/Feather files.

Partitioned datasets

While reading/writing a Parquet file is nice, the real power of arrow comes from splitting big datasets into multiple files, or partitions, based on criteria that make it faster to query. There is currently basic support in sfarrow for multi-file spatial datasets. For additional dataset querying options, see the arrow documentation.

Querying and reading Datasets

sfarrow accesses arrows’s dplyr interface to explore partitioned, Arrow datasets.

For this example we will use a dataset which was created by randomly splitting the nc.shp file first into three groups and then further partitioning into two more random groups. This creates a nested set of files.

list.files(system.file("extdata", "ds", package = "sfarrow"), recursive = TRUE)
#> [1] "split1=1/split2=1/part-3.parquet" "split1=1/split2=2/part-0.parquet"
#> [3] "split1=2/split2=1/part-1.parquet" "split1=2/split2=2/part-5.parquet"
#> [5] "split1=3/split2=1/part-2.parquet" "split1=3/split2=2/part-4.parquet"

The file tree is showing that the data were partitioned by the variables “split1” and “split2”. Those are the column names that were used for the random splits. This partitioning is in “Hive style” where the partitioning variables are in the paths.

The first step is to open the Dataset using arrow.

ds <- arrow::open_dataset(system.file("extdata", "ds", package="sfarrow"))

For small datasets (as in the example) we can read the entire set of files into an sf object.

nc_ds <- read_sf_dataset(ds)

nc_ds
#> Simple feature collection with 100 features and 16 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
#> Geodetic CRS:  NAD27
#> First 10 features:
#>     AREA PERIMETER CNTY_ CNTY_ID     NAME  FIPS FIPSNO CRESS_ID BIR74 SID74
#> 1  0.097     1.670  1833    1833 Hertford 37091  37091       46  1452     7
#> 2  0.062     1.547  1834    1834   Camden 37029  37029       15   286     0
#> 3  0.109     1.325  1841    1841   Person 37145  37145       73  1556     4
#> 4  0.081     1.288  1880    1880  Watauga 37189  37189       95  1323     1
#> 5  0.044     1.158  1887    1887   Chowan 37041  37041       21   751     1
#> 6  0.086     1.267  1893    1893   Yadkin 37197  37197       99  1269     1
#> 7  0.170     1.680  1903    1903 Guilford 37081  37081       41 16184    23
#> 8  0.118     1.601  1946    1946  Madison 37115  37115       58   765     2
#> 9  0.134     1.755  1958    1958    Burke 37023  37023       12  3573     5
#> 10 0.116     1.664  1964    1964 McDowell 37111  37111       56  1946     5
#>    NWBIR74 BIR79 SID79 NWBIR79 split1 split2                       geometry
#> 1      954  1838     5    1237      1      1 MULTIPOLYGON (((-76.74506 3...
#> 2      115   350     2     139      1      1 MULTIPOLYGON (((-76.00897 3...
#> 3      613  1790     4     650      1      1 MULTIPOLYGON (((-78.8068 36...
#> 4       17  1775     1      33      1      1 MULTIPOLYGON (((-81.80622 3...
#> 5      368   899     1     491      1      1 MULTIPOLYGON (((-76.68874 3...
#> 6       65  1568     1      76      1      1 MULTIPOLYGON (((-80.49554 3...
#> 7     5483 20543    38    7089      1      1 MULTIPOLYGON (((-79.53782 3...
#> 8        5   926     2       3      1      1 MULTIPOLYGON (((-82.89597 3...
#> 9      326  4314    15     407      1      1 MULTIPOLYGON (((-81.81628 3...
#> 10     134  2215     5     128      1      1 MULTIPOLYGON (((-81.81628 3...

With large datasets, more often we will want query them and return a reduced set of the partitioned records. To create a query, the easiest way is to use dplyr::filter() on the partitioning (and/or other) variables to subset the rows and dplyr::select() to subset the columns. read_sf_dataset() will then use the arrow_dplyr_query and call dplyr::collect() to extract and then process the Arrow Table into sf.

nc_d12 <- ds %>% 
            filter(split1 == 1, split2 == 2) %>%
            read_sf_dataset()

nc_d12
#> Simple feature collection with 20 features and 16 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -83.36472 ymin: 34.71101 xmax: -75.45698 ymax: 36.58965
#> Geodetic CRS:  NAD27
#> First 10 features:
#>     AREA PERIMETER CNTY_ CNTY_ID        NAME  FIPS FIPSNO CRESS_ID BIR74 SID74
#> 1  0.114     1.442  1825    1825        Ashe 37009  37009        5  1091     1
#> 2  0.143     1.630  1828    1828       Surry 37171  37171       86  3188     5
#> 3  0.153     2.206  1832    1832 Northampton 37131  37131       66  1421     9
#> 4  0.118     1.421  1836    1836      Warren 37185  37185       93   968     4
#> 5  0.114     1.352  1838    1838     Caswell 37033  37033       17  1035     2
#> 6  0.143     1.663  1840    1840   Granville 37077  37077       39  1671     4
#> 7  0.108     1.483  1900    1900     Forsyth 37067  37067       34 11858    10
#> 8  0.111     1.392  1904    1904    Alamance 37001  37001        1  4672    13
#> 9  0.104     1.294  1907    1907      Orange 37135  37135       68  3164     4
#> 10 0.122     1.516  1932    1932    Caldwell 37027  37027       14  3609     6
#>    NWBIR74 BIR79 SID79 NWBIR79 split1 split2                       geometry
#> 1       10  1364     0      19      1      2 MULTIPOLYGON (((-81.47276 3...
#> 2      208  3616     6     260      1      2 MULTIPOLYGON (((-80.45634 3...
#> 3     1066  1606     3    1197      1      2 MULTIPOLYGON (((-77.21767 3...
#> 4      748  1190     2     844      1      2 MULTIPOLYGON (((-78.30876 3...
#> 5      550  1253     2     597      1      2 MULTIPOLYGON (((-79.53051 3...
#> 6      930  2074     4    1058      1      2 MULTIPOLYGON (((-78.74912 3...
#> 7     3919 15704    18    5031      1      2 MULTIPOLYGON (((-80.0381 36...
#> 8     1243  5767    11    1397      1      2 MULTIPOLYGON (((-79.24619 3...
#> 9      776  4478     6    1086      1      2 MULTIPOLYGON (((-79.01814 3...
#> 10     309  4249     9     360      1      2 MULTIPOLYGON (((-81.32813 3...
plot(sf::st_geometry(nc_d12), col="grey")

When using select() to read only a subset of columns, if the geometry column is not returned, the default behaviour of sfarrow is to throw an error from read_sf_dataset. If you do not need the geometry column for your analyses, then using arrow and not sfarrow should be sufficient. However, setting find_geom = TRUE in read_sf_dataset will read in any geometry columns in the metadata, in addition to the selected columns.

# this command will throw an error
# no geometry column selected for read_sf_dataset
# nc_sub <- ds %>% 
#             select('FIPS') %>% # subset of columns
#             read_sf_dataset()

# set find_geom
nc_sub <- ds %>%
            select('FIPS') %>% # subset of columns
            read_sf_dataset(find_geom = TRUE)

nc_sub
#> Simple feature collection with 100 features and 1 field
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
#> Geodetic CRS:  NAD27
#> First 10 features:
#>     FIPS                       geometry
#> 1  37091 MULTIPOLYGON (((-76.74506 3...
#> 2  37029 MULTIPOLYGON (((-76.00897 3...
#> 3  37145 MULTIPOLYGON (((-78.8068 36...
#> 4  37189 MULTIPOLYGON (((-81.80622 3...
#> 5  37041 MULTIPOLYGON (((-76.68874 3...
#> 6  37197 MULTIPOLYGON (((-80.49554 3...
#> 7  37081 MULTIPOLYGON (((-79.53782 3...
#> 8  37115 MULTIPOLYGON (((-82.89597 3...
#> 9  37023 MULTIPOLYGON (((-81.81628 3...
#> 10 37111 MULTIPOLYGON (((-81.81628 3...

Writing to Datasets

To write an sf object into multiple files, we can again construct a query using dplyr::group_by() to define the partitioning variables. The result is then passed to sfarrow.

world %>%
  group_by(continent) %>%
  write_sf_dataset(file.path(tempdir(), "world_ds"), 
                   format = "parquet",
                   hive_style = FALSE)
#> Warning: This is an initial implementation of Parquet/Feather file support and
#> geo metadata. This is tracking version 0.1.0 of the metadata
#> (https://github.com/geopandas/geo-arrow-spec). This metadata
#> specification may change and does not yet make stability promises.  We
#> do not yet recommend using this in a production setting unless you are
#> able to rewrite your Parquet/Feather files.

In this example we are not using Hive style. This results in the partitioning variable not being in the folder paths.

list.files(file.path(tempdir(), "world_ds"))
#> [1] "Africa"                  "Antarctica"             
#> [3] "Asia"                    "Europe"                 
#> [5] "North America"           "Oceania"                
#> [7] "Seven seas (open ocean)" "South America"

To read this style of Dataset, we must specify the partitioning variables when it is opened.

arrow::open_dataset(file.path(tempdir(), "world_ds"), 
                    partitioning = "continent") %>%
  filter(continent == "Africa") %>%
  read_sf_dataset()
#> Simple feature collection with 51 features and 5 fields
#> Geometry type: GEOMETRY
#> Dimension:     XY
#> Bounding box:  xmin: -17.62504 ymin: -34.81917 xmax: 51.13387 ymax: 37.34999
#> Geodetic CRS:  WGS 84
#> First 10 features:
#>     pop_est            name iso_a3 gdp_md_est continent
#> 1  53950935        Tanzania    TZA   150600.0    Africa
#> 2    603253       W. Sahara    ESH      906.5    Africa
#> 3  83301151 Dem. Rep. Congo    COD    66010.0    Africa
#> 4   7531386         Somalia    SOM     4719.0    Africa
#> 5  47615739           Kenya    KEN   152700.0    Africa
#> 6  37345935           Sudan    SDN   176300.0    Africa
#> 7  12075985            Chad    TCD    30590.0    Africa
#> 8  54841552    South Africa    ZAF   739100.0    Africa
#> 9   1958042         Lesotho    LSO     6019.0    Africa
#> 10 13805084        Zimbabwe    ZWE    28330.0    Africa
#>                          geometry
#> 1  POLYGON ((33.90371 -0.95, 3...
#> 2  POLYGON ((-8.66559 27.65643...
#> 3  POLYGON ((29.34 -4.499983, ...
#> 4  POLYGON ((41.58513 -1.68325...
#> 5  POLYGON ((39.20222 -4.67677...
#> 6  POLYGON ((24.56737 8.229188...
#> 7  POLYGON ((23.83766 19.58047...
#> 8  POLYGON ((16.34498 -28.5767...
#> 9  POLYGON ((28.97826 -28.9556...
#> 10 POLYGON ((31.19141 -22.2515...