[GIS] Converting data.frame with no lat-long info to raster using R

rraster

I need to read in a .txt file into R and later do some spatial analyses on it.
I’m working in R. 3.0.2 on a Windows server.
I load the following libraries (probably more than I need):

library(raster)
library(rgdal)
library(rgeos)
library(sp)
I then read in the text file:
v2<-read.delim("E:\GIT\GEO_04_Research\03_Fire\03_FireIASI\Data\fire\GigGFEDRedownload\GFED3.1_199902_BA.txt", header=FALSE, sep = "", dec = ".", strip.white=TRUE)
This works fine and
str(v2) #tells me 'data.frame': 360 obs. of 720 variables.

However, there is no xy information. The documentation on this data layer states this data runs from a top-left of 179.75W, 89.75N to the lower right corner at 179.75E, 89.75S, with each reading in this regular ‘grid’ occurring at 0.5 degrees.

So I created a raster to match these specifications:

gfed<-raster(ncol=720,nrow=360,xmn=-89.75,xmx=89.75,ymn=-179.75,ymx=179.75)
gfed… # tells me:
class : RasterLayer
dimensions : 360, 720, 259200 (nrow, ncol, ncell)
resolution : 0.2493056, 0 (x, y)
extent : -89.75, 89.75, -179.75, 179.75 (xmin, xmax, ymin, ymax)
coord. ref. : NA

But now how do I get the xy information from the ‘gfed’ raster onto the values in the data ‘v2’?

I tried ‘stack’ hoping it would take the ‘square’ data.frame onto to the gfed raster, but I get 0 layers. I’m new to coding and I suspect that this requires some loop that reads row 1 from 1-720 and then goes to the next row until it reads all 360 rows.

I think part of the problem is that the data.frame ‘v2’ doesn’t have any spatial information associated with it, so SpatialPixels and SpatialDataFrame options don’t work – or at least, I can’t get them to work.

I just can’t get passed this point.

I'm almost there with your solution of using 'scan' into a newly created raster. It is a bit embarrassing to admit – I understood that my co-author on this was sending ‘world-wide’ files at half-degree scale, and previous data have been just that, and so I accepted that even though I looked at the data and saw co-ordinates that didn’t start at -180, +90. The study area is actually southern Africa (approximately xmn=13, xmx=36, ymn=-36, ymx=-18) and it is at a quarter-degree (0.25) scale!
I thought I’d be very clever and apply your solution (make grid and then scan the data in to the new grid) with these southern Africa parameters. I couldn’t get it to work and after spending time staring at the data, I realised two things:

  1. The points don’t start neatly at round numbers. I applied shifts separately to the x and y.
    1. Then I realised that the data that has been sent to me is an irregular grid. It is not consistently 0.25-degrees. It is in general has a 0.3-degree sized grid cell followed by a 0.2- degree sized grid cell followed by another 0.3- degree sized grid cell.

Any idea how I read/scan data from an irregular ‘grid’ into a nice neat 0.25 grid that I have made. I have now made

southernA<-raster(nrow=72, ncol=92, xmn=13, xmx=36, ymn=-36, ymx=-18)

I guess I would like R to ‘resample’ the ‘nearest-neighbour’ grid from the irregular grid into my newly created, regular grid called southernA.

Someone asked what the data is all about: we are looking to see how much biomass burning contributes to CO in the atmosphere. I used to do a lot of fire ecology work for a conservation agency that I worked for. I still dabble in the field. I’d be happy to post a copy of the irregular grid data if that would help find the solution.

Best Answer

So I've created a file that reads into a data frame in the same way as yours:

> str(v2)
'data.frame':   360 obs. of  720 variables:

BUT data.frame isn't really the right thing here. Its really meant for record-oriented data, where each row is a record and each column is a potentially different variable for that record (eg each row is a person, the columns are name, age, height, etc).

So you really only need to scan the data in as one long vector and feed it to a raster.

Step 1, define an empty raster of the right size and shape (note I'm assuming the raster covers the whole world, so the limits are not the cell centres):

> m2=raster(nrow=360,ncol=720,xmn=-180,xmx=180,ymn=-90,ymx=90)

Step 2, read numeric values into the raster data slot:

> m2[]=scan("d.txt",what=1)
Read 259200 items

And give it a projection if needed:

> projection(m2)="+init=epsg:4326"
> plot(m2)

If you want to check that the resolution and the cell centres are as expected, use these functions:

> res(m2)
[1] 0.5 0.5
> xFromCol(m2,1:10)
 [1] -179.75 -179.25 -178.75 -178.25 -177.75 -177.25 -176.75 -176.25 -175.75
[10] -175.25
> yFromRow(m2,1:10)
 [1] 89.75 89.25 88.75 88.25 87.75 87.25 86.75 86.25 85.75 85.25

which shows the resolution is half a degree and the cell centres (or at least the first 10) are at those specified coordinates.

Related Question