CODEWITHSUNDEEP
rgeospatialr-raster
Larusson
Larusson

Reputation: 267

plot LAT/LON coordinates on geotiff in R

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)

Upvotes: 5

Views: 2425

Answers (2)

Robert Hijmans
Robert Hijmans

Reputation: 47601

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

Upvotes: 5

Geo-sp
Geo-sp

Reputation: 1704

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)

Upvotes: 1

Related Questions