mapPolygon(coastlineWorld) has horiz. lines in some projections #388

Closed
richardsc opened this Issue Jan 27, 2014 · 14 comments

Comments

Projects
None yet
2 participants
@richardsc
Collaborator

richardsc commented Jan 27, 2014

Example:

library(oce)
data(coastlineWorld)
mapPlot(coastlineWorld)  # default mollweide
mapPolygon(coastlineWorld, col='red')

screen shot 2014-01-27 at 3 44 25 pm

@dankelley

This comment has been minimized.

Show comment
Hide comment
@dankelley

dankelley Jan 27, 2014

Owner

Thanks for posting this. It's long been a known problem. It may have come and gone, as various schemes have been tested for fixing the bug. It needs to be an open issue.

The problem arises at the "cut points" as a land mass goes off one edge of the plot and reenters the other.

I've tried some solution methods based on locating line segments that are a lot longer than others, but that is silly (e.g. what if a person wants to draw a line along the equator). A better method would perhaps look for points that are (in some sense) far apart in map-projected space, but (in some sense) near to each other on the sphere.

Exactly what to do with the wild points is a question. Inserting NA will keep lines from crossing the plot but it doesn't help with filled curves. For little islands, filling is not an issue but a filled continent will be a problem, and it should really be solved by adding fake points that lie along the edge of the graph. But here's the thing. For some projections (like those illustrated above) we could figure out where the edge is, but that is tricky in general.

The problem exists in mapproj, as illustrated below. I wonder if it would make sense to go back to the code upon which mapproj is based (proj4, I think) and see if there is more info that could be used. For example, I am pretty sure that proj4 has inverse functions for all projections, and that could prove handy.

Test code showing mapproj has the problem:

library(oce)
library(mapproj)
xy <- mapproject(coastlineWorld[['longitude']], coastlineWorld[['latitude']], proj="mollweide")
plot(xy$x, xy$y, type='l', asp=1)

mapproj

Owner

dankelley commented Jan 27, 2014

Thanks for posting this. It's long been a known problem. It may have come and gone, as various schemes have been tested for fixing the bug. It needs to be an open issue.

The problem arises at the "cut points" as a land mass goes off one edge of the plot and reenters the other.

I've tried some solution methods based on locating line segments that are a lot longer than others, but that is silly (e.g. what if a person wants to draw a line along the equator). A better method would perhaps look for points that are (in some sense) far apart in map-projected space, but (in some sense) near to each other on the sphere.

Exactly what to do with the wild points is a question. Inserting NA will keep lines from crossing the plot but it doesn't help with filled curves. For little islands, filling is not an issue but a filled continent will be a problem, and it should really be solved by adding fake points that lie along the edge of the graph. But here's the thing. For some projections (like those illustrated above) we could figure out where the edge is, but that is tricky in general.

The problem exists in mapproj, as illustrated below. I wonder if it would make sense to go back to the code upon which mapproj is based (proj4, I think) and see if there is more info that could be used. For example, I am pretty sure that proj4 has inverse functions for all projections, and that could prove handy.

Test code showing mapproj has the problem:

library(oce)
library(mapproj)
xy <- mapproject(coastlineWorld[['longitude']], coastlineWorld[['latitude']], proj="mollweide")
plot(xy$x, xy$y, type='l', asp=1)

mapproj

@ghost ghost assigned dankelley Jan 27, 2014

@dankelley

This comment has been minimized.

Show comment
Hide comment
@dankelley

dankelley Jan 27, 2014

Owner

A quick test indicates that proj4 does a better job with projections than mapproj. Another nice thing is that proj4 offers inverse calculations, which I think could be very handy. (At the moment, I kludge an inverse calculation by doing an optimization problem, which is expensive and can fail.)

Below is test code

library(oce)
library(proj4)
library(mapproj)
data(coastlineWorld)
lon <- coastlineWorld[['longitude']]
lat <- coastlineWorld[['latitude']]

par(mfrow=c(2,1), mar=c(3, 3, 1, 1), mgp=c(2, 0.7, 0))
xy <- project(cbind(lon,lat), "+proj=moll")
plot(xy[,1], xy[,2], type='l', asp=1)

xy <- mapproject(coastlineWorld[['longitude']], coastlineWorld[['latitude']], proj="mollweide")
plot(xy$x, xy$y, type='l', asp=1)

unknown

Owner

dankelley commented Jan 27, 2014

