geospatial - plot LAT/LON coordinates on geotiff in R -


i have sea ice map "geotiff". eventual goal extract sea-ice concentrations @ specific lat/lon coordinates.

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

what trying load geotiff using raster() , overlay locations , using extract() function acquire values raster file @ specific location.

however lat/lon points accumulate in centre of map. go wrong? or input 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))   ### long , lat 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 reference sea-ice concentration  seaice_conc = extract(r1, spdf) 

geo-sp's solution work, sub-optimal (slow , imprecise). should (re-)project vector (points in case) data , not raster data. projecting raster data changes values, while not case vector data. projecting raster data more computationally intensive.

thus, should 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 important note made 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 spatialpointsdataframe , assign wrong coordinate reference system (crs). must assign crs matches data, "+proj=longlat +datum=wgs84". after that, can transform data crs have (using sptransform).


Comments