[R-sig-Geo] Rgshhs for getting maps around longitude 0

[Karl Ove Hufthammer]

Roger Bivand wrote:

>> Is it necessary to fetch the map twice? If so, which xlim values should
>> be used, and how should the two maps be merged into one SpatialPolygons?
>> Note that if I change ‘xlim’ to c(0,30), Great Britain is missing ... :-(
> 
> As things stand, I suspect that there is no alternative, but will take a
> look (and encourage others to do so too!). Merging with
> unionSpatialPolygons() as usual.

OK. I’ve created a function for doing this automatically. I haven’t tested
it *very* much, but it seems to work fine for both negative coordinates,
positive coordinates and mixed coordinates.

Having to change the IDs (row names) like this just to merge (rbind) two 
polygons seems a bit too complicated, but I guess it’s necessary?

getRgshhsMap = function (xl, yl, fn = system.file("share/gshhs_c.b", package = "maptools"), 
        shift = TRUE, level = 1, no.clip = FALSE) 
{
    # First try fetching the map directly, possibly with negative coordinates.
    # Note that some polygons with longitude < 0 use negative coordinates 
    # (e.g., Great Britain), and some use positve coordinates (e.g., Ireland).
    #    
    # (Must use 'try' here, because for example xl=c(-40,-10)
    # results in an error, while xl=c(-40,-5) does not.)
    map1 = try(Rgshhs(fn, xlim = xl, ylim = yl, shift = shift, 
                    level = level, no.clip = no.clip)$SP)
    
    # Now try fetching the same area using positive coordinates.
    xl.west = (xl + 360)%%360
    if (xl.west[2] < xl.west[1])
        xl.west[2] = 360
    map2 = Rgshhs(fn, xlim = xl.west, ylim = yl, shift = shift, 
            level = level, no.clip = no.clip)$SP
    
    # If there where no polygons with negative coordinates, just
    # use the positive coordinates.
    if (class(map1) == "try-error") 
        map.union = map2 else { # Else merge the two maps into one.
        
        # First store the original polygon IDs in data frames.
        df1 = data.frame(polyID = row.names(map1))
        row.names(df1) = df1$polyID
        map1.spdf = SpatialPolygonsDataFrame(map1, df1)
        
        df2 = data.frame(polyID = row.names(map2))
        row.names(df2) = df2$polyID
        map2.spdf = SpatialPolygonsDataFrame(map2, df2)
        
        # Generate new polygon IDs to avoid duplicate IDs when
        # rbinding the two maps.
        row.names(map1.spdf) = as.character(seq_along(map1 at polygons))
        row.names(map2.spdf) = as.character(length(map1 at polygons) + 
                        seq_along(map2 at polygons))
        map.merged = rbind(map1.spdf, map2.spdf)
        
        # Finally, combine all the polygons, using the
        # original polyon IDs.
        map.union = unionSpatialPolygons(map.merged, map.merged$polyID)
    }
    map.union
}

-- 
Karl Ove Hufthammer