Say that i have two "data sources";

  1. A global map of average wind speeds,
  2. A shape file (or similar information) of all country borders in the world.

How do i combine these two sources and extract the average wind speeds of a country?

The wind speed data is structured such that each lat-lon combination has an average wind speed:

Lat, Lon, AvgSpeed
-90, -180, 4.54
-90, -179, 4.55
-90, -178, 4.57

The output should be somewhat similar to the input, but with "an extra column" for country:

Lat, Lon, AvgSpeed, Country
-90, -180, 4.54, Country 1
-90, -179, 4.55, Country 1
-90, -178, 4.57, Country 2

I have Grass GIS, R and Python installed - but your answer need not be restricted to these software packages/programming languages.

  • Can we use postgis in the answer? – nickves Oct 14 '12 at 19:29

You first need to import the CSV file with v.in.ascii to create a vector points map. Add a new column "Country" of varchar(25) with v.db.addcol. Then simply populate the new column with v.what.vect (see also example in that manual page).

  • This is great - but not exactly what i was looking for. This gives me and min, max, average, etc of each country. But what i would like is all the AvgSpeed of each lat-lon inside a country. The output should be similar to the input but with an extra column with countries (i updated the question) – bonna Oct 14 '12 at 17:07
  • Ah, I read to fast, I have rewritten my answer and moved the previous example into the GRASS Wiki at grass.osgeo.org/wiki/Zonal_statistics – markusN Oct 14 '12 at 19:36

In R, the gIntersection() tool from the rgeos package should intersect and match the data, but doesn't seem to be working (gives only the shapes, not the dataframes (see @ari-b-friedman's question).

QGIS will work for you if you use the Vector|Geoprocessing|Intersect tool, giving you a set of all the points with the overlapping country data merged; this is odd as both R and QGIS should be using the same method.

Using R, you can extract a single country's data with overlay():

# Load the maptools and a shapefile of all countries
world <- readShapePoly("~/workspace/World_countries.shp")

# Select a single country from the names:
kenya <- world[world$CNTRY_NAME == "Kenya",]

Kenya outline

# Create a wind dataset and make it SpatialPoints
wind <- data.frame(lon = rep(30:45, 11), lat = rep(-5:5, each = 16), 
        AvgSpeed = rnorm(11*16, mean = 5, sd = 2))
coordinates(wind) <- ~lon+lat
plot(wind, add = T, pch = 20, cex = 0.4, col = 'red')

enter image description here

# Extract the vector of overlapping points
kenya.wind.over <- overlay(wind, kenya)

# Extract just those points from the wind data
kenya.wind <- wind[!is.na(kenya.wind.over),]

# And plot to see if it looks correct
plot(wind, col = 'gray', add = T, cex = 0.2)
plot(kenya.wind, col = 'red', cex = kenya.wind$AvgSpeed/10, add = T, pch = 20)

enter image description here

You could loop over the countries data frame (for country in CNTRY_NAME) and cbind the resulting rows, but seems much easier just to do it in QGIS.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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