A quick test indicates that proj4 does a better job with projections than mapproj. Another nice thing is that proj4 offers inverse calculations, which I think could be very handy. (At the moment, I kludge an inverse calculation by doing an optimization problem, which is expensive and can fail.)

Below is test code

library(oce)
library(proj4)
library(mapproj)
data(coastlineWorld)
lon <- coastlineWorld[['longitude']]
lat <- coastlineWorld[['latitude']]

par(mfrow=c(2,1), mar=c(3, 3, 1, 1), mgp=c(2, 0.7, 0))
xy <- project(cbind(lon,lat), "+proj=moll")
plot(xy[,1], xy[,2], type='l', asp=1)

xy <- mapproject(coastlineWorld[['longitude']], coastlineWorld[['latitude']], proj="mollweide")
plot(xy$x, xy$y, type='l', asp=1)

unknown

@richardsc

This comment has been minimized.

Show comment
Hide comment
@richardsc

richardsc Jan 30, 2014

Collaborator

A downside to switching to proj4 is that the syntax (for projections, etc) will change. I wonder if you should create a branch for testing using proj4 instead of mapproj so we can thoroughly examine all the potential consequences of switching.

Collaborator

richardsc commented Jan 30, 2014

A downside to switching to proj4 is that the syntax (for projections, etc) will change. I wonder if you should create a branch for testing using proj4 instead of mapproj so we can thoroughly examine all the potential consequences of switching.

@dankelley

This comment has been minimized.

Show comment
Hide comment
@dankelley

dankelley Jan 30, 2014

Owner

Yes, the plan is for a separate branch. (Indeed, with the name you suggest).

Also, my plan was to convert from mapproj syntax to proj4 syntax, because the latter is sort of weird looking. There are only a few dozen projections in consideration, so the names will be easy. It's just a matter of figuring out the "params" and "orientation" parameters.

We are on the exact same wavelength, as so often happens.

I may play with this on Saturday morning.

Owner

dankelley commented Jan 30, 2014

Yes, the plan is for a separate branch. (Indeed, with the name you suggest).

Also, my plan was to convert from mapproj syntax to proj4 syntax, because the latter is sort of weird looking. There are only a few dozen projections in consideration, so the names will be easy. It's just a matter of figuring out the "params" and "orientation" parameters.

We are on the exact same wavelength, as so often happens.

I may play with this on Saturday morning.

@dankelley

This comment has been minimized.

Show comment
Hide comment
@dankelley

dankelley Feb 15, 2014

Owner

I did some tests, and proj4 doesn't look any better than mapproject. This can be seen in the following sample code. However, I am not done with proj4 yet, since I like the fact that it can do inverses. So I could imagine a scheme where bad points are found (as below) and then maybe a scheme involving the inverse could be used to alter the data. (Right now, Oce just inserts NA when there are problems, which is of no help for filled land masses.)

if (!interactive()) png("01.png",
                        width=7, height=4.5, unit="in", res=150, pointsize=8)
## Row 1: mollweide (default)
## Row 2: mollweide (shifted to dateline)
showbad <- function(x, y, col='red', pch=16)
{
    span <- 0.5 * diff(range(x, na.rm=TRUE))
    dx <- c(0, diff(x))
    bad <- abs(dx) > span / 2
    bad[is.na(bad)] <- FALSE
    points(x[bad], y[bad], col=col, pch=pch)
    cat("span:", span, "\n")
}

library(oce)
library(proj4)
library(mapproj)
data(coastlineWorld)
lon <- coastlineWorld[['longitude']]
lat <- coastlineWorld[['latitude']]

par(mfrow=c(2,2), mar=c(3, 3, 1, 1), mgp=c(2, 0.7, 0), cex=2/3)

xymapproj <- mapproject(lon, lat, proj="mollweide")
plot(xymapproj$x, xymapproj$y, type='l', asp=1)
showbad(xymapproj$x, xymapproj$y)

xyproj4 <- project(cbind(lon,lat), "+proj=moll")
plot(xyproj4[,1], xyproj4[,2], type='l', asp=1)
showbad(xyproj4[,1], xyproj4[,2])

## oriented with dateline at centre
xymapproj <- mapproject(lon, lat, proj="mollweide", orientation=c(90,-180,0))
plot(xymapproj$x, xymapproj$y, type='l', asp=1)
showbad(xymapproj$x, xymapproj$y)
xyproj4 <- project(cbind(lon, lat),
                   list(proj="moll", lon_0=-180))
