[R-sig-Geo] Accessing gshhs database of world coastlines

[Roger Bivand]

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