1

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.

  • If you don't have spatial information in the df or supporting information about which row/cell represents which location, there is no way to put this into a raster. Please post sth like head(v2) to give us a clue what this data shows. – Curlew Jan 21 '14 at 11:28
  • head(v2) will give us 6 rows of 720 numbers... – Spacedman Jan 22 '14 at 8:18
  • 1
    Sure you got x and y right? That's surely degrees, and x goes from -180 to 180... (longitude) – Spacedman Jan 22 '14 at 8:22
4

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.

  • Dear Spacedman! Thanks so much for your solution - I think it is the correct approach, but I had some glitches due to my own fault - not looking at the data thoroughly enough.. due to space, I'm going to try report back and ask for new help (pretty please) in a new 'Answer Your Question' below...Tx – Helen Jan 24 '14 at 9:44
  • If its just an adjustment to your question you should edit your question above rather than add an answer below - unless you are really answering your question! – Spacedman Jan 25 '14 at 21:54
1

You started right with the:

gfed<-raster(ncol=720,nrow=360,xmn=-89.75,xmx=89.75,ymn=-179.75,ymx=179.75)

Then you can define the projection and set the resolution right:

projection(gfed) <- "+proj=utm +zone=48 +datum=WGS84" # for example
res(gfed) <- 100

Finally, you can use:

gfed[]<-v2

to put your data in the raster.

  • 3
    A 'gotcha' here is that the xmn,xmx etc parameters to raster specify the outer edges of the cells and not the cell centre points. So when I see a raster that has a limit of '179.75W' I suspect that's probably a cell centre, and the cell extent is actually 180W to 179.5W, so ymn should be -180 so the grid covers the whole world. Similarly for the other extents. – Spacedman Jan 22 '14 at 8:15
  • 1
    I'm also not sure where you got that projection from. This looks like lat-long degrees, most likely WGS84, so "+init=epsg:4326". – Spacedman Jan 22 '14 at 8:17

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.