Short intro: I’m using raster data for part of my dissertation, but I’m storing the paleoenvironment data in a database as the latitude and longitude of the center of each grid cell. This lets me query my environment localities, add new ones, and make specific datasets.
Sometimes I take existing maps of other people’s environmental interpretations and digitize them, leaving me with vector polygon data. I need to convert the polygons into the same kind of data I’ve been storing. This is surprisingly hard in both QGIS and GDAL alone, so I finally figured out how to do it in R today courtesy of Simbamangu over at StackExchange.
See the attached images for how things look at each stage.
Below are the general terminal commands, first in GDAL and then R. I will try to add more rationalization and explanation later. My operating system is Mac OS X but these should be modifiable.
GDAL in your terminal
gdal_rasterize -at -where 'freshwater = 1' -a freshwater -ts 360 180 -te -180 -90 180 90 -l "source_layer_name_in_shapefile" "shapefile.shp" "output_raster.tif"
Produces a geoTiff of your input polygons where the attribute ‘freshwater’ is equal to ‘1’. For more detail I invite you to read the GDAL documentation, it’s quite good.
R
To run through a single file:
# load necessary libraries
library(maptools);
library(raster);
library(rgdal);
# Load the raster from a file
r <- raster("rastername.tif")
# Convert to spatial points
p <- as(r, "SpatialPointsDataFrame")
# Set name of column of raster values for later use
names(p)[1]<-"value"
# Get only the values we want (values above zero) in that column
psubset<-subset(p,p$value>0);
# Save as a shapefile (extension appended automatically)
writeSpatialShape(psubset, "outputfile");
You can also loop through a directory:
# load necessary libraries
library(maptools);
library(raster);
library(rgdal);
# Set working directory
setwd("/my/working/directory/");
# list all .tif files in working directory
for (file in dir(pattern="*.tif")) {
r <- raster(file);
p <- as(r, "SpatialPointsDataFrame");
# Set name of column of raster values
names(p)[1]<-"value";
# Get only the values we want (values above zero)
psubset<-subset(p,p$value>0);
# Make a new filename
oldname<-substr(file, start = 1, stop = nchar(file) -4);
newname<-paste(oldname,"pts",sep=" ");
# Make a shapefile
writeSpatialShape(psubset, newname) ;
# Make a CSV file
# Convert to data frame so you can drop extra columns
psubsetdf<-as.data.frame(psubset)
write.table(psubsetdf[,c(2,3)],file=paste(newname,".csv",sep=""),sep=",", row.names=FALSE);
}