Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I'm trying to recreate the Gridded Mercury emissions 2000 (0.5 x 0.5 degree grid) designed by the Arctic Monitoring and Assessment Programme. The file is in ASCII format, but without any coordinates. The arithmetic to produce the grid is provided, but I don't quite understand how to create the actual grid with the information provided. For you math junkies this is probably a fairly simple problem. I have access to ArcGIS 10 to complete this work. This is the arithmetic used, and the cellcodes are on the downloadable dataset.

259200 (720 x 360) cells
Z05 Cellcode = (j*1000) + i
j = row number starting at 1 for 90S to 89.5S latitude, to 360 for 89.5N to 90N latitude
i = column number starting at 1 for 180W to 179.5W longitude, to 720 for 179.5E to 180E longitude.
(coordinates represent the center of the gridcell)
The latitude and longitude of the center of a gridcell is given by:
latitude = ((j-181)/2) + 0.25
longitude = ((i-361)/2) + 0.25
share|improve this question
add comment

1 Answer 1

up vote 3 down vote accepted

We can avoid math by using appropriate software. Because R is free and handy, here is a solution in that platform. It uses some math purely in a data-processing capacity to split the two parts of the [Z05Code] field into row and column indexes.

# Read the file (details depend on the location and name on your system)
data <- read.csv("f:/temp/hg/hg.txt")

# Extract the `i` and `j` indexes from the [Z05Code] field
# Reference: http://www.amap.no/Resources/HgEmissions/HgInventoryDocs.html#Grids
j <- data$Z05Code %/% 1000 # Prefix: latitude index
i <- data$Z05Code %% 1000  # Three-digit suffix: longitude index

The challenge is to convert what is effectively an (i,j,value) format into an actual array of values. R makes this easy, because it allows indexing into an array using the (i,j) pairs.

# Convert one of the data fields to an array.
# For this example, [HgT_00_a] will be used.  Execute `names(data)` for a list
# .. of all field names.
a <- matrix(NaN, ncol=720, nrow=360)  # Start with numeric NoData values
a[cbind(361-j,i)] <- data$HgT_00_a    # Notice how rows need to be reversed!

The rest is routine. We probably want to create a disk dataset in a standard grid format. But first we have to specify the coordinate system: the math gets done here.

# Convert to raster format in order to exploit the `raster` library capabilities.
library(raster) 
a.raster <- raster(a, crs="+proj=longlat +datum=WGS84")
extent(a.raster) <- c(-180,180, -90,90)

It's always a good idea to visually check one's work. (This check originally showed me I had to reverse the rows in the array.)

# Double-check by making a map.
# (To see detail, we plot the log Hg emissions.)
plot(log(a.raster), main="Log(HtT_00_a)")

Map

Once you're happy, create the output.

# Use ?writeRaster to see more formatting options.
writeRaster(a.raster, "F:/temp/HgT_00_a.asc", format="ascii")

The "ascii" output format is the ESRI ASCII grid format, which is readily imported into ArcGIS.

Don't go through this manually for each of the 48 value fields! R has many iteration constructs that can automate producing disk datasets for every field.

share|improve this answer
1  
Very cool. It didn't occur to me to use R, but then again I only have a basic understanding of the software. I'll give that a try. Cheers. –  dchaboya May 25 '12 at 15:46
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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