Search code examples
rmapsgeospatialprojection

Assigning spatial coordinates to an array in R

I have downloaded a text file of data from the following link: http://radon.unibas.ch/files/us_rn_50km.zip

After unzipping I use the following lines of code to plot up the data: # load libraries library(fields) 

# function to rotate a matrix (and transpose)
rotate <- function(x) t(apply(x, 2, rev))

# read data
data <- as.matrix(read.table("~/Downloads/us_rn_50km.txt", skip=6))
data[data<=0] <- NA

# rotate data
data <- rotate(data)

# plot data
mean.rn <- mean(data, na.rm=T)
image.plot(data, main=paste("Mean Rn emissions =", sprintf("%.3f", mean.rn)) )

This all looks OK, but I want to be able to plot the data on a lat-long grid. I think I need to convert this array into an sp class object but I don't know how. I know the following (from the web site): "The projection used to project the latitude and longitude coordinates is that used for the Decade of North American Geology (DNAG) maps. The projection type is Spherical Transverse Mercator with a base latitude of zero degrees and a reference longitude of 100 degrees W. The scale factor used is 0.926 with no false easting or northing. The longitude-latitude datum is NAD27 and the units of the xy-coordinates are in meters. The ellipsoid used is Clarke 1866. The resolution of the map is 50x50km". But don't know what to do with this data. I tried:

proj4string(data)=CRS("+init=epsg:4267")
data.sp <- SpatialPoints(data, CRS("+proj=longlat+ellps=clrk66+datum=NAD27") )

But had various problems (with NA's) and fundamentally I think that the data isn't in the right format.. I think that the SpatialPoints function wants a data on location (in 2-D) and a third array of values associated with those locations (x,y,z data - I guess my problem is working out the x and the y's from my data!)

Any help greatly appreciated!

Thanks, Alex

Solution

  • The file in question is an ASCII raster grid. Coordinates are implicit in this format; a header describes the position of the (usually) lower left corner, as well as the grid dimensions and resolution. After this header section, values separated by white space describe how the variable varies across the grid, with values given in row-major order. Open it in a text editor if you're interested.

    You can import such files to R with the fantastic raster package, as follows:

    download.file('http://radon.unibas.ch/files/us_rn_50km.zip', 
          destfile={f <- tempfile()})
unzip(f, exdir=tempdir())
r <- raster(file.path(tempdir(), 'us_rn_50km.txt'))

    You can plot it immediately, without assigning the projection:

    no_crs

    If you didn't want to transform it to another CRS, you wouldn't necessarily need to assign the current coordinate system. But since you do want to transform it to WGS84 (geographic), you need to first assign the CRS:

    proj4string(r) <- '+proj=tmerc +lon_0=-100 +lat_0=0 +k_0=0.926 +datum=NAD27 +ellps=clrk66 +a=6378137 +b=6378137 +units=m +no_defs'

    Unfortunately I'm not entirely sure whether this proj4string correctly reflects the info given at the website that provided the data (it would be great if they actually provided the definition in a standard format).

    After assigning the CRS, you can project the raster with projectRaster:

    r.wgs84 <- projectRaster(r, crs='+init=epsg:4326')

    And if you want, write it out to a raster format of your choice, e.g.:

    writeRaster(r.wgs84, filename='whatever.tif')