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(m...@polygons))
row.names(map2.spdf) = as.character(length(m...@polygons) +
seq_along(m...@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
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
