Date last run: 29Sep2020
Introduction
This document shows how to use a query to select only a few features of a WFS data set. With a query on non-spatial attributes I had no problems but when I tried to do something with the spatial attribute (the geometry) I ran into troubles. Trying to solve these I learned that the queries depend on the version of WFS (first I used the latest (2.0.0) because the latest version is the best (??) but I could get working some queries only by using version 1.1.0) and that queries can be coded no only as GET requests but also as POST requests.
Load packages
HOQCutil::silent_library(c('sf','httr','xml2','purrr','glue','jsonlite','tibble'))
GetCapabilities
Retrieve the GetCapabilities information.
base_url <- "https://geodata.nationaalgeoregister.nl/wijkenbuurten2019/wfs"
getcapabilities <- function(version,output_name=NULL){
url <- parse_url(base_url)
url$query <- list(service = "WFS",
version = version,
request = "GetCapabilities")
request <- build_url(url)
request
xml_doc <- xml2::read_xml((request),options="NOWARNING")
xml_ns_strip(xml_doc)
if (!is.null(output_name)) {
unlink(output_name)
write_xml(xml_doc,output_name)
}
xml_doc
}
xml_cap1 <- getcapabilities('1.1.0',r'(d:\data\magweg\getcapabilities_110.xml)' )
xml_cap2 <- getcapabilities('2.0.0',r'(d:\data\magweg\getcapabilities_200.xml)' )
FeatureTypes
Show which FeatureTypes the WFS data set contains with the default CRS. Because there is no difference in the following output for version 1.1.0 and 2.0.0 we show only the output for the first:
featuretypes <- function(xml_doc){
version <- xml_attr(xml_find_first(xml_doc,"//wfs:WFS_Capabilities"),"version")
if (version == '2.0.0') {
crs_clause1 <- ".//FeatureTypeList//FeatureType"
crs_clause2 <- ".//DefaultCRS"
} else {
crs_clause1 <- ".//FeatureType"
crs_clause2 <- ".//DefaultSRS"
}
map_dfr(xml_find_all(xml_doc, crs_clause1),
function(ft) {
c(layer=xml_text(xml_find_first(ft, ".//Name")),
defaultcrs=xml_text(xml_find_first(ft, crs_clause2)))
})
}
FeatureTypes <- featuretypes(xml_cap1)
knitr::kable(FeatureTypes,caption='FeatureTypes with default CRS')
Table: FeatureTypes with default CRS
| layer | defaultcrs |
|---|---|
| wijkenbuurten2019:cbs_buurten_2019 | urn:x-ogc:def:crs:EPSG:28992 |
| wijkenbuurten2019:cbs_wijken_2019 | urn:x-ogc:def:crs:EPSG:28992 |
| wijkenbuurten2019:gemeenten2019 | urn:x-ogc:def:crs:EPSG:28992 |
Operations and parameters
Show which operations (types of request) exist and the allowed values for the parameters. Not quite sure why the version argument is not mentioned for the GetFeature request. Because there is a difference between the available operations and the structure of the parameters we show the output for both versions.
op_parms <- function (xml_doc) {
version <- xml_attr(xml_find_first(xml_doc,"//wfs:WFS_Capabilities"),"version")
if (version == '2.0.0') {
value_clause <- "//ows:AllowedValues//ows:Value"
} else {
value_clause <- "//ows:Value"
}
ops <- map(xml_find_all(xml_doc, ".//ows:OperationsMetadata//ows:Operation"),
function(op) {
op1 = xml_find_all(op,
paste0(".//ows:Parameter[contains(@name, 'AcceptVersions')]",
value_clause) )
op1 = list(AcceptVersions=map_chr(op1,xml_text))
op2 = xml_find_all(op,
paste0(".//ows:Parameter[contains(@name,'AcceptFormats')]",
value_clause))
op2 = list(AcceptFormats=map_chr(op2,xml_text))
op3 = xml_find_all(op,
paste0(".//ows:Parameter[contains(@name, 'outputFormat')]",
value_clause))
op3 = list(outputFormat=map_chr(op3,xml_text))
op4 = xml_find_all(op,
paste0(".//ows:Parameter[contains(@name,'resultType')]",
value_clause))
op4 = list(resultType=map_chr(op4,xml_text))
res = c( Operation= xml_attr(op, "name",default=" ") ,
op1, op2, op3, op4)
})
names(ops) = map_chr(ops,~pluck(.,"Operation"))
ops = map(ops,~discard(.,c(T,F,F,F,F)))
HOQCutil::cap.out( str(ops,nchar.max = 500,vec.len = 20),width=95 )
}
op_parms(xml_cap1)
#> List of 4
#> $ GetCapabilities :List of 4
#> ..$ AcceptVersions: chr [1:2] "1.0.0" "1.1.0"
#> ..$ AcceptFormats : chr "text/xml"
#> ..$ outputFormat : chr(0)
#> ..$ resultType : chr(0)
#> $ DescribeFeatureType:List of 4
#> ..$ AcceptVersions: chr(0)
#> ..$ AcceptFormats : chr(0)
#> ..$ outputFormat : chr "text/xml; subtype=gml/3.1.1"
#> ..$ resultType : chr(0)
#> $ GetFeature :List of 4
#> ..$ AcceptVersions: chr(0)
#> ..$ AcceptFormats : chr(0)
#> ..$ outputFormat : chr [1:13] "text/xml; subtype=gml/3.1.1" "GML2" "KML" "application/gml+x
#> ml; version=3.2" "application/json" "application/vnd.google-earth.kml xml" "application/vnd.goo
#> gle-earth.kml+xml" "csv" "gml3" "gml32" "json" "text/xml; subtype=gml/2.1.2" "text/xml; subtype
#> =gml/3.2"
#> ..$ resultType : chr [1:2] "results" "hits"
#> $ GetGmlObject :List of 4
#> ..$ AcceptVersions: chr(0)
#> ..$ AcceptFormats : chr(0)
#> ..$ outputFormat : chr(0)
#> ..$ resultType : chr(0)
op_parms(xml_cap2)
#> List of 6
#> $ GetCapabilities :List of 4
#> ..$ AcceptVersions: chr [1:3] "1.0.0" "1.1.0" "2.0.0"
#> ..$ AcceptFormats : chr "text/xml"
#> ..$ outputFormat : chr(0)
#> ..$ resultType : chr(0)
#> $ DescribeFeatureType :List of 4
#> ..$ AcceptVersions: chr(0)
#> ..$ AcceptFormats : chr(0)
#> ..$ outputFormat : chr "text/xml; subtype=gml/3.2"
#> ..$ resultType : chr(0)
#> $ GetFeature :List of 4
#> ..$ AcceptVersions: chr(0)
#> ..$ AcceptFormats : chr(0)
#> ..$ outputFormat : chr [1:13] "text/xml; subtype=gml/3.2" "GML2" "KML" "application/gml+xml
#> ; version=3.2" "application/json" "application/vnd.google-earth.kml xml" "application/vnd.googl
#> e-earth.kml+xml" "csv" "gml3" "gml32" "json" "text/xml; subtype=gml/2.1.2" "text/xml; subtype=g
#> ml/3.1.1"
#> ..$ resultType : chr [1:2] "results" "hits"
#> $ GetPropertyValue :List of 4
#> ..$ AcceptVersions: chr(0)
#> ..$ AcceptFormats : chr(0)
#> ..$ outputFormat : chr(0)
#> ..$ resultType : chr(0)
#> $ ListStoredQueries :List of 4
#> ..$ AcceptVersions: chr(0)
#> ..$ AcceptFormats : chr(0)
#> ..$ outputFormat : chr(0)
#> ..$ resultType : chr(0)
#> $ DescribeStoredQueries:List of 4
#> ..$ AcceptVersions: chr(0)
#> ..$ AcceptFormats : chr(0)
#> ..$ outputFormat : chr(0)
#> ..$ resultType : chr(0)
FeatureType wijkenbuurten2019:cbs_buurten_2019
Show the fields for this FeatureType. Because there is no difference in the available fields we only show the output for the version 1.1.0 code.
describefeaturetype <- function(version,typename,output_name=NULL){
url <- parse_url(base_url)
url$query <- list(service = "WFS",
version = version,
request = "DescribeFeatureType",
typename = typename)
request <- build_url(url)
xml_doc <- xml2::read_xml((request),options="NOWARNING")
xml_ns_strip(xml_doc)
# write_xml(xml_doc,r'(d:\data\magweg\describe_feature_type.txt)' )
if (!is.null(output_name)) {
unlink(output_name)
write_xml(xml_doc,output_name)
}
fieldnames <- xml_children(
xml_find_first(xml_doc,
paste0(".//xsd:extension[contains(@base,'gml:AbstractFeatureType')]",
"//xsd:sequence") )
)
map_chr(fieldnames,~xml_attr(.,'name'))
}
# version 1.1.0
describefeaturetype(version='1.1.0',
typename = "wijkenbuurten2019:cbs_buurten_2019",
output_name = r'(d:\data\magweg\describefeaturetype_110.txt)')
#> [1] "buurtcode"
#> [2] "jrstatcode"
#> [3] "buurtnaam"
#> [4] "wijkcode"
#> [5] "wijknaam"
#> [6] "gemeentecode"
#> [7] "gemeentenaam"
#> [8] "ind_wbi"
#> [9] "water"
#> [10] "meest_voorkomende_postcode"
#> [11] "dekkingspercentage"
#> [12] "omgevingsadressendichtheid"
#> [13] "stedelijkheid_adressen_per_km2"
#> [14] "bevolkingsdichtheid_inwoners_per_km2"
#> [15] "aantal_inwoners"
#> [16] "mannen"
#> [17] "vrouwen"
#> [18] "percentage_personen_0_tot_15_jaar"
#> [19] "percentage_personen_15_tot_25_jaar"
#> [20] "percentage_personen_25_tot_45_jaar"
#> [21] "percentage_personen_45_tot_65_jaar"
#> [22] "percentage_personen_65_jaar_en_ouder"
#> [23] "percentage_ongehuwd"
#> [24] "percentage_gehuwd"
#> [25] "percentage_gescheid"
#> [26] "percentage_verweduwd"
#> [27] "aantal_huishoudens"
#> [28] "percentage_eenpersoonshuishoudens"
#> [29] "percentage_huishoudens_zonder_kinderen"
#> [30] "percentage_huishoudens_met_kinderen"
#> [31] "gemiddelde_huishoudsgrootte"
#> [32] "percentage_westerse_migratieachtergrond"
#> [33] "percentage_niet_westerse_migratieachtergrond"
#> [34] "percentage_uit_marokko"
#> [35] "percentage_uit_nederlandse_antillen_en_aruba"
#> [36] "percentage_uit_suriname"
#> [37] "percentage_uit_turkije"
#> [38] "percentage_overige_nietwestersemigratieachtergrond"
#> [39] "oppervlakte_totaal_in_ha"
#> [40] "oppervlakte_land_in_ha"
#> [41] "oppervlakte_water_in_ha"
#> [42] "geom"
# version 2.0.0
# describefeaturetype(version='2.2.0',
# typename = "wijkenbuurten2019:cbs_buurten_2019",
# output_name = r'(d:\data\magweg\describefeaturetype_220.txt)')
Do a non-spatial query on wijkenbuurten2019:cbs_buurten_2019
Use the gemeentenaam (name of municipality) field to retrieve the buurten (neighbourhoods) of one municipality and plot these buurten.
A non-spatial query with cql_filter
Using both version 1.1.0 and 2.0.0 we see that the results are the same. Therefore we only do the plot for the code with the first version. Because we did not specify a CRS the default one (EPSG::28992) is used.
cqlfilter <- function(version,plotdata=F){
url <- parse_url(base_url)
url$query <- list(service = "WFS",
version = version,
request = "GetFeature",
# cql_filter = URLencode("gemeentenaam='Amstelveen'"),
cql_filter = "gemeentenaam='Amstelveen'",
typename = "wijkenbuurten2019:cbs_buurten_2019",
outputFormat = "application/json")
request <- build_url(url)
res_data <- NULL ;
HOQCutil::cap.out( res_data <- st_read(request),width=95 )
if (plotdata) plot(res_data[,1],main='Gemeente Amstelveen')
res_data
}
res_data <- cqlfilter('1.1.0',plotdata=T)
#> Reading layer `OGRGeoJSON' from data source `https://geodata.nationaalgeoregister.nl/wijkenbuu
#> rten2019/wfs?service=WFS&version=1.1.0&request=GetFeature&cql_filter=gemeentenaam%3D%27Amstelve
#> en%27&typename=wijkenbuurten2019%3Acbs_buurten_2019&outputFormat=application%2Fjson' using driv
#> er `GeoJSON'
#> Simple feature collection with 47 features and 42 fields
#> geometry type: MULTIPOLYGON
#> dimension: XY
#> bbox: xmin: 114538.3 ymin: 472779.6 xmax: 122404.5 ymax: 482651.8
#> projected CRS: Amersfoort / RD New
res_data <- cqlfilter('2.0.0',plotdata=F)
#> Reading layer `OGRGeoJSON' from data source `https://geodata.nationaalgeoregister.nl/wijkenbuu
#> rten2019/wfs?service=WFS&version=2.0.0&request=GetFeature&cql_filter=gemeentenaam%3D%27Amstelve
#> en%27&typename=wijkenbuurten2019%3Acbs_buurten_2019&outputFormat=application%2Fjson' using driv
#> er `GeoJSON'
#> Simple feature collection with 47 features and 42 fields
#> geometry type: MULTIPOLYGON
#> dimension: XY
#> bbox: xmin: 114538.3 ymin: 472779.6 xmax: 122404.5 ymax: 482651.8
#> projected CRS: Amersfoort / RD New
In the output of the DescribeFeatureType operation we saw that in the WFS environment there were 42 attributes present and that the ‘geometry’ attribute was named geom . When we look at the attributes of the Simple Feature result we see that an id attribute has been added and that the ‘geometry’ attribute is now named geometry :
names(res_data)
#> [1] "id"
#> [2] "buurtcode"
#> [3] "jrstatcode"
#> [4] "buurtnaam"
#> [5] "wijkcode"
#> [6] "wijknaam"
#> [7] "gemeentecode"
#> [8] "gemeentenaam"
#> [9] "ind_wbi"
#> [10] "water"
#> [11] "meest_voorkomende_postcode"
#> [12] "dekkingspercentage"
#> [13] "omgevingsadressendichtheid"
#> [14] "stedelijkheid_adressen_per_km2"
#> [15] "bevolkingsdichtheid_inwoners_per_km2"
#> [16] "aantal_inwoners"
#> [17] "mannen"
#> [18] "vrouwen"
#> [19] "percentage_personen_0_tot_15_jaar"
#> [20] "percentage_personen_15_tot_25_jaar"
#> [21] "percentage_personen_25_tot_45_jaar"
#> [22] "percentage_personen_45_tot_65_jaar"
#> [23] "percentage_personen_65_jaar_en_ouder"
#> [24] "percentage_ongehuwd"
#> [25] "percentage_gehuwd"
#> [26] "percentage_gescheid"
#> [27] "percentage_verweduwd"
#> [28] "aantal_huishoudens"
#> [29] "percentage_eenpersoonshuishoudens"
#> [30] "percentage_huishoudens_zonder_kinderen"
#> [31] "percentage_huishoudens_met_kinderen"
#> [32] "gemiddelde_huishoudsgrootte"
#> [33] "percentage_westerse_migratieachtergrond"
#> [34] "percentage_niet_westerse_migratieachtergrond"
#> [35] "percentage_uit_marokko"
#> [36] "percentage_uit_nederlandse_antillen_en_aruba"
#> [37] "percentage_uit_suriname"
#> [38] "percentage_uit_turkije"
#> [39] "percentage_overige_nietwestersemigratieachtergrond"
#> [40] "oppervlakte_totaal_in_ha"
#> [41] "oppervlakte_land_in_ha"
#> [42] "oppervlakte_water_in_ha"
#> [43] "geometry"
From now on I will, to make this blog post somewhat more readable, omit the st_read output by specifying the quiet argument. Instead of that I will use a function that shows only the dimensions and CRS of the result:
show_res <- function (sf_object) {
cat(
paste0('dim: (', nrow(sf_object),', ',ncol(sf_object),
') crs: ', st_crs(sf_object)$input)
,'\n' )
}
A non-spatial query with a WFS filter
Do the same with a WFS filter. It looks more complicated but it is more powerful than the cql_filter method because the latter method is not available for all WFS data sources. The MapServer WFS Filter Encoding documentation page gives some examples of WFS filters. The page OGC® Filter Encoding 2.0 Encoding Standard – With Corrigendum seems to be the authoritative source about filter.
A difference in the sf object compared with the previous example is that the geometry is in terms of CRS EPSG:4326 (or WGS 84 the ‘standard’ coordinates in longitude and latitude degrees) because we specified that in the srsName parameter in the GET request list.
We can also use outputFormat = "text/xml; subtype=gml/3.2", but in that case we have a CRS with value NA.
my_filter = paste0(
"<Filter>",
"<PropertyIsEqualTo><PropertyName>wijkenbuurten2019:gemeentenaam</PropertyName>",
"<Literal>Amstelveen</Literal></PropertyIsEqualTo>",
"</Filter>")
wfsfilter <- function(version,plotdata=F){
url <- parse_url(base_url)
url$query <- list(service = "WFS",
version = version,
request = "GetFeature",
typename = "wijkenbuurten2019:cbs_buurten_2019",
srsName = "EPSG:4326" ,
outputFormat = "application/json",
filter = my_filter
)
request <- build_url(url)
res_data <- NULL ; res_data <- st_read(request,quiet=T)
}
res_data = wfsfilter(version='1.1.0') ; show_res(res_data)
#> dim: (47, 43) crs: WGS 84
res_data = wfsfilter(version='2.0.0') ; show_res(res_data)
#> dim: (47, 43) crs: WGS 84
A spatial query on wijkenbuurten2019:cbs_buurten_2019
I tried to use a WFS filter (as describe in the previous paragraph) to do a spatial query but could not get it working: I got an error indicating a java.lang.NullPointerException or the filter was ignored (i.e. all features were returned and not only the expected ones). During my ‘experiments’ I became aware of the differences between version 1.1.0 and 2.0.0 and while searching the internet I also saw the use of POST (instead of the GET that is more familiar to me) requests.
For testing the filter possibilities I tried to select only buurten within 500 meters within my own neighborhood by using the DWITHIN function. Because I normally use coordinates in degrees longitude and latitude I specify them in EPSG:4326 and convert them to CRS EPSG::28992 (because that is the default for this WFS data source) :
my_point = c (4.8652, 52.3091) # OL, NB (long, lat)
my_point = st_sfc(st_point(my_point))
my_point = st_sf(geometry=my_point,crs=4326)
st_transform(my_point,crs=28992)
#> Simple feature collection with 1 feature and 0 fields
#> geometry type: POINT
#> dimension: XY
#> bbox: xmin: 119400.6 ymin: 480254.4 xmax: 119400.6 ymax: 480254.4
#> projected CRS: Amersfoort / RD New
#> geometry
#> 1 POINT (119400.6 480254.4)
In the following subsections I show various query methods that got me the expected results.
Using the CQL_FILTER GET method.
cqlfilter2 <- function (version) {
url <- parse_url(base_url)
url$query <- list(service = "WFS",
version = version,
request = "GetFeature",
typename = "wijkenbuurten2019:cbs_buurten_2019",
outputFormat = "application/json",
cql_filter= "dwithin(wijkenbuurten2019:geom, point(119400 480254), 500, meters)"
)
request <- build_url(url)
res_data <- NULL ; res_data <- st_read(request,quiet=T)
}
res_data <- cqlfilter2('1.1.0') ; show_res(res_data)
#> dim: (5, 43) crs: Amersfoort / RD New
res_data <- cqlfilter2('2.0.0') ; show_res(res_data)
#> dim: (5, 43) crs: Amersfoort / RD New
So we see that there 5 neighborhoods within 500 meters of the starting point. That means: within a circle with radius 500 meters and with the starting point as centre you can find (part of) 5 different neighborhoods.
Using the WFS Filter POST method
I tried to create a spatial WFS filter with the GET method, but did not succeed.
I found a lot of examples with a POST request but only when I saw the article by Andy Teucher (ateucher) that I could get it working. And then first only for version 1.1.0 ! Only after carefully reading and comparing with a contribution by Wouter.Visscher in a GeoForum article I got it working.
Schemas for version 2.0.0 with examples can be found in http://schemas.opengis.net/wfs/2.0/ .
The following functions and variables are created:
handle_body: this does the actualPOSTrequestfilter110: an example of a spatial filter for version1.1.0. Contents wise the same asfilter200filter200: an example of a spatial filter for version2.0.0. Contents wise the same asfilter110getfeature_body: a function that builds thePOSTrequest. It includes (just copies) the relevant ‘filter’ and select and sort on feature attributes.
sep = '\n'
handle_body <- function (body,base_url,to_sf=T,returnres=F,verbose=F) {
res <- NULL
if (verbose ==T)
res <- POST(base_url,body=body,content_type("text/xml"),verbose() )
else
res <- POST(base_url,body=body,content_type("text/xml") )
res_data <- content(res,encoding = 'UTF-8',as='text')
if (grepl('application/json',res$headers$`content-type`,fixed = T)) {
if (to_sf)
r <- sf::read_sf(res_data,quiet=(verbose==F))
else
r <- jsonlite::fromJSON(res_data)
} else {
r<- res_data
}
if (returnres==T){
list(df=r,res=res)
} else {
r
}
}
filter110 <- paste0( c(
'<ogc:Filter>',
'<ogc:DWithin>',
'<ogc:PropertyName >wijkenbuurten2019:geom</ogc:PropertyName>',
'<gml:Point srsName="EPSG:28992">',
'<gml:coordinates cs="," decimal="." ts=" ">119400,480254</gml:coordinates>',
"</gml:Point>",
"<ogc:Distance ogc:units='meters'>500.0</ogc:Distance>",
"</ogc:DWithin>",
"</ogc:Filter>"
),collapse = sep)
filter110 <- paste0(filter110,sep)
filter200 <- paste0( c(
'<fes:Filter>',
'<fes:DWithin>',
'<fes:ValueReference>wijkenbuurten2019:geom</fes:ValueReference>',
'<gml:Point srsName="EPSG:28992">',
'<gml:coordinates cs="," decimal="." ts=" ">119400,480254</gml:coordinates>',
"</gml:Point>",
"<fes:Distance fes:units='meters'>500.0</fes:Distance>",
"</fes:DWithin>",
"</fes:Filter>"),
collapse=sep)
filter200 <- paste0(filter200,sep)
getfeature_body <- function(version,
outputFormat="application/json",
count = NULL,
hits = F,
maxFeatures = NULL,
startIndex = NULL,
includeschema =T,
sel_cols = c('omgevingsadressendichtheid',
'aantal_inwoners',
'wijknaam',
'buurtnaam',
'geom'),
sortby = c('wijknaam', 'DESC',
'aantal_inwoners', 'ASC'),
filter = NULL,
prefix="wijkenbuurten2019",
layer="cbs_buurten_2019",
srsName = "EPSG:28992",
verbose = F,
sep = '\n') {
if (version == '1.1.0') {
if (!is.null(maxFeatures)) count=maxFeatures
if (!is.null(count)) count=glue('maxFeatures="{count}"') else count=''
startIndex = ''
ogcfes ='ogc';sufwfs='';suffes='';sufgml=''
schemeinfo = ''
typeNames = 'typeName' ; propref = "PropertyName"
extraogc = ''
} else {
if (!is.null(maxFeatures)) count=maxFeatures
if (!is.null(count)) count=glue('count="{count}"') else count=''
if (!is.null(startIndex)) startIndex=glue('startIndex="{startIndex}"')
else startIndex=''
ogcfes ='fes';sufwfs='/2.0';suffes='/2.0';sufgml='/3.2'
if (includeschema==T && version == '2.0.0') {
schemeinfo = c(
'xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" ',sep,
'xsi:schemaLocation=',sep,
' "http://www.opengis.net/wfs',sufwfs,sep,
' http://schemas.opengis.net/wfs/2.0/wfs.xsd"',sep
)
} else schemeinfo = ''
typeNames = 'typeNames' ; propref = "ValueReference" # "PropertyValue" #
extraogc = paste0('xmlns:ogc="http://www.opengis.net/ogc" ',sep)
}
# select columns
if (!is.null(sel_cols) ) {
sel_cols <- glue::glue_collapse(
glue::glue('<PropertyName>{prefix}:{sel_cols}</PropertyName>'),
sep=sep
)
} else sel_cols <- ''
# sortby columns :
if (!is.null(sortby) ) {
sortby1 <- matrix(sortby,ncol=2,byrow=T)
sortby2 <- paste0(
'<{ogcfes}:SortProperty>',
'<{ogcfes}:{propref}>{prefix}:{naam}</{ogcfes}:{propref}>',
'<{ogcfes}:SortOrder>{order}</{ogcfes}:SortOrder>',
'</{ogcfes}:SortProperty>',collapse = ''
)
sortby3 <- glue::glue(sortby2,naam=sortby1[,1],order=sortby1[,2],)
sortby <- paste0(
glue::glue('<{ogcfes}:SortBy>{sep}'),
glue::glue_collapse(sortby3,sep=sep),
glue::glue('{sep}</{ogcfes}:SortBy>'), sep=sep, collapse=sep)
} else sortby <- ''
# filter
if (is.null(filter) ) { filter <- '' }
if (hits==T) hits=glue('resultType="hits"') else hits=''
x= glue::glue(
'<?xml version="1.0" encoding="ISO-8859-1"?>',sep,
'<wfs:GetFeature ',sep,
'service="WFS" version="{version}" {count} {startIndex} ',
'outputFormat="{outputFormat}" {hits} ',sep,
'xmlns:{prefix}="http://wijkenbuurten2019.geonovum.nl" ',sep,
'xmlns:wfs="http://www.opengis.net/wfs{sufwfs}" ',sep,
'xmlns:{ogcfes}="http://www.opengis.net/{ogcfes}{suffes}" ',sep,
extraogc,
'xmlns:gml="http://www.opengis.net/gml{sufgml}" ',sep,
paste0(schemeinfo,collapse = ''),
'>',sep,
'<wfs:Query {typeNames}="{prefix}:{layer}" srsName="{srsName}">',sep,
sel_cols,sep,
sortby,
filter,
'</wfs:Query>',sep,
'</wfs:GetFeature>'
)
x= gsub('\n\n','\n',as.character(x),fixed=T)
if (verbose) cat('\n',x,'\n')
x
}
Generate a POST request by calling the getfeature_body with all the default arguments except for the filter argument. Show the generated POST request and the resulting sf object. In the first example we use version 1.1.0 syntax and we show as much output as can be produced. However a large part of that is formed by the coordinates of the buurten so we cap that part to 300 characters.
my_body110 <- getfeature_body('1.1.0',filter= filter110)
HOQCutil::cap.out(
x110 <- handle_body(my_body110,base_url,to_sf = T,returnres = F,verbose=T),
se =300,
width=95 )
#> -> POST /wijkenbuurten2019/wfs HTTP/1.1
#> -> Host: geodata.nationaalgeoregister.nl
#> -> User-Agent: libcurl/7.64.1 r-curl/4.3 httr/1.4.2
#> -> Accept-Encoding: deflate, gzip
#> -> Accept: application/json, text/xml, application/xml, */*
#> -> Content-Type: text/xml
#> -> Content-Length: 1321
#> ->
#> >> <?xml version="1.0" encoding="ISO-8859-1"?>
#> >> <wfs:GetFeature
#> >> service="WFS" version="1.1.0" outputFormat="application/json"
#> >> xmlns:wijkenbuurten2019="http://wijkenbuurten2019.geonovum.nl"
#> >> xmlns:wfs="http://www.opengis.net/wfs"
#> >> xmlns:ogc="http://www.opengis.net/ogc"
#> >> xmlns:gml="http://www.opengis.net/gml"
#> >> >
#> >> <wfs:Query typeName="wijkenbuurten2019:cbs_buurten_2019" srsName="EPSG:28992">
#> >> <PropertyName>wijkenbuurten2019:omgevingsadressendichtheid</PropertyName>
#> >> <PropertyName>wijkenbuurten2019:aantal_inwoners</PropertyName>
#> >> <PropertyName>wijkenbuurten2019:wijknaam</PropertyName>
#> >> <PropertyName>wijkenbuurten2019:buurtnaam</PropertyName>
#> >> <PropertyName>wijkenbuurten2019:geom</PropertyName>
#> >> <ogc:SortBy>
#> >> <ogc:SortProperty><ogc:PropertyName>wijkenbuurten2019:wijknaam</ogc:PropertyName><ogc:SortO
#> rder>DESC</ogc:SortOrder></ogc:SortProperty>
#> >> <ogc:SortProperty><ogc:PropertyName>wijkenbuurten2019:aantal_inwoners</ogc:PropertyName><og
#> c:SortOrder>ASC</ogc:SortOrder></ogc:SortProperty>
#> >> </ogc:SortBy>
#> >> <ogc:Filter>
#> >> <ogc:DWithin>
#> >> <ogc:PropertyName >wijkenbuurten2019:geom</ogc:PropertyName>
#> >> <gml:Point srsName="EPSG:28992">
#> >> <gml:coordinates cs="," decimal="." ts=" ">119400,480254</gml:coordinates>
#> >> </gml:Point>
#> >> <ogc:Distance ogc:units='meters'>500.0</ogc:Distance>
#> >> </ogc:DWithin>
#> >> </ogc:Filter>
#> >> </wfs:Query>
#> >> </wfs:GetFeature>
#> <- HTTP/1.1 200
#> <- Date: Tue, 29 Sep 2020 13:57:35 GMT
#> <- Content-Disposition: inline; filename=geoserver-GetFeature.application
#> <- Content-Type: application/json;charset=UTF-8
#> <- X-Cnection: close
#> <- Access-Control-Allow-Origin: *
#> <- Access-Control-Allow-Methods: POST, GET, OPTIONS, HEAD
#> <- Access-Control-Max-Age: 1000
#> <- Access-Control-Allow-Headers: SOAPAction,X-Requested-With,Content-Type,Origin,Authorization
#> ,Accept
#> <- X-Cnection: close
#> <- Transfer-Encoding: chunked
#> <-
#> Reading layer `OGRGeoJSON' from data source `{"type":"FeatureCollection","totalFeatures":5,"fe
#> atures":[{"type":"Feature","id":"cbs_buurten_2019.4024","geometry":{"type":"MultiPolygon","coor
#> dinates":[[[[118819.686,479910.0661],[118819.8686,479911.354],[118820.0718,479912.6387],[118820
#> .296,479913.9201 ...
#> Simple feature collection with 5 features and 5 fields
#> geometry type: MULTIPOLYGON
#> dimension: XY
#> bbox: xmin: 118787.6 ymin: 479085.1 xmax: 119918.3 ymax: 481680.6
#> projected CRS: Amersfoort / RD New
knitr::kable(x110) # show geom column
| id | buurtnaam | wijknaam | omgevingsadressendichtheid | aantal_inwoners | geometry |
|---|---|---|---|---|---|
| cbs_buurten_2019.4024 | Stadshart | Stadshart | 3630 | 3165 | MULTIPOLYGON (((118819.7 47… |
| cbs_buurten_2019.4017 | Randwijck Oost | Randwijck | 3048 | 3720 | MULTIPOLYGON (((119679.6 48… |
| cbs_buurten_2019.4022 | Vredeveldbuurt | Elsrijk | 3303 | 1695 | MULTIPOLYGON (((119282.5 48… |
| cbs_buurten_2019.4021 | Kruiskerkbuurt | Elsrijk | 3343 | 2400 | MULTIPOLYGON (((118817.5 48… |
| cbs_buurten_2019.4023 | Elsrijk Oost | Elsrijk | 3658 | 3825 | MULTIPOLYGON (((119903.7 48… |
In the second example we omit all intermediate output and only show the data.frame in the sf object without the coordinates.
my_body200 <- getfeature_body('2.0.0',filter= filter200)
HOQCutil::cap.out(
x200 <- handle_body(my_body200,base_url,to_sf = T,returnres = F,verbose=F),
se =300,
width=95 )
knitr::kable(st_drop_geometry(x200)) # drop geom column
| id | buurtnaam | wijknaam | omgevingsadressendichtheid | aantal_inwoners |
|---|---|---|---|---|
| cbs_buurten_2019.4024 | Stadshart | Stadshart | 3630 | 3165 |
| cbs_buurten_2019.4017 | Randwijck Oost | Randwijck | 3048 | 3720 |
| cbs_buurten_2019.4022 | Vredeveldbuurt | Elsrijk | 3303 | 1695 |
| cbs_buurten_2019.4021 | Kruiskerkbuurt | Elsrijk | 3343 | 2400 |
| cbs_buurten_2019.4023 | Elsrijk Oost | Elsrijk | 3658 | 3825 |
Session Info
This document was produced on 29Sep2020 with the following R environment:
#> R version 4.0.2 (2020-06-22)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows 10 x64 (build 18363)
#>
#> Matrix products: default
#>
#> locale:
#> [1] LC_COLLATE=English_United States.1252
#> [2] LC_CTYPE=English_United States.1252
#> [3] LC_MONETARY=English_United States.1252
#> [4] LC_NUMERIC=C
#> [5] LC_TIME=English_United States.1252
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] tibble_3.0.3 jsonlite_1.7.1 glue_1.4.2 purrr_0.3.4 xml2_1.3.2
#> [6] httr_1.4.2 sf_0.9-6
#>
#> loaded via a namespace (and not attached):
#> [1] Rcpp_1.0.5 highr_0.8 compiler_4.0.2 pillar_1.4.3
#> [5] class_7.3-17 tools_4.0.2 digest_0.6.25 evaluate_0.14
#> [9] lifecycle_0.2.0 pkgconfig_2.0.3 rlang_0.4.7 DBI_1.1.0
#> [13] curl_4.3 xfun_0.17 e1071_1.7-3 dplyr_1.0.2
#> [17] stringr_1.4.0 knitr_1.30 generics_0.0.2 vctrs_0.3.4
#> [21] classInt_0.4-3 grid_4.0.2 tidyselect_1.1.0 R6_2.4.1
#> [25] rmarkdown_2.3 magrittr_1.5 htmltools_0.5.0 ellipsis_0.3.1
#> [29] units_0.6-6 KernSmooth_2.23-17 stringi_1.4.6 HOQCutil_0.1.24
#> [33] crayon_1.3.4