I've used the following code to make a test raster that prints like yours (but is about 1/100000000 of the size):
testlat = raster(matrix(1:100,10,10))
extent(testlat) = c(xmin=13000,xmax=278000,ymin=307000,ymax=615000)
crs(testlat) = "+proj=longlat +datum=WGS84 +no_defs"
> testlat
class : RasterLayer
dimensions : 10, 10, 100 (nrow, ncol, ncell)
resolution : 26500, 30800 (x, y)
extent : 13000, 278000, 307000, 615000 (xmin, xmax, ymin, ymax)
crs : +proj=longlat +datum=WGS84 +no_defs
source : memory
names : layer
values : 1, 100 (min, max)
Now that crs is wrong - the TIFF you got this from must have the wrong crs encoded in it. Numbers that big (in the extent) can't be lat-long (which must be degrees).
So you can tell R that the CRS isn't that:
> crs(testlat) = "+init=epsg:28992"
Warning message:
In showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj = prefer_proj) :
Discarded datum Amersfoort in Proj4 definition
You can probably ignore the warning (if you get it).
This gives you something that has the right CRS, so now anything you do with it should also get the right CRS. But note that "extents" in the raster
package don't have coordinate systems, they are just purely four numbers. So if you want to convert the extent of the raster to a polygon then you have to reassign the CRS:
p=as(testlat@extent, "SpatialPolygons")
crs(p) = crs(testlat)
pts=spsample(p, n=500, "random")
pts$value=runif(500, min=28, max = 51)
which gives you an sp
Spatial Points object with the correct coordinate system:
> pts
class : SpatialPointsDataFrame
features : 500
extent : 13631.68, 277753.9, 307923.2, 614785 (xmin, xmax, ymin, ymax)
crs : +proj=sterea +lat_0=52.1561605555556 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs
variables : 1
names : value
min values : 28.0383359261323
max values : 50.9162634953391
If you want that in lat-long then you can use the rgdal
package:
> ptsLL = spTransform(pts, "+init=epsg:4326")
> ptsLL
class : SpatialPointsDataFrame
features : 500
extent : 3.277139, 7.174634, 50.75612, 53.51738 (xmin, xmax, ymin, ymax)
crs : +proj=longlat +datum=WGS84 +no_defs
variables : 1
names : value
min values : 28.0383359261323
max values : 50.9162634953391
In your question you convert your data frame down to a two-column X,Y data frame and try and use project
. The error message is because you are using inverse=
instead of inv=
, and would also fail because project
needs a matrix argument and not a data frame. So this works:
> dataframe=as.data.frame(pts)
> dataframe$value = NULL
> pj <- project(as.matrix(dataframe), proj4, inv = TRUE)
> head(pj)
x y
[1,] 7.151821 51.43042
[2,] 5.089141 52.88179
[3,] 5.946595 51.90404
But you can directly get the coordinates as a 2-column matrix using the coordinates
function, so this is a lot neater:
pj <- project(coordinates(pts), proj4, inv=TRUE)
but ideally you should probably try and keep your data in spatial data objects unless you really need to get the raw numbers out to do something that can't be done with the spatial functions.
Also you should ideally look at doing as much as possible using the sf
package and its (non-compatible but easily convertible) spatial object classes.
crs(testlat) = "+init=epsg:28992"
thenpts
should also have that coordinate system. Then usespTransform(pts, "+init=epsg:4326")
should convert to a lat-long set of points. Read up on usingsp
andsf
in R - try: geocompr.robinlovelace.net