Month: March 2012

GIS: Raster cell centroids as vector points with GDAL and R

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); }

Ephemeral ladybug trace fossil

I just came across this image of a ladybug making tracks in condensation on a window. As far as I can tell, this is an original photo by user 77dlsd at reddit. Fascinating! This might be a decent way of capturing different kinds of arthropod trackways, although you may have to be something of a Snowflake Bentley to frost a glass plate, drop a critter, catch a critter, and photograph the plate without losing anything.