Accessing gshhs database of world coastlines
On Sun, 16 Oct 2005, Denis Chabot wrote:
Hi Roger, Le 05-10-16 ? 14:09, Roger Bivand a ?crit :
On Sun, 16 Oct 2005, Denis Chabot wrote:
...
I thought I did not need to because a shapefile version of the
database exists, and I have managed (with help from Roger in
particular) to import such files in R before. But this one is huge
(154 MB for the shp file alone). The function read.shape refuses to
open it, telling me right way
Erreur dans read.shape("gshhs_land.shp") : unable to open file
This means read.shape() cannot see the file.
Maybe my copy is correupted? I tried again giving the full pathname to the file to insure I had not screwed up with defining my working directory, and got yet the same (or close to):
> test <- read.shape("Macintosh
HD:Users:dchabot:Documents:Cartographie:Shapefiles:gshhs:gshhs_land.shp"
)
Erreur dans read.shape("Macintosh
HD:Users:dchabot:Documents:Cartographie:Shapefiles:gshhs:gshhs_land.shp"
) :
unable to open file
(I had used a utility to give me the full pathname of the shp file).
I have no access to OSX, I'm afraid. If you do setwd() to the directory
where gshhs_land.shp is placed, and then list.files(pattern="shp"), is the
file visible?`What does getwd() say when you are there? The error message
is from:
hSHP = SHPOpen(CHAR(STRING_ELT(shpnm,0)), "rb" );
if( hSHP == NULL )
error("unable to open file");
in src/Rshapeget.c, where shpnm is the file name given to read.shape(). So
either the compilation failed, or the file isn't there. Try
getinfo.shape("gshhs_land.shp")
the error message there should say "Unable to open: gshhs_land.shp".
Alternative explanation: I just upgraded to R2.2, and wiped out my 2.1 install. Because R2.2 is not officially on line at CRAN yet, I could not download mac binaries for maptools and sp and compiled them myself. I did get a warning about some name being used twice, and although the package seems to work (it certainly imports a smaller shapefile) maybe it is not 100% functional.
I managed to read gshhs_land.shp into R with current maptools:
try1 <- read.shape("gshhs_land.shp")
Shapefile type: PolyLine, (3), # of Shapes: 191811
object.size(try1)
[1] 395387692 but it is a shapefile of lines, not polygons - the base GSHHS is now polygonal.
I don't understand this part. the "gshhs_land" is not "base GSHHS"? I thought I downloaded the most recent version, at the highest resolution. Where did you find the version with polygons?
It is not a shapefile, but in GSHHS format. It is formatted as polygons. The problem with the shapefile as lines is that it cannot be used for filling land or sea, just for drawing coastlines.
I have 1GB under linux, so I think memory management is effectively the same as yours. I've put a PDF (33MB) of your area on: http://spatial.nhh.no/misc/nw_atlantic.pdf please save link as, and open locally, otherwise you'll be sitting looking at a blank screen for a long time. I did this so that you can use the Acroread zoom to find out if your islands are there.
Yes the islands are there! This coastline is an improvement on what I have!
OK, so I'll follow this route as well.
read.shape() was not too slow, it seemed to go quite well. Is there any use in writing an R function on top of the gshhs program for extracting polygons from the gshhs databases directly into sp class objects - it doesn't look too difficult, and would give direct access without having to use the shapefile?
I would love to have it. At a minimum, I'd need clear(er) instructions as to how to compile and use the programs provided on GSHHS's site. But until now I valued being able to work with my dataset in R as I am used to, and produce maps for the right subset of data I'm interested in. So I don't plan to use GMT any time soon and would much prefer R functions to access the very nice gshhs databases and take what I need out of them.
I think it can be done.
As an aside, to you know of a way to do map borders in R that look like those in GMT and many paper maps (thick borders, alternating from white to black from one degree of lat or long to the next)? This is a cosmetic issue only but it would be a nice touch to be able to do this in R.
For unprojected it could be done, for projected like GMT is more trouble. If anyone would like to contribute GMTaxis() to work like GMT, please say, contributions welcome!
Anyway, I'll wait a few days for the full R2.2 Mac to be at CRAN and get fresh binaries of maptools and sp to see if I can open that shapefile, since you were able to. It irritates me to have failed.
Roger
With many thanks, Denis
Roger Bivand Economic Geography Section, Department of Economics, Norwegian School of Economics and Business Administration, Helleveien 30, N-5045 Bergen, Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43 e-mail: Roger.Bivand at nhh.no