Error with eng_search()

[melissaminter16]

Hi, I have been running some code and it has been working fine until today I have tried to re-run something and I'm coming up with an error from the eng_search() function. Any idea why this might be?

local_catalog <- eng_search(tile) Error in while (jsonRespParsed$jobStatus %in% c("esriJobSubmitted", "esriJobExecuting")) { : argument is of length zero

[h-a-graham]

Hi @melissaminter16, Thanks for the issue. Yes, you're right and by the looks of things DEFRA have finally ditched the ESRI REST API!

As you can see, if you try to download the data from the new portal(https://environment.data.gov.uk/survey), the links now point to https://api.agrimetrics.co.uk/tiles/collections/survey/... rather than ESRI. This should be great news making access much easier and faster...

The bad news is, I haven't got a clue how to access the API, I can't find any documentation on this at all! I'd be very happy to update this package with the new API (I have in fact been waiting for this change before making any updates to this package).

I have just submitted a request here: https://dsp-support.atlassian.net/servicedesk/ asking for API documentation to be made avaiable but it doesn't seem to be public yet. I will update here once I here back... Perhaps you could also send a message asking for this to be made available stressing its importance from a research perspective?

In the meantime, I have no idea how we can programmatically get at this data (at least not without some quite hacky code that is unlikely to remain stable over time). If you are able to find any further info, please let me know!

Cheers, Hugh

[melissaminter16]

Thanks Hugh for getting back to me so quickly! I've added a request too so fingers crossed they can sort this soon!

Thanks Melissa

[melissaminter16]

This is the info I got from the request - is any of this useful? If not I will reply again

LiDAR Data Survey Catalogue: https://protect-eu.mimecast.com/s/eBP7C5yNyt3vM95FzahrX

Announcement of ESRI services ending on the DSP: https://protect-eu.mimecast.com/s/p1EJC68O8t7DPmjT6lMq8

Also, the LiDAR Geomatics site can also hold some useful information: Environment Agency Geomatics Hub (arcgis.com)

[h-a-graham]

Okay, so I heard back from DEFRA and actually they helped point me in the right direction. I had totally misunderstood the way the new OGC format worked. It is both really amazing and also a little tricky at the same time. Amazing because the data is crazy fast. For example, open up QGIS, add a WCS layer with the following url: https://environment.data.gov.uk/spatialdata/lidar-composite-digital-surface-model-first-return-dsm-1m/wcs?request=GetCapabilities&service=WCS&version=2.0.1 Which is from this page: https://environment.data.gov.uk/dataset/df4e3ec3-315e-48aa-aaaf-b5ae74d7b2bb

then you can render the layers in QGIS virtually instantly - it's really impressive! unlike the old WMS layers, this also includes the actual elevation data which is really amazing.

The tricky thing is that in order to make the most of this format, we need to rethink how we query it - the best thing is to use gdalwarp via probably {sf} which we already use but It may take me a while to figure it all out and I suspec tthere will be breaking changes. I'll try to find time to get this working in the next few weeks.

Cheers, Hugh

[h-a-graham]

This is the info I got from the request - is any of this useful? If not I will reply again

LiDAR Data Survey Catalogue: https://protect-eu.mimecast.com/s/eBP7C5yNyt3vM95FzahrX

Announcement of ESRI services ending on the DSP: https://protect-eu.mimecast.com/s/p1EJC68O8t7DPmjT6lMq8

Also, the LiDAR Geomatics site can also hold some useful information: Environment Agency Geomatics Hub (arcgis.com)

They must have replied at the same time thanks!

[melissaminter16]

Great thanks for the information that is really helpful :)

[h-a-graham]

Okay, so I've managed to make a bit of progess with testing. There are now two different ways to access the data. The first is for the most recent mosaic products, we can hit the WCS layer and it is crazy fast - I was able to download 10 km² of 1 m data in 130 seconds.

I've also managed to figure out the actual tile API which is equivalent to the previous ESRI one. this seems really good and allows access to all the data available from the catalog in the interactive app.

So the question now is how do we integrate this with the gblidar API? I like the current workflow of building a dataframe, filtering and then merging but not really sure yet how well this works for the WCS layers.

Do you think it would be useful or more confusing to have additional functions enabling access specifically to english composites?

I can imagine a scenario where everthing carries on working the way things are now (but with the new API), allowing access to specific tiles, and then specialised functions for getting the composite data much more quickly.

If you have any thoughts I'd be interested to know as I suspect you may well be the only person using this package at the moment :)

[melissaminter16]

That's so fast! Was this download done in R still or through the WCS layer in QGIS? For me, the current workflow of splitting a region into chunks worked well for me, basically I'm downloading lidar data across quite large areas in tiles at a time, then processing each tile to calculate tree canopy/height with the ForestTools package. This is quite intensive so tile at a time works well for this, and had this in a loop working through every tile. Yes I think your scenario sounds good! I haven't had a go at loading the WCS data in QGIS yet but will do very soon. Thanks for your help on this!

[h-a-graham]

So that download was done with R but really it's gdal wrapped in R so I'd expect similar performance with QGIS. The nice thing with QGIS is if you're just interested in viewing all the data, adding a WCS layer renders on the fly so you can look at everything without actually downloading anything.

THanks for the context of your use case - which products do you tend to use, is it mainly the composites you're interested in?

No worries - I've been wanting to get better access to these data for a long time and its very helpful to talk to others using it.

[melissaminter16]

I usually just use the DSM and DTM at 1m. Have loaded into QGIS but unfortunately can't do much with them - was hoping to clip to the areas I was missing in the mean time but looks like it won't let me!

[h-a-graham]

Okay so I've updated the ESRI API to the new OGC API. it's not ready to merge into the main branch yet but your current workflow should now work again:

install the DEV branch with:

# install.packages("pak")
pak::pkg_install("h-a-graham/gblidar@DEV")

This needs more work but any feedback you have would be very welcome.

let's leave this issue open and I'll add updates.

Cheers, Hugh

[melissaminter16]

Great thanks for your help Hugh -it works!! Life saver :)

