2

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?

2
  • What's your desired lat and long? Have you looked at the crop function in the raster package? You might need to define an extent using the extent function but you'll have to convert your lat-lon to the coordinates of the raster. – Spacedman Sep 5 '19 at 20:03
  • @Spacedman : I have updated the question to include a workable example. Thanks for pointing me in the right direction. – JForsythe Sep 12 '19 at 17:42
2

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)

enter image description here

1
  • Fantastic. Thank you this worked perfectly for my project. I appreciate all of the background information you have provided as well. – JForsythe Sep 16 '19 at 14:38
0

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')

enter image description here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.