plot(xyproj4[,1], xyproj4[,2], type='l', asp=1)
showbad(xyproj4[,1], xyproj4[,2])
if (!interactive()) dev.off()

01

Owner

dankelley commented Feb 15, 2014

I did some tests, and proj4 doesn't look any better than mapproject. This can be seen in the following sample code. However, I am not done with proj4 yet, since I like the fact that it can do inverses. So I could imagine a scheme where bad points are found (as below) and then maybe a scheme involving the inverse could be used to alter the data. (Right now, Oce just inserts NA when there are problems, which is of no help for filled land masses.)

if (!interactive()) png("01.png",
                        width=7, height=4.5, unit="in", res=150, pointsize=8)
## Row 1: mollweide (default)
## Row 2: mollweide (shifted to dateline)
showbad <- function(x, y, col='red', pch=16)
{
    span <- 0.5 * diff(range(x, na.rm=TRUE))
    dx <- c(0, diff(x))
    bad <- abs(dx) > span / 2
    bad[is.na(bad)] <- FALSE
    points(x[bad], y[bad], col=col, pch=pch)
    cat("span:", span, "\n")
}

library(oce)
library(proj4)
library(mapproj)
data(coastlineWorld)
lon <- coastlineWorld[['longitude']]
lat <- coastlineWorld[['latitude']]

par(mfrow=c(2,2), mar=c(3, 3, 1, 1), mgp=c(2, 0.7, 0), cex=2/3)

xymapproj <- mapproject(lon, lat, proj="mollweide")
plot(xymapproj$x, xymapproj$y, type='l', asp=1)
showbad(xymapproj$x, xymapproj$y)

xyproj4 <- project(cbind(lon,lat), "+proj=moll")
plot(xyproj4[,1], xyproj4[,2], type='l', asp=1)
showbad(xyproj4[,1], xyproj4[,2])

## oriented with dateline at centre
xymapproj <- mapproject(lon, lat, proj="mollweide", orientation=c(90,-180,0))
plot(xymapproj$x, xymapproj$y, type='l', asp=1)
showbad(xymapproj$x, xymapproj$y)
xyproj4 <- project(cbind(lon, lat),
                   list(proj="moll", lon_0=-180))
plot(xyproj4[,1], xyproj4[,2], type='l', asp=1)
showbad(xyproj4[,1], xyproj4[,2])
if (!interactive()) dev.off()

01

@dankelley dankelley added the mapping label May 13, 2014

@dankelley dankelley changed the title from mapPolygon(coastlineWorld) has horizontal lines in some projections to mapPolygon(coastlineWorld) has horiz. lines in some projections Jun 5, 2014

@dankelley

This comment has been minimized.

Show comment
Hide comment
@dankelley

dankelley Jul 17, 2014

Owner

Those interested in this might want to pull the oce-issues repo and do

    cd 388
    make
    make view

and check out the “C” plot.

What I’ve done is to write a C function to shove a NA between any two points that are far apart. This seems useful for coastlines drawn with lines. I assume it will be messed up in some polygon-fill cases, but maybe not so bad because it’s just some bits that are on the edge of the graph anyhow.

The graph of that C part looks like the following. Note that the top-right panel lacks those nasty horizontal lines.
screen shot 2014-07-17 at 11 02 57 am

Owner

dankelley commented Jul 17, 2014

Those interested in this might want to pull the oce-issues repo and do

    cd 388
    make
    make view

and check out the “C” plot.

What I’ve done is to write a C function to shove a NA between any two points that are far apart. This seems useful for coastlines drawn with lines. I assume it will be messed up in some polygon-fill cases, but maybe not so bad because it’s just some bits that are on the edge of the graph anyhow.

The graph of that C part looks like the following. Note that the top-right panel lacks those nasty horizontal lines.
screen shot 2014-07-17 at 11 02 57 am

@dankelley

This comment has been minimized.

Show comment
Hide comment
@dankelley

dankelley Jul 17, 2014

Owner

Well, there is one wonky line left, this time at a 45 deg angle (to NW) from a spot in Africa. But this probably just means some tweaking is required. For example, perhaps should only look for weird line segments that have one end at high abs(x) value.

Naturally, there is basically nothing to be done about Antarctica. Having looked at a few datasets, I know that some odd things are generally done there. Sigh.

Owner

