Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Large .img file processing in R (GIS)

Tags:

r

gis

I am working with U.S. National Landcover Dataset (NLCD) to classify habitat type at more than 150 sites across the northeast U.S. The dataset is very large (15GB) so I cannot upload it here, but it comes in .img format at 30m resolution. I have GPS coordinates for the center point of all the sites. I would like to be able to extract the proportion of landcover classes in a 1 square kilometer around the point. My questions are:

1) How do I upload .img files into r? 2) How do I extract the information from around the GPS coordinates as proportions of the different habitat classes?

Has anyone worked with this dataset in r before? If so I could really use the help. Cheers, Israel

like image 792
I Del Toro Avatar asked Feb 17 '23 22:02

I Del Toro


2 Answers

Use the raster package which can process the files from disk, only reading in chunks at a time.

the raster package has an extract function with a buffer argument. set your buffer to the appropriate value (1000 if your map units are metres and you want a km radius)

like image 112
mnel Avatar answered Feb 27 '23 23:02

mnel


Thanks to mnel. I have gotten the basic idea to work (code below). Now if anyone could give me a pointer on how to calculate the proportion of each category for every coordinate. The extract function gives me matrices of values for each set of coordinates. Is there a way to summarize this data?

#load in map and locality data
NLCD<-raster ('NLCD2006/NLCD2006.img')
sites<-read.csv('sites.csv', header=T)
#crop site data to just latitude and longitude
sites<-sites[,4:5]
#convert lat/lon to appropirate projection
str (sites)
coordinates(sites)  <-  c("Longitude",  "Latitude")
proj4string(sites)  <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")
sites_transformed<-spTransform(sites, CRS("+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"))

#plot the map
plot (NLCD)
#add the converted x y points
points (sites_transformed, pch=16, col="red", cex=.75)
#extract values to poionts
Landcover<-extract (NLCD, sites_transformed, buffer=1000)
like image 45
I Del Toro Avatar answered Feb 27 '23 23:02

I Del Toro