1

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
2

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)

enter image description here

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)

enter image description here

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.