Image Service Pixel Size Cannot be Controlled with arc.raster

[Nova-Scotia]

I'm trying to produce a map using R in RStudio that uses an image service as the baselayer. I'd like to use the arcgisbinding package to clip a raster from an image service provided in ArcGIS Online.

So far, I've "cheated" - I bring the image service into ArcGIS, export the area of interest as a georeferenced .TIF, and bring it into R:

library(raster)
library(sf)
library(tidyverse)
library(RStoolbox)

basemap <- brick("MyBaseLayer.tif")

Then, I map my points onto that map, using ggplot2 and RStoolbox:

ggRGB(basemap, r = 1, g = 2, b = 3) +
  geom_sf(data = my_points, size = 5) 

Ideally, I'd like to be able to specify an extent and then use tools within R to clip and use part of a raster from an image service. Is this possible currently? I'm having trouble finding any resources on how to do this. The closest I got was how to access raster imagery already on a local machine.

What if I'm starting with a study area polygon (defined below), and I want to use it to define the extent of a basemap layer for a map I'm creating in R? For example, what if I wanted to use "Aerial Photography Image Service (Orthophoto) - 2019", an image service in ArcGIS online, as the basemap?

View imagery
https://www.arcgis.com/home/webmap/viewer.html?url=https%3A%2F%2Fimagery.dcgis.dc.gov%2Fdcgis%2Frest%2Fservices%2FOrtho%2FOrtho_2019%2FImageServer&source=sd

Image service url
https://imagery.dcgis.dc.gov/dcgis/rest/services/Ortho/Ortho_2019/ImageServer

library(raster)
library(sf)
library(tidyverse)
library(RStoolbox)

study_area <- tibble(latitude = c(38.871738, 38.869475, 38.869884),
                    longitude = c(-77.046352, -77.048959, -77.04486))

study_area <- study_area %>%
  st_as_sf(coords = c("longitude", "latitude"), 
           crs = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0") %>%
  summarise(geometry = st_combine(geometry)) %>%
  st_cast("POLYGON")

plot(study_area) # show polygon

[orhuna]

You can use the arc.raster object directly and read in the remote raster. I also used a boundary polygon from DC

library(arcgisbinding)
library(raster)

arc.check_product()

## Read in a District Polygon from D.C.
r.bound <- arc.open('https://maps2.dcgis.dc.gov/dcgis/rest/services/DCGIS_DATA/Planning_Landuse_and_Zoning_WebMercator/MapServer/31')

## Set Boundary to Image Service
r.obj <- arc.open('https://imagery.dcgis.dc.gov/dcgis/rest/services/Ortho/Ortho_2019/ImageServer')
r.arc <- arc.raster(r.obj, extent = r.bound@extent)

r.R <- as.raster(r.arc)
plotRGB(r.R)

I ran this once successfully and second time I tried to run it got a read error and currently this is what the Ortho Image Service is returning:

Once the server is back up the above should work.

[Nova-Scotia]

Thanks so much @orhuna, this is much closer!

When I ran your code (and it looks like their servers are back up) I got an error after the line r.R:
NAs produced by integer overflowError in matrix(no_data_val, nrow = nrow * ncol, ncol = length(bands)) : invalid 'nrow' value (too large or NA).

I thought maybe the raster created by the operation was too big, so I tried reducing the size of the extent object, but it still outputs the same error:

library(arcgisbinding)
library(raster)
library(sf)
library(tidyverse)
library(RStoolbox)
arc.check_product()

## Read in a District Polygon from D.C.
r.bound <- arc.data2sf(arc.select(arc.open(
'https://maps2.dcgis.dc.gov/dcgis/rest/services/DCGIS_DATA/Planning_Landuse_and_Zoning_WebMercator/MapServer/31')))
Blagdon <- r.bound %>% filter(NAME == "Blagden Alley Residential Transition Area")
plot(Blagdon)

## Set Boundary to Image Service
r.obj <- arc.open('https://imagery.dcgis.dc.gov/dcgis/rest/services/Ortho/Ortho_2019/ImageServer')
r.arc <- arc.raster(r.obj, extent = extent(Blagdon))

r.R <- as.raster(r.arc)

NAs produced by integer overflowError in matrix(no_data_val, nrow = nrow * ncol, ncol = length(bands)) : 
  invalid 'nrow' value (too large or NA)

Maybe it has to do with my version of arcgisbinding?

R version 4.0.4 (2021-02-15)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 18362)
arcgisbinding_1.0.1.244

