Converting points to latitude and longitude values in R

Question

I have a .tiff file with some coordinate data associated with temperatures that I need to convert into latitude and longitude, using R. I used the following code to read the original file:

testlat<-raster("lol.tif")
p=as(testlat@extent, "SpatialPolygons")
pts=spsample(p, n=500, "random")
pts$value=runif(500, min=28, max = 51)
dataframe=as.data.frame(pts)

If I open testlat, this is what I see

If I open dataframe, this is what I see:

According to ArcGIS, the CRS name is EPSG:28992 - Amersfoort / RD New and the method is: Oblique Stereographic Alternative. I know these points all belong to The Netherlands. To try and convert the x and y into latitude and longitude, I did the following:

proj4 <-"+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +towgs84=565.417,50.3319,465.552,-0.398957,0.343988,-1.8774,4.0725 +units=m +no_defs"
### remove temperature column so we only have x and y
dataframe$value<-NULL
pj <- project(dataframe, proj4, inverse = TRUE)

Here, I get the following error:

Error in project(dataframe, proj4, inverse = TRUE) : unused argument (inverse = TRUE)

Does anyone know how I can fix this?

I have no experience with ArcGIS and very limited experience with R. I tried following the instructions on similar questions, but could not find anything that could help my case.

Answer 1

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.