Thanks Melissa

[h-a-graham]

Okay so I've added some functionality for using the new wcs composites to the DEV branch. If you reinstall again you should have access. There are 4 different dataets available, one of which is the vom model (which looks to be a canopy height model) this could save you some processing time.

In general though this seems to be the way to go for the english composite products. Still a bit rough round the edges but I thought it might be cool for you to test if you're interested. Here is a reprex:

library(sf)
#> Linking to GEOS 3.12.1, GDAL 3.6.4, PROJ 9.1.1; sf_use_s2() is TRUE
#> WARNING: different compile-time and runtime versions for GEOS found:
#> Linked against: 3.12.1-CAPI-1.18.1 compiled against: 3.11.1-CAPI-1.17.1
#> It is probably a good idea to reinstall sf, and maybe rgeos and rgdal too
library(terra)
#> terra 1.7.71
library(gblidar)

search_box <- st_point(c(370126.5, 538567.1)) |>
  st_buffer(1000) |>
  st_sfc() |>
  st_set_crs(27700)

options(gblidar.out_raster_type = "SpatRaster")
options(gblidar.progress = FALSE) # for this reprex for readability.
fz_dsm <- eng_composite(search_box, product = "fz_dsm")
dsm <- eng_composite(search_box, "dsm")
dtm <- eng_composite(search_box, "dtm")
vom <- eng_composite(search_box, "vom", "elevation")

dsm_hs <- eng_composite(search_box, product = "dsm", product_type = "hillshade")
dtm_hs <- eng_composite(search_box, "dtm", "hillshade")

plot_func <- function(x, label) plot(x, main = label, col = hcl.colors(150, "SunsetDark"))
par(mfrow = c(3, 2))
invisible(mapply(
  plot_func, list(fz_dsm, dsm, dtm, vom, dsm_hs, dtm_hs),
  c("fz_dsm", "dsm", "dtm", "vom", "dsm_hs", "dtm_hs")
))

par(mfrow = c(1, 1))

^{Created on 2024-02-14 with reprex v2.0.2}

[melissaminter16]

Hi Hugh,

Sorry for the delay! This works and was extremely quick for those tiles :) Is the default to download at 1m resolution? And can I check what fz_dsm is - is this first return?

Thanks! Melissa

[h-a-graham]

No worries! So this is the fz_dsm: https://environment.data.gov.uk/dataset/df4e3ec3-315e-48aa-aaaf-b5ae74d7b2bb - looks like you're right that it's the first return.

So the default resolution is 1 m but you can switch this to 2 m with the prod_res argument; changing this will load a different data source.

But there is also another argument warp_res where you can specify the desired output resolution of your raster - this argument is provided to gdalwarp so you can down or upsample as you like and specify the resmapling method with the resample argument.

work is a bit turbo at the moment but hoping to be able to add some more testing and better documentation in the not too distant future.

Cheers and glad its looking useful!