[orhuna]

@Nova-Scotia apologies for the late response. Here is a "fix" for the problem above. It is not a sustainable fix and I will walk you through my brooding over process:

Get the Rest Definition for Image Service

library(arcgisbinding)
library(raster)
library(sf)
library(tidyverse)

arc.check_product()

## Read in a District Polygon from D.C.
r.bound <- arc.data2sf(arc.select(arc.open(
  'https://maps2.dcgis.dc.gov/dcgis/rest/services/DCGIS_DATA/Planning_Landuse_and_Zoning_WebMercator/MapServer/31')))
Blagdon <- r.bound %>% filter(NAME == "Blagden Alley Residential Transition Area")

## Set Boundary to Image Service
r.obj <- arc.open('https://imagery.dcgis.dc.gov/dcgis/rest/services/Ortho/Ortho_2019/ImageServer')
print(r.obj)

This returns:

dataset_type    : RasterDataset
path            : https://imagery.dcgis.dc.gov/dcgis/rest/services/Ortho/Ortho_2019/ImageServer 
format          : Image Service
pixel_type      : U8 (8bit)
compression_type: NA
nrow            : 300000
ncol            : 240000
extent          : xmin=389400, ymin=124200, xmax=408600, ymax=148200
WKT             : PROJCS["NAD_1983_StatePlane_Maryland_FIPS_1900",GEOGCS["GCS_...
WKID            : 26985
bands           : 4
         ncol   nrow min max     mean   stddev
Band_1 240000 300000   0 255 90.10798 55.68015
Band_2 240000 300000   0 255  92.6319 52.41673
Band_3 240000 300000   0 255 76.96512  50.9368
Band_4 240000 300000   0 255 138.3582 59.66786

Subset the Image Service with Blagdon

r.arc <- arc.raster(r.obj, extent=extent(Blagdon))
print(r.arc)

This returns

type            : Raster
pixel_type      : U8 (8bit)
nrow            : 300000
ncol            : 240000
resample_type   : NearestNeighbor
cellsize        : 0.000720728385417412, 0.000548166250002881
nodata          : NA
extent          : xmin=397746, ymin=137367.9, xmax=397919, ymax=137532.3
WKT             : PROJCS["NAD_1983_StatePlane_Maryland_FIPS_1900",GEOGCS["GCS_...
WKID            : 26985 
bands[1:4]      : Band_1, Band_2, Band_3, Band_4

Note that extent is updated by Blagdon BUT nrow and ncol are the same. It seems we cannot control the resolution of a remote imager service from our arc object and it automatically does resampling to fit nrow and ncol .

As for my "fix" (for lack of a better word)

r.arc <- arc.raster(r.obj, extent=extent(Blagdon), nrow = r.arc$nrow/10, ncol = r.arc$ncol/10)
print(r.arc)

Set the number of rows and columns explicitly. This is not a good solution since it changes the original resolution but this is the size I can fit into an R matrix.

We will be working on functionality for on-the-fly clipping so that we don't need this "fix".

[Nova-Scotia]

Hey, thanks for this. If there's a place that I can follow the progress on this functionality (on the fly clipping), let me know. This is fantastic though, I got the raster to appear (after dividing by 1000, 10 caused R to hang for me). Really appreciate it!

[orhuna]

@Nova-Scotia glad it is working with the fix above. We will track this internally but will keep this issue open and you can track our progress here.

[Nova-Scotia]

@orhuna any progress on this?

[orhuna]

@Nova-Scotia thanks for following up!

Unfortunately, I am not with the R-ArcGIS team anymore. I am adding @JosiahParry who would be a great person to follow up on the latest arc.raster functionality and whether there are other functionality in the pipeline that might help.