Use Lat and Long to Zoom and Crop A Color GEOTIFF Based Map In R

Question

I have a strong background in R but have only done limited GIS work.

I have a color Geotiff file that I downloaded from the USGS Earthexplorer. It is a "LandsatLook Image with Geographic Reference". https://earthexplorer.usgs.gov/

Given:

download.file("https://cdn.filesend.jp/private/4DxFm0MU9rS3nl-VdLqX1CKUGUBhasI0wQL4wg_solTUoS3QizjrJRQvub0IM8dl/LC08_L1TP_011031_20190731_20190819_01_T1.tif", destfile = "myraster.tif")
LL_Limits<-list(x=c(-70.4,-69.9), y=c(42.1,41.5)) # Lat/Long of Desired Plotting Location
POI<-c(-70.2,42.05) #Point Of Interest I Want To Mark

So far in my simple example I can load raster and look at it in color with,

library(raster)
x<-brick("myraster.tif")
plotRGB(x)

But I am unsure how to proceed with utilizing CROP and EXTENT to get the output I need. I want to crop to my desired lat and long, have the image be at the full resolution of the file, and put a marker on my point of interest. For reference the output of GDALinfo is:

library(rgdal)
GDALinfo("myraster.tif")

rows        7931 
columns     7801 
bands       3 
lower left origin.x        319785 
lower left origin.y        4503885 
res.x       30 
res.y       30 
ysign       -1 
oblique.x   0 
oblique.y   0 
driver      GTiff 
projection  +proj=utm +zone=19 +datum=WGS84 +units=m +no_defs 
file        myraster.tif 
apparent band summary:
  GDType hasNoDataValue NoDataValue blockSize1 blockSize2
1   Byte          FALSE           0        256        256
2   Byte          FALSE           0        256        256
3   Byte          FALSE           0        256        256
apparent band statistics:
  Bmin Bmax Bmean Bsd
1    0  255    NA  NA
2    0  255    NA  NA
3    0  255    NA  NA
Metadata:
AREA_OR_POINT=Area 

Any suggestions?

Answer 1

In GIS operations you must care about coordinate reference system (CRS) of each file. You must learn a little bit about it before doing this kind of processes.

Landsat images are distributed in a cartesian CRS, not a geographic CRS. What does this mean? The coordinate unit of a Landsat scene is meter, not degree. GDALinfo("myraster.tif") is showing you implicitly:

lower left origin.x        319785 
lower left origin.y        4503885 
res.x       30 
res.y       30

and explicitly:

projection  +proj=utm +zone=19 +datum=WGS84 +units=m +no_defs 
file        myraster.tif

The most used geographic reference system is based in WGS 84 datum. So, it's imperative to convert coordinates to a unique CRS before crop the raster by a boundary box.

The process is to create an extent object, convert it to a polygon and reproject data to Landsat scene's CRS. The proj4 string of geographic WGS 84 CRS is in spatialreference page.

ext <- extent(c(-70.4,-69.9,41.5,42.1))

poly <- as(ext,'SpatialPolygons')

proj4string(poly) <- '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs '
poly <- spTransform(poly,proj4string(x))

Then apply crop function:

x2 <- crop(x,poly)

plotRGB(x2)

Answer 2

Thank you to @aldo_tapia. I have expanded on @aldo_tapia's answer to completely achieve my desired output the code is:

library(raster)
library(rgdal)
library(sp)

x<-brick("myraster.tif")
ext <- extent(c(-70.4,-69.9,41.5,42.1))
poly <- as(ext,'SpatialPolygons')

proj4string(poly) <- '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs '
poly <- spTransform(poly,proj4string(x))

x2 <- crop(x,poly)
plotRGB(x2)

POI<-c(-70.2,42.05)
xy <- data.frame(ID = 1, X = POI[1], Y = POI[2])
coordinates(xy) <- c("X", "Y")
proj4string(xy) <- '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs '
POI_Loc <- spTransform(xy,proj4string(x))
points(POI_Loc@coords,pch=9,col='red')