# How to get latitude and longitude for a matrix data using specific projection in R?

I have matrix of temperature data and I would like to get lat and long for this data set in a particular projection

e.g.

PROJCS["DHDN_3_Degree_Gauss_Zone_3",GEOGCS["GCS_Deutsches_Hauptdreiecksnetz",DATUM["D_Deutsches_Hauptdreiecksnetz",SPHEROID["Bessel_1841",6377397.155,299.1528128]]

Does anyone know how do this in R?

The matrix data is like:

NCOLS 654   654
NROWS 866   866
XLLCORNER   3280414.    711633467
YLLCORNER 5237500.  62890625
CELLSIZE 1000
NODATA_VALUE -999
-999    -999    -999    -999    -999    -999    -999    -999    -999    -999 -999 15 18 12 13 20 15 16 18 21
• Does the matrix have the numbers for the coordinates in it? Could you show us a few rows? – Spacedman Nov 28 '18 at 15:45
• NCOLS 654 654 NROWS 866 866 XLLCORNER 3280414. 711633467 YLLCORNER 5237500. 62890625 CELLSIZE 1000 NODATA_VALUE -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 15 18 12 13 20 15 16 18 21 – Sasha Kosanic Nov 28 '18 at 16:05

Searching for that coordinate reference gets me here: https://epsg.io/31463

Your file appears to be an ESRI ASCII GRID header (but with the NCOLS and NROWS repeated, which is odd, but lets ignore).

If I create a similar but scaled down version:

NCOLS 35
NROWS 46
XLLCORNER 3280414. 711633467
YLLCORNER 5237500. 62890625
CELLSIZE 20000
[ 35 * 46 numbers....]

then I can read it in using the raster package:

library(raster)
r = raster("./test.grd")

and assign it that projection string:

crs = "+init=epsg:31463"
projection(r) = crs

Now r is a raster in that coordinate system:

> r
class       : RasterLayer
dimensions  : 46, 35, 1610  (nrow, ncol, ncell)
resolution  : 20000, 20000  (x, y)
extent      : 3280414, 3980414, 5237500, 6157500  (xmin, xmax, ymin, ymax)
coord. ref. : +init=epsg:31463 +proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +datum=potsdam +units=m +no_defs +ellps=bessel +towgs84=598.1,73.7,418.2,0.202,0.045,-2.455,6.7
plot(r)

Let's check this is in the right place, by transforming the extent to the WGS84 GPS lat-long coordinate system, epsg code 4326:

> projectExtent(r, "+init=epsg:4326")
class       : RasterLayer
dimensions  : 46, 35, 1610  (nrow, ncol, ncell)
resolution  : 0.3156736, 0.1834158  (x, y)
extent      : 5.523803, 16.57238, 47.1081, 55.54523  (xmin, xmax, ymin, ymax)
coord. ref. : +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0

So the lat-long of the corners of the object are about (5.52, 47.11) for the bottom left and (16.57, 55.54) at the top-right which corresponds to Germany.

However you can't simply scale the grid to a lat-long coordinate system, because a grid of 1000m squares in that Gauss-Kruger system is not a set of regular squares in degrees lat-long. You need to reproject, which will stretch your data:

rll = projectRaster(r, crs="+init=epsg:4326")
plot(rll)

The other thing you can do is get the Gauss-Kruger coordinates of each cell centre and transform those to lat-long, and get a matrix of three columns, lat, long, value. But that will not form a regular grid.

So the final process depends on precisely what you plan to do with this data. But hopefully there's enough bits here to get you on the way.