plot LAT/LON coordinates on geotiff in R

Question

I have a sea ice map as a "geotiff". The eventual goal is to extract sea-ice concentrations at specific lat/lon coordinates.

The geotiff can be found at: http://www.iup.uni-bremen.de:8084/amsr2data/asi_daygrid_swath/n6250/2015/jan/asi-AMSR2-n6250-20150101-v5.tif

What i am trying to do is load the geotiff using raster() and then overlay it by my locations and then using the extract() function to acquire values from the raster file at the specific location.

However my lat/lon points accumulate in the centre of the map. Where do I go wrong? Any help or input is greatly appreciated!

library(raster)
library(sp)

r1 = raster("test.tif")

##check plot
plot(r1)

## check projection
projection(r1)

mydf <- structure(list(longitude = rep(22,7), latitude = seq(60,90,5)),.Names = c("longitude","latitude"), class = "data.frame", row.names = c(NA, -7L)) 

### Get long and lat from data.frame.
xy <- mydf[,c(1,2)]
spdf <- SpatialPointsDataFrame(coords = xy, data = mydf,
                           proj4string = CRS("+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"))

 points(spdf, col="red", lwd=2)

 ##extract RGB values as reference for sea-ice concentration
 seaice_conc = extract(r1, spdf)

Your tiff file is not in geographic projection; you have to re-project it first. — Geo-sp, Apr 06 '16 at 14:45

score 5 · Accepted Answer · answered Apr 06 '16 at 16:20

Geo-sp's solution would work, but is sub-optimal (slow and imprecise). You should always (re-)project your vector (points in this case) data and not your raster data. Projecting raster data changes values, while this is not the case with vector data. Projecting raster data is also much more computationally intensive.

Thus, you should do something like this:

library(raster)

r <- raster("asi-AMSR2-n6250-20150101-v5.tif")
crs(r)
# CRS arguments:
# +proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 

df <- data.frame(longitude = rep(22,7), latitude = seq(60,90,5), ID=1:7)

spdf <- SpatialPointsDataFrame(coords = df[,1:2], data = df,
         proj4string = CRS("+proj=longlat +datum=WGS84"))


library(rgdal)
p <- spTransform(spdf, crs(r))       

extract(r, p)

It is important to note that you made a very common mistake:

spdf <- SpatialPointsDataFrame(coords = xy, data = mydf, proj4string =
          CRS("+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m"))

You create a SpatialPointsDataFrame and assign the wrong coordinate reference system (crs). You must assign the crs that actually matches your data, which is "+proj=longlat +datum=WGS84". After that, you can transform the data to the crs that you would like to have (using spTransform).

great! that works...thanks for the detailed explanation. Now i just struggle plotting it, as if I try to plot `plot(r)` i get the following error message `Error in colnames<-(*tmp*, value = "asi.AMSR2.n6250.20150101.v5") : length of 'dimnames' [2] not equal to array extent` — Larusson, Apr 06 '16 at 16:50

score 1 · Answer 2 · answered Apr 06 '16 at 14:51

1

You can use projectRaster to reproject your raster file. Then you can overlay the points and extract the values.

newproj <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0")
p <- projectRaster(r1, newproj)

answered Apr 06 '16 at 14:51

Geo-sp

1,704
3
19
42

plot LAT/LON coordinates on geotiff in R

2 Answers2