Recreating A Lat/Long Grid Using a Specified Formula

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

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)

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

• 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