Convert matrix to raster with assigning LCC projection

Question

I would like to convert data-only matrix (not a xyz type of table) to raster layer. Each column and row represents longitude and latitude.

However, while rasterizing the matrix, I'm struggling with assigning CRS on it.

My data should be projected in LCC projection. And all the parameters that I know are below.

DX: 7500.f (meter)
DY: 7500.f (meter)
CEN_LAT: 37.99787f
CEN_LON: 127.4592f
TRUELAT1: 30.f
TRUELAT2: 60.f

So I assigned mycrs with CRS() function, and assigned it as below. I could not assign res=7500 argument in raster() because it returned an error.

mx <- matrix(rep(1, 13770), nrow=90, ncol=153) # reproducible example of mx
mycrs <- CRS("+proj=lcc +lat_1=30 +lat_2=60 +lat_0=37.99787 +lon_0=127.4592 +datum=WGS84 +units=m +no_defs")
r <- raster(mx, crs=mycrs)

class       : RasterLayer 
dimensions  : 90, 153, 13770  (nrow, ncol, ncell)
resolution  : 0.006535948, 0.01111111  (x, y)
extent      : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=lcc +lat_1=30 +lat_2=60 +lat_0=37.99787 +lon_0=127.4592 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
data source : in memory
names       : layer 
values      : 1, 1  (min, max)

Well, the console did not return any warning, but the result was awkward.

The dimensions seems ok but resolution and extent must be wrong.

The domain of the layer should be like

(123.25, 43.23) ------- (131.78, 43.17)
       :                        :
       :                        :
(123.80, 32.73) ------- (131.01, 32.68)

Is there any other method to convert mx into raster with the given CRS parameters?

Answer 1

You need to know the extent in terms of the coordinate reference system (crs), llc in this case. The data source should provide that for the data to be usable. You have the extent in lon/lat, but not in llc. We can estimate it:

library(raster)         
e <- extent(123.25, 131.78, 32.68, 43.23)
p <- as(e, "SpatialPolygons")
crs(p) <- "+proj=longlat +datum=WGS84"

library(rgdal)
llc <- "+proj=lcc +lat_1=30 +lat_2=60 +lat_0=37.99787 +lon_0=127.4592 +datum=WGS84"
pp <- spTransform(p, llc)

And create a RasterLayer

r <- raster(pp)

Now set the resolution (from what you provide I understand that it is 7500 m)

res(r) <- 7500
r
#class      : RasterLayer 
#dimensions : 152, 105, 15960  (nrow, ncol, ncell)
#resolution : 7500, 7500  (x, y)
#extent     : -390412.6, 397087.4, -567400.8, 572599.2  (xmin, xmax, ymin, ymax)
#crs        : +proj=lcc +lat_1=30 +lat_2=60 +lat_0=37.99787 +lon_0=127.4592 +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 

This does not match the number of rows and columns that you provided (but where did you get these?).

The extent could be approximately correct, but that it is not guaranteed at all. The data provider should really give it to you to be exact. If you actually know the number of rows and columns, then it is obviously wrong.