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 am very new at using GIS data and only modestly experienced with R. I've been reading about how to analyze spatial data using the spatial-analyst.net PDF book, so I'm not completely lost, but I thought I could describe my problem and people might suggest ideas.

I have a dataset with about 2000 measurements at different lat/long coordinates, although I will probably subdivide this dataset as the data was collected over 3 years and conditions changed over time. Let's call the variable being measured "IP."

I want to create a map of IP in the full area under question using Kriging or some other interpolation method on the sample data. Then I want to create a histogram measuring the amount of land in various IP buckets. I will also need to create a histogram which shows the number of samples in each bucket (note a sample could have a higher or lower actual IP than what kriging predicts for its land).

I follow how to load the data into a SpatialPointsDataFrame and run a kriging analysis, where I'm having trouble is how to convert that data into a gridded dataframe so I can do the histogram analysis. Any suggestions for converting points into grids? I figure this is something that is easy if you know what you're doing, but unfortunately I don't!

Any other suggestions would be much appreciated.

Thanks,

Zach

share|improve this question

1 Answer 1

up vote 5 down vote accepted

You're right ... it is pretty easy! The "raster" package has some pretty straightforward ways of dealing with creating and manipulating rasters.

library(maptools)
library(raster)

# Load your point shapefile (with IP values in an IP field):
pts <- readShapePoints("pts.shp")

# Create a raster, give it the same extent as the points
# and define rows and columns:

rast <- raster()
extent(rast) <- extent(pts) # this might be unnecessary
ncol(rast) <- 20 # this is one way of assigning cell size / resolution
nrow(rast) <- 20

# And then ... rasterize it! This creates a grid version 
# of your points using the cells of rast, values from the IP field:
rast2 <- rasterize(pts, rast, pts$IP, fun=mean) 

You can assign grid size and resolution in a number of ways - have a good look at the raster package documentation.

The values of the raster cells from rasterize can be calculated with a function - 'mean' in the example above. Make sure you put this in: otherwise it just uses the value of IP from the last point it comes across!


From a CSV:

pts <- read.csv("IP.csv")
coordinates(pts) <- ~lat+lon
rast <- raster(ncol = 10, nrow = 10)
extent(rast) <- extent(pts)
rasterize(pts, rast, pts$IP, fun = mean)
share|improve this answer
    
Hey this is very useful, but how would the code look if I started with the points in a simple CSV with lat/longs rather than a shapefile? The columns in the CSV would be IP, Lat, Long, etc, etc, etc. –  user1080253 May 2 '12 at 17:03
    
You did indicate that you've already loaded the data into a SpatialPointsDataFrame ... which is exactly what pts is in my example above. Just run the code on your SpatialPointsDataFrame object! –  Simbamangu May 2 '12 at 18:33
2  
This reply, although excellent, does not seem to address what is needed. (It appears to offer a solution to gis.stackexchange.com/questions/20018 instead.) The challenge is to interpolate 2000 or so points, not just assign their values to raster cells. Given that the OP claims already to have "run a kriging analysis," this question comes down to extracting the values of a raster (say, r) for using in a hist-like procedure, which is simply a matter of an expression like hist(getValues(r)). –  whuber May 2 '12 at 19:46
    
@whuber - Looks like OP asks "where I'm having trouble is how to convert that data into a gridded dataframe so I can do the histogram analysis ... any suggestions for converting points into grids" as the actual question, and knows how to make a SpatialPointsDataFrame and run the kriging. But you're right, it seems to be a duplicate of 20018 (except for gridded input). –  Simbamangu May 3 '12 at 3:37
    
Apologies, @user1080253 ... I read 'grid' as 'raster' which isn't correct, and not helpful for kriging; see here for a better idea on creating a regular grid and interpolating your data to that grid. –  Simbamangu May 3 '12 at 3:51

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.