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
Post a Comment