Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I have a data set of values over a km grid in the continental U.S. The columns are "latitude", "longitude", and "observation", e.g.:

"lat"    "lon"     "yield"
 25.567  -120.347  3.6 
 25.832  -120.400  2.6
 26.097  -120.454  3.4
 26.363  -120.508  3.1
 26.630  -120.562  4.4

or, as an R data frame:

mydata <- structure(list(lat = c(25.567, 25.832, 26.097, 26.363, 26.63), 
lon = c(-120.347, -120.4, -120.454, -120.508, -120.562), 
yield = c(3.6, 2.6, 3.4, 3.1, 4.4)), .Names = c("lat", 
"lon", "yield"), class = "data.frame", row.names = c(NA, -5L))

(the full data set can be downloaded as csv here)

The data are output from a crop model (intended to be on) a 30km x 30km grid (from Miguez et al 2012).

enter image description here

How can I convert these to a raster file with GIS - related metadata such as map projection?

Ideally the file would be a text (ASCII?) file because I would like for it to be platform and software independent.

share|improve this question
    
As CSV, this already is a "text file" in ASCII. Also, as it uses no projection at all, there may be little relevant metadata to add (datum, mostly). Could you be a little more specific about what kind of output you seek and what you intend to do with it? –  whuber Feb 8 '12 at 22:42
    
I would like to make it as easy as possible for someone to use the data to with a variety of mapping software (ArcGIS, Google Maps, Grass, R, etc.) so as to facilitate reuse, e.g. by not requiring additional conversion steps. Based on the Wikipedia page of GIS file formats, I infer 1) a "raster" file should have rownames with latitude and column names of longitude, like an image and that 2) metadata should include geographical information (location of a corner, area covered by data). –  Abe Feb 8 '12 at 23:30

2 Answers 2

up vote 21 down vote accepted

Several steps required:

  1. You say it's a regular 1km grid, but that means the lat-long aren't regular. First you need to transform it to a regular grid coordinate system so the X and Y values are regularly spaced.

    a. Read it into R as a data frame, with columns x, y, yield.

    pts = read.table("file.csv",......)
    

    b. Convert the data frame to a SpatialPointsDataFrame using the sp package and something like:

    library(sp)
    library(rgdal)
    coordinates(pts)=~x+y
    

    c. Convert to your regular km system by first telling it what CRS it is, and then spTransform to the destination.

    proj4string(pts)=CRS("+init=epsg:4326") # set it to lat-long
    pts = spTransform(pts,CRS("insert your proj4 string here"))
    

    d. Tell R that this is gridded:

    gridded(pts) = TRUE
    

    At this point you'll get an error if your coordinates don't lie on a nice regular grid.

  2. Now use the raster package to convert to a raster and set its CRS:

    r = raster(pts)
    projection(r) = CRS("insert your proj4 string here")
    
  3. Now have a look:

    plot(r)
    
  4. Now write it as a geoTIFF file using the raster package:

    writeRaster(r,"pts.tif")
    

This geoTIFF should be readable in all major GIS packages. The obvious missing piece here is the proj4 string to convert to: this will probably be some kind of UTM reference system. Hard to tell without some more data...

share|improve this answer
    
+1 Thanks for laying out the workflow. Note that the data are available at the link provided in the question: take a look. You will discover, alas, that some of your assumptions about them are incorrect. (In particular, I hunted for any documentation about the projection used to create the grid but found none. And it's a strange projection, as you can see by plotting the points.) –  whuber Feb 9 '12 at 13:38
    
Its very close to being a UTM system, but none of the ones I've tried are close enough to a regular grid for R to grid them. I'm half tempted to loop through R's entire epsg database.... –  Spacedman Feb 9 '12 at 14:23
    
That would be a real tour de force if you could discover the projection that way! The key is to find an effective and efficient criterion to determine when these 7,000+ points are close enough to lying on a regular grid (because it's possible they might not form a perfect grid in any standard projection at all). For a quick run through the database, it should be enough to compare a small number of distances, such as an east-west distance in the north part of the grid to an east-west distance in the south part. That ought to eliminate the vast majority of candidates quickly. –  whuber Feb 9 '12 at 14:32
3  
I ran through all the (default) projections supported by Mathematica 8. It found a projection in which the points really do seem to fall on a grid: Alaska State Plane (1983) Zone 10! This is a Lambert Conformal Conic projection. I believe it is EPSG 26940. If you modify this to center it approximately at longitude -106, the points form a pretty good grid. –  whuber Feb 10 '12 at 15:32
1  
Abe, do you mean to read the Web page? It was r = Import[ "https://ebi-forecast.igb.illinois.edu/bety/miscanthusyield.csv", "Data"];. You can get a quick plot of the points afterwards via data = Rest[r]; ListPlot[data[[;; , {3, 2}]]] (or ListPointPlot3D[data[[;; , {3, 2, 4}]]]). For reprojections, start with the help on GeoGridPosition, then make some intelligent guesses and cross-references to figure out what's going on :-). BTW, @Spacedman's explanation really is relevant: the metric distortion from 25 to 49 degrees equals cos(25)/cos(49) = 1.38; that's substantial. –  whuber Feb 10 '12 at 22:14

Since the question was last answered, a much easier solution exists by using the raster package's rasterFromXYZ function that encapsulates all of the steps necessary (including specification of the CRS string).

share|improve this answer
    
Apologies to the tireless @Spacedman, who has assisted me often, but I think this answer deserves to inherit the jolly green tick. –  geotheory yesterday

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.