2

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.

bonna
  • 121
  • 4

2 Answers2

4

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).

markusN
  • 12,976
  • 31
  • 48
  • 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 http://grass.osgeo.org/wiki/Zonal_statistics – markusN Oct 14 '12 at 19:36
4

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
library(maptools)
world <- readShapePoly("~/workspace/World_countries.shp")

# Select a single country from the names:
kenya <- world[world$CNTRY_NAME == "Kenya",]
plot(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(kenya)
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.

Simbamangu
  • 14,773
  • 6
  • 59
  • 93