[R-sig-Geo] Datum/CRS for the worldHires map data

[Karl Ove Hufthammer]

On Mon, 11 Jan 2010 14:45:42 +0100 Karl Ove Hufthammer <karl at huftis.org> 
wrote:
>> The reason is obviously that the worldHires map is not in the WGS 84 
>> datum (which is reasonable, considering the age of the data). But 
>> which datum is it in?

In case anybody is interested in this information: I have contacted Paul 
Wessel, one of the developers of the GSHHS data. It looks like the 
original data is poorly documented, but the WDB data *might* be in
WGS 72. However, there are probably also other systematic effects such 
as bad digitisation, format conversions &c. I guess this is the reason 
behind the inaccuracy I observed (in *practice* there are not very big 
differences between WGS 72 and WGS 84, even though the two models are 
completely different).

However, the story doesn't end here. It looks the 'worldHires' data in 
'mapdata' is actually different from the version at
http://rimmer.ngdc.noaa.gov/ (select 'World Data Bank II from the drop-
down menu). It's not better, just different.

For instance, compare the area at long = c(25,26), lat = c(70.8, 71.05).

Here is some code that does this. It assumes you have already downloaded 
the GSHHS data (see the Rgshhs help page) and the WDB data from
http://rimmer.ngdc.noaa.gov/ to the currenct directory.

As you can see, the WDB data is shifted from the (correct) GSHHS data. 
But also *some* of the WDB polygons differ between the worldHires 
version and NGDC version. For instance, the polygon for the upper island 
is identical in the two data sets, but the lower land area polygon is 
shifted.

library(maps)
library(mapdata)
library(maptools)

xlim=c(25, 26)
ylim=c(70.8, 71.05)

map1 = Rgshhs("gshhs/gshhs_f.b", xlim=xlim, ylim=ylim, no.clip=TRUE)
plot(map1$SP, xlim=xlim, ylim=ylim, col="black")
map("worldHires", xlim=xlim, ylim=ylim, add=TRUE, col="red", lw=3)
map2 = read.table("24513.dat") # Change the filename here
lines(map2,col="green", lw=2)

-- 
Karl Ove Hufthammer