dankelley commented Jul 17, 2014

Well, there is one wonky line left, this time at a 45 deg angle (to NW) from a spot in Africa. But this probably just means some tweaking is required. For example, perhaps should only look for weird line segments that have one end at high abs(x) value.

Naturally, there is basically nothing to be done about Antarctica. Having looked at a few datasets, I know that some odd things are generally done there. Sigh.

@dankelley

This comment has been minimized.

Show comment
Hide comment
@dankelley

dankelley Jan 28, 2015

Owner

The problem is not unique to oce. I just tried mapdata and got similar.

screen shot 2015-01-28 at 11 47 53 am

Owner

dankelley commented Jan 28, 2015

The problem is not unique to oce. I just tried mapdata and got similar.

screen shot 2015-01-28 at 11 47 53 am

@dankelley

This comment has been minimized.

Show comment
Hide comment
@dankelley

dankelley May 26, 2015

Owner

Possibly sp:: nowrapSpatialLines() will help, although the sp author suggests that existing solutions are not ideal.

Owner

dankelley commented May 26, 2015

Possibly sp:: nowrapSpatialLines() will help, although the sp author suggests that existing solutions are not ideal.

@dankelley

This comment has been minimized.

Show comment
Hide comment
@dankelley

dankelley Jun 1, 2015

Owner

I wrote a new function that gets fairly close to a solution. It has some "opposite-edge shadowing" (see e.g. the edges of Antarctica, which have lines going north to the northernmost point of the continent). I may do some work to remove that later, but that can be done without changing the UI.

The test code only works on the branch named project_with_rgdal, and I do not plan on merging that into develop for several days at the earliest. Below is the test code named 640/09.R.

Notice that the new function coastlineCut() may be changed without notice.

After the code is a snap of the first page. If you build the whole PDF, you can page through and get a mini-movie of a spinning globe.

library(oce)
data(coastlineWorld)
if (!interactive()) pdf("09.pdf", width=7, height=7, pointsize=8)
par(mfrow=c(2,1), mar=c(2, 2, 1, 1), mgp=c(2, 0.7, 0))
xlim <- ylim <- NULL # yields identical map scales on successive pages
for (lon_0 in seq(-180, 180, 10)) {
    proj <- "robin"
    mod <- coastlineCut(coastlineWorld, lon_0=lon_0)
    lon <- mod[["longitude"]]
    lat <- mod[["latitude"]]
    plot(lon, lat, xlim=c(-180,180), ylim=c(-90,90), type='l')
    lines(c(-180, 180, 180, -180, -180), c(-90, -90, 90, 90, -90), col='gray')
    polygon(lon, lat, col='gray')
    proj <- sprintf("+proj=%s +lon_0=%.0f", proj, lon_0)
    if (is.null(xlim)) {
        mapPlot(mod, fill='gray', proj=proj)
        xlim <- par('usr')[1:2]
        ylim <- par('usr')[3:4]
    } else {
        mapPlot(mod, fill='gray', proj=proj, xlim=xlim, ylim=ylim, xaxs="i", yaxs="i")
    }
    mtext("BUG: opposite-edge 'shadow'", line=0, col='magenta', font=2, adj=0)
    mtext(proj, side=3, line=0, adj=1)
}
if (!interactive()) dev.off()

screen shot 2015-06-01 at 11 24 56 am

Owner

dankelley commented Jun 1, 2015

I wrote a new function that gets fairly close to a solution. It has some "opposite-edge shadowing" (see e.g. the edges of Antarctica, which have lines going north to the northernmost point of the continent). I may do some work to remove that later, but that can be done without changing the UI.

The test code only works on the branch named project_with_rgdal, and I do not plan on merging that into develop for several days at the earliest. Below is the test code named 640/09.R.

Notice that the new function coastlineCut() may be changed without notice.

After the code is a snap of the first page. If you build the whole PDF, you can page through and get a mini-movie of a spinning globe.

