Geospatial data in R: I have a matrix of latitudes, a matrix of longitudes, and a matrix of values-how combine?

I have three matrices: the first one is a matrix of longitude points, the second is a matrix of latitude points, and the third is a matrix of air quality values, measured at each latitude and longitude. In R, I would like to combine all of these into a raster. Following the discussion here, and using the madeup data below, I thought I would be able to combine the matrices in this way:

mat.lat=matrix(rep(41:50,10),10)
mat.long=matrix(rep(91:100,10),10)
mat.aq=matrix(rnorm(100),10)

r=raster(mat.aq,
    xmn=min(mat.long),
    xmx=max(mat.long),
    ymn=min(mat.lat),
    ymx=max(mat.lat)
)

However, when I plot the data, they are a little bit off. I suspect it's because my data is actually more complicated than the one in the madeup data, and the actual latitude and longitude grid aren't evenly spaced, and evenly spaced is what the raster command expects.

But I feel like there has to be a smarter way to combine the matrices than just plowing along the minimum and maximum latitude and longitude values. I actually have the values of the lat/lon, so I don't want to interpolate them. But after a lot of searching, I can't figure out how. This question came close, the poster found a solution to their problem that doesn't solve mine.

How do you layer together an two matrices of latitudes and longitudes with a matrix of data points into a raster (or other spatial data frame)?

For background, I received two netcdf files from a colleague: one with a grid of latitude and longitude points that cover our area of interest. In this first netcdf file, the latitude and longitude are variables, not dimensions. The second netcdf file contains the data we're interested in, but no latitude or longitude. I would like to map the data in these files in R, but the standard functions for reading on netcdf files assume that the dimension of the netcdf file has the latitude and longitude already. Following the discussion here, I figured that I could extract the three arrays of information, and layer them back together, but I'm having a tough time doing that. REVISION I found a way to do it, although I'm sure it's a horribly, embarrassingly inefficient way to skin the cat, I share in case someone else google-stumbles onto this page.

pts=cbind(lon=as.vector(mat.lon),
        lat=as.vector(mat.lat),
        aq=as.vector(mat.aq))
    inter1=data.frame(pts)
    inter1 <- cbind(inter1, 
        cat = rep(1L, nrow(inter1)),
        stringsAsFactors = FALSE)

    #convert to spatial points
    coordinates(inter1) = ~lon + lat
    proj4string(inter1)<-CRS("+init=epsg:4269")

    r.grid <- raster(inter1,crs=CRS("+init=epsg:4269"),
        nrows=315,
        ncols=288)
    r=rasterize(x=pts[,1:2],y=r,field=pts[,3])

Solution

I would start with using the dimensions that are provided. Presumably coordinates in a different spatial reference system? Once you have these as a raster, you can use projectRaster to transform to lon/lat.

If the data are not on a regular grid, treat the data as points. You could do this

mlat=matrix(rep(41:50,10),10)
mlong=matrix(rep(91:100,10),10, byrow=TRUE)
maq=matrix(rnorm(100),10)

pts <- cbind(lat=as.vector(mlat), lon=as.vector(mlong), aq=as.vector(maq))

If the data were on a regular grid (as they are in the example data), you could do

library(raster)
r <- rasterFromXYZ(pts)
plot(r)

If they are not, and you do want them on a regular grid, you need interpolate. See ?raster::interpolate