HanOostdijk/HOQCwfs
Web Feature Service (WFS) read-only interface

Package index
Search the HanOostdijk/HOQCwfs package
Vignettes
Functions
60
Source code
20
Man pages
21
Home
/
GitHub
/
HanOostdijk/HOQCwfs
/

In HanOostdijk/HOQCwfs: Web Feature Service (WFS) read-only interface 

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(HOQCwfs)

Introduction {#introduction}

This package contains functions to retrieve information from a Web Feature Service (WFS) resource. The Open Geospatial Consortium WFS Interface Standard provides an interface allowing requests for geographical features across the web using platform-independent calls. This package generates such calls in the form of HTML requests in GET or POST format. Such a request is normally invisible to the user but can be made visible by adding httrverbose=T to the function call. See Debugging for details.

In this vignette we use a WFS data source that contains information about the trees in the municipality where I live. Apart from the name of the species each record also contains the exact location of the tree.

We will start with an example that retrieves trees of a certain species in a certain bounding box.
To do this we use the GetFeature operation of WFS that is called by the package function WFS_getfeature. Assuming we did not do a misspecification the result will be a simple features (sf) object. An sf object (managed by the R package sf) is a data.frame with spatial attributes. Instead of reading the whole data set and doing a filtering on the resulting data.frame we can add the filter to the GetFeature operation and retrieve only the part that we need. For large files this is obviously more efficient.

From the contents you can deduce that we introduce the package with a typical use case for the WFS_getfeature function (GetFeature operation). After that, we show how to find out what can be done with the WFS resource by using the WFS_getcapabilities (GetCapabilities operation) function to see which 'sub datasets' (featuretypes) can be used and with the WFS_describetype (DescribeFeatureType operation) to see which fields a featuretype contains.
After that we gradually introduce the possibilities to fine-tune the WFS_getfeature results in the section Retrieve the actual data.
The section Debugging describes some ways to extract more information when a query delivers no result or an unexpected one.
The package contains some tools to build spatial queries that are described in the section Spatial Tools and this vignette ends with some References to resources on the internet.

Introductory example {#intro-example}

In this first example we combine a text filter with a spatial filter: we want to retrieve some attributes of all trees of the species Alnus glutinosa 'Laciniata' in a certain area. So we need the call to WFS_getfeature that is done in the lines #11 and #12. I will start by mentioning the arguments that are not specified there: 

# devtools::install_github("HanOostdijk/HOQCwfs",build_vignettes=TRUE)       #1

library(HOQCwfs)                                                             #2
typename  <- 'topp:gidw_groenbomen'                                          #3
fields    <- 'boom_omschrijf,objec_omschrijf'                                #4
tree_name <- "Alnus glutinosa ''Laciniata''"                                 #5 
( tree_exp  <- glue::glue("boom_omschrijf='{tree_name}'") )                  #6
bbox_wgs84 <- c(4.860793, 52.313319, 4.861587, 52.316493 )                   #7
bbox_spec <- paste0(glue::glue_collapse(bbox_wgs84,sep=','),",'EPSG:4326'")  #8
( bbox_spec <- paste0('bbox(geometrie,',bbox_spec,')') )                     #9
cql_spec  <- paste0(tree_exp,' and ',bbox_spec )                             #10
wfs_ie    <-  WFS_getfeature(typename,cql_filter=cql_spec,                   #11
                             propertyname=fields,srsName='EPSG:4326' )       #12
dim(wfs_ie)                                                                  #13
print(wfs_ie,n=2,digits=5)                                                   #14

To see the expression that will be sent to the WFS server the argument httrverbose = T can be added. This will give some lines of information about the request. Only the first two lines are shown here:

wfs_ie    <-  WFS_getfeature(typename,cql_filter=cql_spec,                   
                             propertyname=fields,srsName='EPSG:4326',
                             httrverbose = T) 

x = HOQCutil::capture.output.both(
              WFS_getfeature(typename,cql_filter=cql_spec,                   
                             propertyname=fields,srsName='EPSG:4326',
                             httrverbose = T)  )
x = purrr::map_chr(x[1:2],URLdecode) 
HOQCutil::cap.out(cat(x))

[back to 'contents']

Workflow {#workflow}

A typical workflow looks like this :

Initialize {#initialize}

Load the package and indicate which WFS resource you want to access and the version that you want to use for the requests and their answers. In this way you can avoid specifying this in the separate function calls.
Also specify the separator that by default will be used in generated XML fragments.

library(HOQCwfs)
WFS_set_url()      # default https://geoweb.amstelveen.nl/geoserver/topp/wfs will be used
WFS_set_version()  # default 1.1.0 will be used
WFS_set_sep()      # default '\n' will be used

Look at the possibilities the WFS resource has to offer {#look-at-possibilities}

Retrieve the GetCapabilities document that describes the various features of the WFS resource and the actions that can be done. Because I have only read-access to WFS resources I am only interested in the query actions. An 'impression' of the rather extensive GetCapabilities xml document is given below.

xmlcap <- WFS_getcapabilities() # with the defaults for version and url 
class(xmlcap)

HOQCutil::cap.out(print(xmlcap),width=110)

The xml document can be written to a local dataset with the out_path argument:

xmlcap <- WFS_getcapabilities(out_path='mycap.xml')

[back to 'contents']

Look at the various features of the WFS resource {#look-at-features}

From the GetCapabilities document we can extract the features that are available:

ft1 <- WFS_featuretypes(xmlcap)
dim(ft1)

ft2 <- WFS_featuretypes(xmlcap,
       filternames=stringr::fixed("bomen", ignore_case = T))

When there are many features (r dim(ft1)[1] in this case) and we are just interested in feature names that contain the characters 'bomen' (i.e. 'trees') we can filter the names in the following way that results in only r dim(ft2)[1] cases



class(ft2)
dim(ft2)
ft2$layer

str(ft2[3,]) # show only featuretype 'topp:gidw_groenbomen'

Note that per featuretype the default Coordinate Reference System is indicated and coordinates of its bounding box. Here these coordinates are denoted not in this default CRS but in terms of WGS 84 (the 'standard' CRS in degrees). 

bomen_fields_1_3<- WFS_describefeaturetype(ft2$layer[c(1,3)])

[back to 'contents']

Look at the data fields of featuretypes {#look-at-data-fields}

For this we use the function WFS_describefeaturetype that executes the WFS DescribeFeatureType request.

Because we later will use the fields of topp:gidw_groenbomen we will show all fields for this featuretype. Because we request only one featuretype we can save the result of DescribeFeatureType in an xml file. Then we have the standard output:

WFS_describefeaturetype("topp:gidw_groenbomen",out_path = 'desc.xml')

and the desc.xml file with its contents:

HOQCutil::cap.out(cat(readLines("desc.xml",warn = F),sep='\n'), width=110) 

unlink('desc.xml')

[back to 'contents']

Retrieve the actual data {#retrieve}

Now that we know the metadata of our resource we can actually retrieve the attributes and geometrical data of the features. We do this with WFS_getfeature to:

If we are only interested in the number of features that would be returned we can add resultType=hits to any of the above. This only returns this number without any other data.

In the following sections we will give some examples.

[back to 'contents']

retrieve all features with all attributes and the geometry {#retrieve-all}
typename  <- 'topp:gidw_groenbomen'
wfs1      <-  WFS_getfeature(typename)
class(wfs1)
print(wfs1,n=5) # show only first 5 features

[back to 'contents'] [back to 'Retrieve the actual data']

retrieve some features {#retrieve-some-features}

For non-complex queries the cql_filter can be used. The features for the tree Prunus serrulata 'Kanzan' can be queried with the following code. The number of quotes in a query is often confusing (for me at least): here a set of quotes is needed because boom_omschrijf is a character field and quotes around 'Kanzan' have to be duplicated. 

typename  <- 'topp:gidw_groenbomen'
tree_name <- "Prunus serrulata ''Kanzan''"  
tree_exp  <- glue::glue("boom_omschrijf='{tree_name}'")
print(tree_exp)
wfs2      <-  WFS_getfeature(typename,cql_filter=tree_exp )
print(wfs2,n=5,digits=5) # show only first 5 features, 5 digits in geometry

[back to 'contents'] [back to 'Retrieve the actual data']

retrieve some attributes {#retrieve-attributes}

When we are only interested in the 'omschrijvingen' (descriptions) fields we can add the propertyname argument to request the fields boom_omschrijf and objec_omschrijf. Note that the id and geometry fields are always included.

typename  <- 'topp:gidw_groenbomen'
tree_name <- "Prunus serrulata ''Kanzan''"  
tree_exp  <- glue::glue("boom_omschrijf='{tree_name}'")
fields    <- 'boom_omschrijf,objec_omschrijf'
wfs3      <-  WFS_getfeature(typename,cql_filter=tree_exp,propertyname=fields )
print(wfs3,n=5,digits=5) # show only first 5 features, 5 digits in geometry

[back to 'contents'] [back to 'Retrieve the actual data']

retrieve the number of results {#number-of-features}

If we want to know beforehand the number of features that a query would return we can add the resultType=hits argument:

typename  <- 'topp:gidw_groenbomen'
tree_name <- "Prunus serrulata ''Kanzan''"  
tree_exp  <- glue::glue("boom_omschrijf='{tree_name}'")
fields    <- 'boom_omschrijf,objec_omschrijf'
  WFS_getfeature(typename,cql_filter=tree_exp,propertyname=fields,resultType='hits' )

[back to 'contents'] [back to 'Retrieve the actual data']

change the Spatial Reference System {#change-crs}

In the output of the previous examples we saw that the coordinates in the geometry were denoted in the Coordinate Reference System (CRS) that is used in the Netherlands: 'Amersfoort / RD New' also known as 'EPSG:28992'. If we want to have these coordinates denoted in the more commonly used 'WGS 84' (EPSG:4326) in terms of longitude and latitude we can add the srsName argument. SRS (Spatial Reference System) is a synonym for CRS.

typename  <- 'topp:gidw_groenbomen'
tree_name <- "Prunus serrulata ''Kanzan''"  
tree_exp  <- glue::glue("boom_omschrijf='{tree_name}'")
fields    <- 'boom_omschrijf,objec_omschrijf'
wfs4      <-  WFS_getfeature(typename,cql_filter=tree_exp,propertyname=fields,srsName='EPSG:4326' )
print(wfs4,n=5,digits=5) # show only first 5 features, 5 digits in geometry

[back to 'contents'] [back to 'Retrieve the actual data']

retrieve a subset of the query results with startindex, maxFeatures and count and using verbose {#retrieve-subset}

With the startindex and maxFeatures (1.1.0) or count (2.0.0) argument a subset of a query is produced. The WFS__getfeature function will use the correct version. By adding the echo_request=T clause the generated request url will be displayed so that you can verify this: count is translated to maxFeatures . Another way to display the generated url is setting httrverbose argument to FALSE.
Notice that startIndex=2 indicates that the first feature to be retrieved is number 3. In other words startIndex indicates the number of records to skip.

typename  <- 'topp:gidw_groenbomen'
tree_name <- "Prunus serrulata ''Kanzan''"  
tree_exp  <- glue::glue("boom_omschrijf='{tree_name}'")
fields    <- 'boom_omschrijf,objec_omschrijf'
wfs5a     <-  WFS_getfeature(typename,cql_filter=tree_exp,startIndex=2,count=6 
                             ,propertyname=fields,srsName='EPSG:4326'
                             ,echo_request=T,version='1.1.0' )

HOQCutil::cap.out({
wfs5a     <-  WFS_getfeature(typename,cql_filter=tree_exp,startIndex=2,count=6 
                             ,propertyname=fields,srsName='EPSG:4326'
                             ,echo_request=T,version='1.1.0' )  
}, width=93)

print(wfs5a,n=5,digits=5) # show only first 5 features, 5 digits in geometry

wfs5b     <-  WFS_getfeature(typename,cql_filter=tree_exp,startIndex=2,maxfeatures=6 
                             ,propertyname=fields,srsName='EPSG:4326'
                             ,echo_request=T,version='2.0.0' )

HOQCutil::cap.out({
wfs5b     <-  WFS_getfeature(typename,cql_filter=tree_exp,startIndex=2,maxfeatures=6 
                             ,propertyname=fields,srsName='EPSG:4326'
                             ,echo_request=T,version='2.0.0' )
}, width=93)

identical(wfs5a,wfs5b)

[back to 'contents'] [back to 'Retrieve the actual data']

retrieve features in a bounding box {#retrieve_boundingbox}

When you want to retrieve all features in a bounding box the bbox= argument can be used. When specifying the coordinates of the bounding box it is assumed that the default CRS is used. The default CRS for typename is EPSG:28992 (see Look at the various features of the WFS resource ). Therefore when using WGS 84 (EPSG:4326) the EPSG code has to be added to the specification. Note that this only impacts the features that are selected and not their coordinates. Use srsName='EPSG:4326' when you want WGS 84 coordinates.
In the two examples below we are interested in the same area but in the first case we use EPSG:28992 coordinates and in the second WGS 84 coordinates. Therefore the same features are selected

bbox_wgs84 <- c(4.860793, 52.313319, 4.861587, 52.316493 )
bbox_28992 <- convert_bbox(bbox_wgs84,'EPSG:4326','EPSG:28992')  

(bbox_spec <- glue::glue_collapse(bbox_28992,sep=','))
typename  <- 'topp:gidw_groenbomen'
fields    <- 'boom_omschrijf,objec_omschrijf'
wfs6a      <-  WFS_getfeature(typename,bbox=bbox_spec,
                             propertyname=fields )

(bbox_spec <- paste0(glue::glue_collapse(bbox_wgs84,sep=','),',EPSG:4326'))
wfs6b      <-  WFS_getfeature(typename,bbox=bbox_spec,
                             propertyname=fields )
identical(wfs6a,wfs6b)

[back to 'contents'] [back to 'Retrieve the actual data']

retrieve features in a bounding box with additional filtering {#retrieve_boundingbox-cqlfilter}

A combination of cql_filter= or filter= with bbox= is not allowed. In that case the bbox filtering has to be included in the cql_filter= or filter= construct by using the bbox function. First we show this for the cql_filter= construct for the two CRS cases.
Note the quotes around EPSG:4326 that were not necessary in the bbox= case.

typename  <- 'topp:gidw_groenbomen'
fields    <- 'boom_omschrijf,objec_omschrijf'
tree_name <- "Alnus glutinosa ''Laciniata''"  
tree_exp  <- glue::glue("boom_omschrijf='{tree_name}'")
bbox_spec <- glue::glue_collapse(bbox_28992,sep=',')
cql_spec<- paste0(tree_exp,' and bbox(geometrie,',bbox_spec,')') 

cat("\n",cql_spec,"\n\n")

HOQCutil::cap.out({
cat("\n",cql_spec,"\n\n")
}, width=93)

wfs7a      <-  WFS_getfeature(typename,cql_filter=cql_spec,
                             propertyname=fields )
print(wfs7a,n=5,digits=5)

bbox_spec <- paste0(glue::glue_collapse(bbox_wgs84,sep=','),",'EPSG:4326'")
cql_spec<- paste0(tree_exp,' and bbox(geometrie,',bbox_spec,')') 

cat("\n",cql_spec,"\n\n")

HOQCutil::cap.out({
cat("\n",cql_spec,"\n\n")
}, width=93)

wfs7b      <-  WFS_getfeature(typename,cql_filter=cql_spec,
                             propertyname=fields,srsName='EPSG:4326' )
print(wfs7b,n=5,digits=5)

Instead of cql_filter= we can also use filter= . In that case we have to specify the query in xml terms. To do this we can use the build_filter with the fg and bg functions. In the xml_query1 filter we specify the filter in the default CRS; when we use it in wfs7c we specify that we want the coordinates in longitude/latitude (EPSG:4326). In xml_query2 and wfs7d we do it the other way around.
Note that xml_query1 and xml_query2 are also different because we use version '1.1.0' in the first and '2.0.0' in the second case.

xml_query1=build_filter(version='1.1.0',
     fg("And"
           , propeq_xml('topp:boom_omschrijf',"Alnus glutinosa 'Laciniata'")
           , bbox_xml("geometrie","EPSG:28992",bbox_28992)
         )
  )
cat(xml_query1)
typename  <- 'topp:gidw_groenbomen'
fields    <- 'boom_omschrijf,objec_omschrijf'
wfs7c    <-  WFS_getfeature(typename
             ,version='1.1.0'
             ,filter=xml_query1
             ,propertyname=fields
             ,srsName='EPSG:4326')
print(wfs7c,n=5,digits=5)

xml_query2=build_filter(version='2.0.0',
     fg("And"
           , propeq_xml('topp:boom_omschrijf',"Alnus glutinosa 'Laciniata'")
           , bbox_xml("geometrie","EPSG:4326",bbox_wgs84)
         )
  )
cat(xml_query2)
wfs7d    <-  WFS_getfeature(typename
             ,version='2.0.0'
             ,filter=xml_query2
             ,propertyname=fields
             ,srsName='EPSG:28992')
print(wfs7d,n=5,digits=5)

[back to 'contents'] [back to 'Retrieve the actual data']

retrieve features with a spatial function {#retrieve-spatial-function}

There are number of operators in WFS that can be used to filter spatial features. As an example in this section we show how to filter the features within a distance of 50 meters from a given point using the R functions build_filter, spat_xml and spat_feature and WFS operator DWithin.
At the end of the example we use the sf::st_distance to check if the distance of the filtered features to the point is indeed less or equal to 50 meters.

The spat_xml function also works for the WFS operators Disjoint, Equals, DWithin, Beyond, Intersects, Touches, Crosses, Within, Contains, Overlaps and BBOX but in those cases the distance and units arguments are not used.

See Spatial Tools for more information about constructing filter queries.

version   <- '1.1.0'
coords    <- c(119103.4, 480726.0)
my_point  <- sf::st_sfc(sf::st_point(coords),crs='EPSG:28992')
print(my_point)
my_coords <- sf::st_coordinates(my_point)[,c('X','Y')]
print(my_coords)
xml_query <- build_filter(version=version,
     spat_xml('geometrie'
              ,spat_feature('point','EPSG:28992',my_coords) 
              ,spat_fun='DWithin'
              ,distance=50
              ,units='meters')
  )
typename  <- 'topp:gidw_groenbomen'
fields    <- 'boom_omschrijf,objec_omschrijf'
wfs7e     <-  WFS_getfeature(typename
             ,version=version                
             ,filter=xml_query
             ,propertyname=fields)

sf::st_distance(wfs7e,my_point,by_element=T)

HOQCutil::cap.out({
sf::st_distance(wfs7e,my_point,by_element=T)
}, width=95)

[back to 'contents'] [back to 'Retrieve the actual data']

order the features according to a sort key {#order-the-features}

To order features the SortBy arguments can be used. My experiences with Sortby in combination with filtering are the following:

Firstly we show the use of cql_filter for sorting one field: sorting on ascending (ASC) or descending (DESC) 'jaar' (i.e. 'year').

typename  <- 'topp:gidw_groenbomen'
fields    <- 'boom_omschrijf,jaar,objec_omschrijf'
tree_name <- "Alnus glutinosa ''Laciniata''"  
tree_exp  <- glue::glue("boom_omschrijf='{tree_name}'")
wfs8a      <-  WFS_getfeature(typename,cql_filter=tree_exp
                              ,sortby='jaar ASC'
                              ,propertyname=fields
                              ,count=5
)
print(wfs8a,n=5,digits=5)

wfs8b      <-  WFS_getfeature(typename,cql_filter=tree_exp
                              ,sortby='jaar DESC'
                              ,propertyname=fields
                              ,count=5
)
print(wfs8b,n=5,digits=5)

In cases where we have more than one sort field, we have to use filter. As an example we use here version '2.0.0' and the 'POST' method, but both versions and methods can be used in any combination.

xml_query1=build_filter(version='2.0.0'
           , propeq_xml('boom_omschrijf',"Alnus glutinosa 'Laciniata'")
  )
cat(xml_query1)
wfs8c      <-  WFS_getfeature(typename,filter=xml_query1
                              ,sortby='jaar ASC,objec_omschrijf DESC'
                              ,httrType = 'POST'
                              ,count=100
                              ,propertyname=fields
                              ,version='2.0.0'
)
print(wfs8c,n=5,digits=5)

[back to 'contents'] [back to 'Retrieve the actual data']

Debugging {#debugging}

When we retrieve the WFS data successfully, we normally want to end up with an sf object. We generate a request that is sent with the httr package to the WFS server. The server answers with a response object. We extract the json text and convert it to a json object with the jsonlite package or directly to an sf object with the sf package. In all these steps something can go wrong.

The WFS_getfeature function has some arguments that by default are set to FALSE but can be used to get (more) information when something appears to be wrong:

argument httrverbose {#debugging-httrverbose}

The logical values in argument httrverbose are passed to httr::verbose so that information about the connection between server and client can be returned . With the default value (FALSE) no such information is returned but specifying c(TRUE,TRUE,TRUE,TRUE) delivers more information. 

typename  <- 'topp:gidw_groenbomen'
fields    <- 'boom_omschrijf,jaar,objec_omschrijf'
tree_name <- "Alnus glutinosa ''Laciniata''"  
tree_exp  <- glue::glue("boom_omschrijf='{tree_name}'")

wfs8a      <-  WFS_getfeature(typename,cql_filter=tree_exp
                              ,sortby='jaar ASC'
                              ,propertyname=fields
                              ,httrverbose = c(T,T,T,T) # <--- 
                              ,count=2
)

HOQCutil::cap.out({
wfs8a      <-  WFS_getfeature(typename,cql_filter=tree_exp
                              ,sortby='jaar ASC'
                              ,propertyname=fields
                              ,httrverbose = c(T,T,T,T) # <--- 
                              ,count=2
  )
}, width=93)

[back to 'contents'] [back to 'Debugging']

argument echo_request {#debugging-echorequest}

If we want to see the request that is generated (when we have reason to think that it could be incorrect, as often we did when building the package) then we can set the echo_request argument to TRUE. The examples below show the difference between the request for the GET and the POST methods:

wfs8b      <-  WFS_getfeature(typename,cql_filter=tree_exp
                              ,sortby='jaar ASC'
                              ,httrType = 'GET'
                              ,propertyname=fields
                              ,echo_request = T # <--- 
                              ,count=2
)

HOQCutil::cap.out({
wfs8b      <-  WFS_getfeature(typename,cql_filter=tree_exp
                              ,sortby='jaar ASC'
                              ,httrType = 'GET'
                              ,propertyname=fields
                              ,echo_request = T # <--- 
                              ,count=2
  )
}, width=93)

wfs8c      <-  WFS_getfeature(typename,cql_filter=tree_exp
                              ,sortby='jaar ASC'
                              ,httrType = 'POST'
                              ,propertyname=fields
                              ,echo_request = T # <--- 
                              ,count=2
)

HOQCutil::cap.out({
wfs8c      <-  WFS_getfeature(typename,cql_filter=tree_exp
                              ,sortby='jaar ASC'
                              ,httrType = 'POST'
                              ,propertyname=fields
                              ,echo_request = T # <--- 
                              ,count=2
  )
}, width=93)

[back to 'contents'] [back to 'Debugging']

argument debug {#debugging-debug}

By setting the debug argument to TRUE only the 'response object' of the httr request will be returned. When everthing works out fine the response object will be converted to a json or sf object. When this conversion does not succeed, it may be worthwhile to study the response object.
The first example shows that the package sometimes will handle an error by providing an error message. In this telling us that it can not convert the string 'two' to the number 2 : 

wfs8d <- WFS_getfeature(typename,cql_filter=tree_exp
                              ,sortby='jaar ASC'
                              ,propertyname=fields
                              ,count= 'two'
  )

 cat(wfs8d) 

HOQCutil::cap.out({
 cat(wfs8d)
}, width=93)

If this is not the case we can study the response object as in the following example where we display the structure of the object, see that it is not an html error (because the response objects reports a message "Success: (200) OK" ) and by viewing the text we see that the WFS server returns an XML ExceptionReport.

wfs8e <- WFS_getfeature(typename,cql_filter=tree_exp
                              ,sortby='jaar ASC'
                              ,debug=T  # <--- 
                              ,propertyname=fields
                              ,count= 'two'
  )

 str(wfs8e) 

HOQCutil::cap.out({
 str(wfs8e)
}, width=93)

 httr::http_error(wfs8e) 

HOQCutil::cap.out({
 httr::http_error(wfs8e) 
}, width=93)

 httr::http_status(wfs8e)$message 

HOQCutil::cap.out({
 httr::http_status(wfs8e)$message 
}, width=93)

 httr::headers(wfs8e)$`content-type`  

HOQCutil::cap.out({
 httr::headers(wfs8e)$`content-type`  
}, width=93)

 httr::content(wfs8e, encoding = 'UTF-8', as  =  'text')

HOQCutil::cap.out({
 httr::content(wfs8e, encoding = 'UTF-8', as  =  'text')
}, width=93)

[back to 'contents'] [back to 'Debugging']

argument sfverbose {#debugging-sfverbose}

Setting the argument sfverbose to TRUE shows the message that the package sf produces when converting the response object to an sf object. 

wfs8f <- WFS_getfeature(typename,cql_filter=tree_exp
                              ,sortby='jaar ASC'
                              ,sfverbose = T  # <--- 
                              ,propertyname=fields
                              ,count= 2
  )

HOQCutil::cap.out({
wfs8f <- WFS_getfeature(typename,cql_filter=tree_exp
                              ,sortby='jaar ASC'
                              ,sfverbose = T  # <--- 
                              ,propertyname=fields
                              ,count= 2
  )
}, width=93)

[back to 'contents'] [back to 'Debugging']

Spatial Tools {#spatial-tools}

In this section we describe some tools (functions) that can be useful when formulating spatial queries. Even when it is often possible to download the whole data set and filter the resulting sf object with dplyr (in case of a query of the attributes) or sf (in case of a query on the geometry) it is more efficient to have this handled by the WFS server.
Spatial filters can be specified by both the cql_filter and filter argument. My preference is to use the filter option because I think that the resulting xml code is better structured. In the section retrieve features with a spatial function we gave already a simple example of a spatial query.

Basic query blocks {#basic-query-blocks}

Functions build_filter , fg and bg

The function build_filter (often using the functions fg and bg) is used to create a filter. An example of a retrieval of a type of Alnus trees within a given bounding box

bbox_28992 <- c(119103, 480726, 119160, 481078)
f1 <- build_filter(version='1.1.0',
    fg("And"
        , propeq_xml('topp:boom_omschrijf',"Alnus glutinosa 'Laciniata'")
        , bbox_xml("geometrie","EPSG:28992",bbox_28992)
      )
   )
cat(f1)

The build_filter function temporarily sets the default version to the one specified (here 1.1.0) and creates a filter with the xmlns definitions for this version. The inner functions( here propeq_xml and bbox_xml) will use the same version if not overruled.
The function fg accepts a 'tag' as first argument and optionally attributes for the starting tag (ta) and one or more fg or bg constructs as in the following example. And end tag is automatically appended. The function bg is simpler version of fg: no tag attributes and only a tag with its value (see the 'Literal' tag).

xml_phrase = fg(
   'PropertyIsLike'
   , bg('PropertyName', 'topp:objec_omschrijf')
   , bg('Literal', '.unmoolen')
   , ta = 'wildCard="*" singleChar="." escape="!"'
   , sep = "\n"
   )
cat(xml_phrase)
Functions propertyname_xml, propeq_xml and bbox_xml
(propertyname_xml('han','2.0.0'))
(propertyname_xml('han','1.1.0',nopref =F))
(propertyname_xml('han','1.1.0',nopref =T))
Function spat_feature to define a spatial feature

In retrieve features with a spatial function a example with spat_feature was already given. The spat_feature function defines a spatial component in xml terms that corresponds to an sf object: sf::st_point, sf::st_linestring, sf::st_polygon, sf::st_multipoint, sf::st_multilinestring, sf::st_multipolygon orsf::st_bbox (for Envelope) . After showing how to define each of these object types we give another example of usage:

# Envelope (bbox)
coords_wgs84 <- c( 4.86, 52.31,  4.867, 52.316)
coords       <- convert_points(coords_wgs84,"EPSG:4326","EPSG:28992",out_matrix=F)
cat(spat_feature('Envelope','EPSG:28992',coords,version='1.1.0'),'\n')  
cat(spat_feature('Envelope','EPSG:28992',coords,version='2.0.0'),'\n')  

# Point
coords    <- c(119103.4, 480726.0)
cat(spat_feature('Point','EPSG:28992',coords,version='1.1.0'),'\n') 
cat(spat_feature('Point','EPSG:28992',coords,version='2.0.0'),'\n')   

# LineString
coords_wgs84 <- c( 4.86, 52.31, 4.86, 52.316, 4.867, 52.316)
coords       <- convert_points(coords_wgs84,"EPSG:4326","EPSG:28992")
cat(spat_feature('LineString','EPSG:28992',coords,version='1.1.0'),'\n')
cat(spat_feature('LineString','EPSG:28992',coords,version='2.0.0'),'\n')

# Polygon without a hole
coords_wgs84 <- c( 4.86, 52.31, 4.86, 52.316, 4.867, 52.316, 4.86, 52.31)
coords       <- convert_points(coords_wgs84,"EPSG:4326","EPSG:28992")
cat(spat_feature('Polygon','EPSG:28992',list(coords),version='1.1.0'),'\n')
cat(spat_feature('Polygon','EPSG:28992',list(coords),version='2.0.0'),'\n')

# Polygon with a hole
out_wgs84    <- c( 4.86, 52.31, 4.86, 52.316, 4.88, 52.316, 4.88, 52.31, 4.86, 52.31)
out_28992    <- convert_points(out_wgs84,"EPSG:4326","EPSG:28992")
in_wgs84    <- c( 4.87, 52.312, 4.87, 52.314, 4.875, 52.314, 4.875, 52.312, 4.87, 52.312)
in_28992    <- convert_points(in_wgs84,"EPSG:4326","EPSG:28992")
cat(spat_feature('Polygon','EPSG:28992',list(out_28992,in_28992),version='1.1.0'),'\n')
cat(spat_feature('Polygon','EPSG:28992',list(out_28992,in_28992),version='2.0.0'),'\n')

# Multipoint
coords_wgs84 <- c( 4.86, 52.31, 4.86, 52.316, 4.88, 52.316)
coords_28992 <- convert_points(coords_wgs84,"EPSG:4326","EPSG:28992") 
cat(spat_feature('Multipoint','EPSG:28992',coords_28992,version='1.1.0'),'\n')
cat(spat_feature('Multipoint','EPSG:28992',coords_28992,version='2.0.0'),'\n')

# MultiLineString
coords_wgs84a <- c( 4.86, 52.31, 4.86, 52.316, 4.88, 52.316)
coords_28992a <- convert_points(coords_wgs84a,"EPSG:4326","EPSG:28992")
coords_wgs84b <- c( 4.87, 52.31, 4.87, 52.306, 4.89, 52.316)
coords_28992b <- convert_points(coords_wgs84b,"EPSG:4326","EPSG:28992")
coords_28992 <- list(coords_28992a,coords_28992b)
cat(spat_feature('MultiLineString','EPSG:28992',coords_28992,version='1.1.0'),'\n')
cat(spat_feature('MultiLineString','EPSG:28992',coords_28992,version='2.0.0'),'\n')

# MultiPolygon (with hole)
out_wgs84    <- c( 4.86, 52.31, 4.86, 52.316, 4.88, 52.316, 4.88, 52.31, 4.86, 52.31)
out_28992    <- convert_points(out_wgs84,"EPSG:4326","EPSG:28992")
in_wgs84    <- c( 4.87, 52.312, 4.87, 52.314, 4.875, 52.314, 4.875, 52.312, 4.87, 52.312)
in_28992    <- convert_points(in_wgs84,"EPSG:4326","EPSG:28992")

pol1 <- list(out_28992,in_28992) # first polygon (with hole)

out_wgs84    <-  c( 4.86, 52.30, 4.86, 52.306, 4.88, 52.306, 4.88, 52.30, 4.86, 52.30)
out_28992    <- convert_points(out_wgs84,"EPSG:4326","EPSG:28992")

pol2 <- list(out_28992)         # second polygon (without hole)

mpol1 <- list(pol1,pol2)        # multipolygon

cat(spat_feature('MultiPolygon','EPSG:28992',mpol1,version='1.1.0'),'\n')
cat(spat_feature('MultiPolygon','EPSG:28992',mpol1,version='2.0.0'),'\n')

In the readme of this package a filter is used that determines all features that lie in an area with the exception of those that lie in a certain part of the area. In this case the spatial area is a polygon with a hole. This shows the use of the spat_feature function in context:

# define outer en inner rectangle in WGS 84 coordinates
out_wgs84 <- c( 4.86,  52.31,  4.86, 52.316, 4.88,  52.316, 
                 4.88,  52.31,  4.86, 52.31)
in_wgs84  <- c( 4.87,  52.312, 4.87, 52.314, 4.875, 52.314, 
                 4.875, 52.312, 4.87, 52.312)
# convert coordinates form WGS84 (EPSG:4326) to the Dutch CRS EPSG:28992
out_28992 <- convert_points(out_wgs84,"EPSG:4326","EPSG:28992")
in_28992  <- convert_points(in_wgs84,"EPSG:4326","EPSG:28992")

# define the filter that selects all points that intersect the rectangle with hole
# use EPSG:28992 because that is the default CRS for this WFS server (data set)
xml_query <- build_filter(version='1.1.0',
     spat_xml('geometrie',
              spat_feature('Polygon','EPSG:28992',list(out_28992,in_28992)),
              spat_fun='Intersects')
  )
# set WFS server (not necessary here because this is the default)
WFS_set_url("https://geoweb.amstelveen.nl/geoserver/topp/wfs")
# set typename (dataset) and attributes that will be retrieved
typename  <- 'topp:gidw_groenbomen'
fields    <- 'boom_omschrijf,objec_omschrijf'
# retrieve the data from the WFS server using the filter
wfs1      <-  WFS_getfeature(typename
             ,filter=xml_query
             ,propertyname=fields)
# show number of trees and attributes
dim(wfs1)
# show only the first two features
head(wfs1,2)

Utility functions convert_points and convert_bbox {#utility-functions}

The functionsconvert_bbox and convert_points convert (the coordinates of) a bbox or a set of points from one CRS to another. 

bbox_28992 <- c(119103, 480726, 119160, 481078)
convert_bbox(bbox_28992,"EPSG:28992","EPSG:4326")
convert_points(bbox_28992,"EPSG:28992","EPSG:4326")
convert_points(bbox_28992,"EPSG:28992","EPSG:4326",out_matrix = F)

[back to 'contents']

References {#references}

[back to 'contents']



HanOostdijk/HOQCwfs documentation built on March 6, 2023, 8:18 a.m.

R Package Documentation

rdrr.io home R language documentation Run R code online

Browse R Packages

CRAN packages Bioconductor packages R-Forge packages GitHub packages

We want your feedback!

Note that we can't provide technical support on individual packages. You should contact the package authors for that.
GitHub issue tracker
ian@mutexlabs.com
Personal blog

 
Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.

Close