library(oce)
data(coastlineWorld)
if (!interactive()) pdf("09.pdf", width=7, height=7, pointsize=8)
par(mfrow=c(2,1), mar=c(2, 2, 1, 1), mgp=c(2, 0.7, 0))
xlim <- ylim <- NULL # yields identical map scales on successive pages
for (lon_0 in seq(-180, 180, 10)) {
    proj <- "robin"
    mod <- coastlineCut(coastlineWorld, lon_0=lon_0)
    lon <- mod[["longitude"]]
    lat <- mod[["latitude"]]
    plot(lon, lat, xlim=c(-180,180), ylim=c(-90,90), type='l')
    lines(c(-180, 180, 180, -180, -180), c(-90, -90, 90, 90, -90), col='gray')
    polygon(lon, lat, col='gray')
    proj <- sprintf("+proj=%s +lon_0=%.0f", proj, lon_0)
    if (is.null(xlim)) {
        mapPlot(mod, fill='gray', proj=proj)
        xlim <- par('usr')[1:2]
        ylim <- par('usr')[3:4]
    } else {
        mapPlot(mod, fill='gray', proj=proj, xlim=xlim, ylim=ylim, xaxs="i", yaxs="i")
    }
    mtext("BUG: opposite-edge 'shadow'", line=0, col='magenta', font=2, adj=0)
    mtext(proj, side=3, line=0, adj=1)
}
if (!interactive()) dev.off()

screen shot 2015-06-01 at 11 24 56 am

@richardsc

This comment has been minimized.

Show comment
Hide comment
@richardsc

richardsc Jun 1, 2015

Collaborator

Great, I'll try testing this later when I get a chance. Though, from the pick above it looks like we've replaced UHL with some UVL! Still, it's an improvement! :)

Collaborator

richardsc commented Jun 1, 2015

Great, I'll try testing this later when I get a chance. Though, from the pick above it looks like we've replaced UHL with some UVL! Still, it's an improvement! :)

@dankelley

This comment has been minimized.

Show comment
Hide comment
@dankelley

dankelley Jun 1, 2015

Owner

Actually the UVL would go away if mapPlot() traced the "edge" of the earth, which I actually think would make sense because longitude lines don't always get drawn there, depending on the rotation. But tracing is a little bit tricky. (It's easy to trace for some projections, harder in others, so the best is to write general code that essentially searches for spots that can do UHL.)

PS. I feel that I'll need to merge into develop before working on any other map-related issue. I've made lots of changes, including some yesterday that get around an uninvertable zone in wintri projections (something I had fixed within the modified PROJ.4 sources previously).

Owner

dankelley commented Jun 1, 2015

Actually the UVL would go away if mapPlot() traced the "edge" of the earth, which I actually think would make sense because longitude lines don't always get drawn there, depending on the rotation. But tracing is a little bit tricky. (It's easy to trace for some projections, harder in others, so the best is to write general code that essentially searches for spots that can do UHL.)

PS. I feel that I'll need to merge into develop before working on any other map-related issue. I've made lots of changes, including some yesterday that get around an uninvertable zone in wintri projections (something I had fixed within the modified PROJ.4 sources previously).

@richardsc

This comment has been minimized.

Show comment
Hide comment
@richardsc

richardsc Nov 3, 2015

Collaborator

This issue is not gone, but I have found that coastlineCut() works for quite a few cases. Relatedly, while searching for something else I found this blog posting:

http://cameron.bracken.bz/finally-an-easy-way-to-fix-the-horizontal-lines-in-ggplot2-maps

and thought there might be some ideas in there. In particular, it seems that they use the PBSmapping package for do the polygon clipping.

Collaborator

richardsc commented Nov 3, 2015

This issue is not gone, but I have found that coastlineCut() works for quite a few cases. Relatedly, while searching for something else I found this blog posting:

http://cameron.bracken.bz/finally-an-easy-way-to-fix-the-horizontal-lines-in-ggplot2-maps

and thought there might be some ideas in there. In particular, it seems that they use the PBSmapping package for do the polygon clipping.

@dankelley

This comment has been minimized.

Show comment
Hide comment
@dankelley

dankelley Jun 28, 2017

Owner

Closing now, since this has been fixed for quite some time. Certainly, the initial test is fine now. It must be noted that coastlineCut() is the best way to shift the central longitude; using +lon_0=... in the projection argument will cause problems. Also, note that the issue with Antarctica still exists (issue #616) but the present issue (as posed) can be closed. DK+CR decision (finally).

Owner

dankelley commented Jun 28, 2017

Closing now, since this has been fixed for quite some time. Certainly, the initial test is fine now. It must be noted that coastlineCut() is the best way to shift the central longitude; using +lon_0=... in the projection argument will cause problems. Also, note that the issue with Antarctica still exists (issue #616) but the present issue (as posed) can be closed. DK+CR decision (finally).

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment