Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Population grid problem - ring is not simple or complete - change of coordinate system and horizontal line due to split polygons solution #25

Open
YanisChrys opened this issue Apr 23, 2018 · 11 comments

Comments

@YanisChrys
Copy link

Hello,

I previously had a problem where I was trying to use a polygon that incorporated samples from the Chukchi peninsula where their coordinates are negative with others from Siberia having positive coordinates. (For example 179.0 and then -179.0 ) which created problems and I ended up not using that area, but now I must use it.

One way I have resolved this in the past is to add 360 to those negative coordinates(for the environment and the individual coordinates) and have 181 for the second one for example. But in EEMS this gives an error saying the “ring is not simple or complete”.

Do you have a way of including such individuals into EEMS from across Siberia? How can I include those individuals when the coordinates make a quick change from 179 to -179 degrees? Is something like that even possible in EEMS?

@dipetkov
Copy link
Owner

Can you share a (simplified) example? [The coordinates & the habitat would be enough, as a dissimilarity matrix can easily be simulated.]

@YanisChrys
Copy link
Author

Yes,
So, as you can see in the example below there is this issue of the coordinates:

161.545571 71.446935
-168.395835 71.423097
-168.395835 56.865992
158.029946 60.167932
161.545571 71.446935

screenshot-2018-4-23 tool for google maps v3 version 3

@dipetkov
Copy link
Owner

Here is my first idea: map the longitude measurements from [-180,+180] degrees to [0, 360]. This is easy: add 360 to the negative longitude values:

long360 <- long180 + (long180 < 0) *  360

This should deal with the "graph is not connected" error.

Then for plotting, if you want to add a world map, you will have to specify in the PROJ.4 string that the longitude degrees range from 0 to 360. I have to look into how to do this (or even if it is possible for that matter.)

(An aside) Of course this doesn't really address the underlying issue that runeems treats longitudes and latitudes as x,y coordinates in 2D spaces and doesn't care about the curvature of the earth one bit.

For habitats that cover only a relatively small area that doesn't matter too much but it does have an effect on large habitats/habitats close to either pole where the distortions are largest.

In theory, I could fix this as follows: take latitudes/longitudes, project them on 2D plane using a projection like Mercator, triangulate in 2D and then map back onto earth. Seems straightforward as long as there are formulas for reversing the projection.

@YanisChrys
Copy link
Author

Hello again,

I have done this for the individual degrees and for the outer.txt file and I use the following command to plot:

eems.plots(mcmcpath,
plotpath = paste("tester",myrun,sep="_"),
longlat = TRUE,
add.outline = TRUE, col.outline = "purple", lwd.outline = 2,
lwd.map = 1.3, remove.singletons = F,
add.map = TRUE,xpd = T,
projection.in = "+proj=longlat +datum=WGS84 +lon_wrap=180",
projection.out = "+proj=merc +lon_0=90 +k=1 +x_0=10018754.171395 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs", #as suggested here for mapping Russia: link
add.title = TRUE, add.colbar = TRUE,
m.plot.xy = { text(coords_merc, pch = labels, font = 2) },
q.plot.xy = { text(coords_merc, pch = labels, font = 2) })

For the projection.in, from what I understand, +lon_wrap is supposed to indicate that the coordinates are 0,360 but removing it doesn't change the output from what I can tell. The solution is given by projection.out.
However, this produces several horizontal lines that span the whole map. From what I've understood this is due to proj4 trying to connect points on opposite sides of the map, but please correct me if I'm wrong.
I have found people trying to correct for this by removing those points but I think the issue here might be the mercator transformed demes not the 5 input individuals I use.

The output looks like this:

image

I have tried using other projections like Lambert Equal distance or the stereographic projection which reduced the number of lines running through the graph but in those cases the EEMS weren't plotted correctly, so I assume they are calculated to match with the mercator projection.

Do you have a solution to this? Is there a way I can get rid of the horizontal lines?

@YanisChrys
Copy link
Author

It appears that the following code was enough to fix this problem:

projection.out = "+proj=merc +lon_0=90 +k_0=2 +x_0=10018754.171395 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs +over"

@YanisChrys YanisChrys changed the title Population grid problem - ring is not simple or complete Population grid problem - ring is not simple or complete - change of coordinate system and horizontal line due to split polygons SOLUTION May 14, 2018
@YanisChrys YanisChrys changed the title Population grid problem - ring is not simple or complete - change of coordinate system and horizontal line due to split polygons SOLUTION Population grid problem - ring is not simple or complete - change of coordinate system and horizontal line due to split polygons solution May 14, 2018
@YanisChrys
Copy link
Author

YanisChrys commented May 14, 2018

It turns out this line of code does avoid the wrapping issues caused by points crossing the dateline however, it doesn't display the part of the map after lon +180.
Many people have tried to solve this, as I understand (1, Geoserver, Mapserver, 4, gdalwarp ),
but, no one can manage to do that with a PROJ4 line of code, so I think it's not possible right now to obtain the above image without the horizontal lines.

Although, I think it's possible to transform from WGS84 to EPSG4326 in a way similar to what they do here unless something more complicated is needed for EEMS.

@YanisChrys YanisChrys reopened this May 14, 2018
@A33miller
Copy link

@Johnchrys I am also having this problem, however I am trying to create a .outer file that widely surrounds New Zealand. Were you ever able to figure it out?

@YanisChrys
Copy link
Author

@Johnchrys I am also having this problem, however I am trying to create a .outer file that widely surrounds New Zealand. Were you ever able to figure it out?

Hey @A33miller, I solved it by adding 180 to all longitude values so they range from 0 to 360. I also had other issues so just tell me if the problem persists and maybe the solutions I found will help you.

@A33miller
Copy link

@Johnchrys I am also having this problem, however I am trying to create a .outer file that widely surrounds New Zealand. Were you ever able to figure it out?

Hey @A33miller, I solved it by adding 180 to all longitude values so they range from 0 to 360. I also had other issues so just tell me if the problem persists and maybe the solutions I found will help you.

Thanks @Johnchrys ! I am also trying to create a habitat/.outer file with holes (i.e., to include only the water surrounding New Zealand) any chance you did this as well? The problem I have is that the birdtheme (kml creator webpage) only creates polygons with holes by listing the points in a counterclockwise and then clockwise fashion, and eems only wants these to be listed in one direction. If you have experience with this I would appreciate your recommendations, if not, dipetkov very nicely said he would offer me some help.

@dipetkov
Copy link
Owner

@Johnchrys I am also having this problem, however I am trying to create a .outer file that widely surrounds New Zealand. Were you ever able to figure it out?

Hey @A33miller, I solved it by adding 180 to all longitude values so they range from 0 to 360. I also had other issues so just tell me if the problem persists and maybe the solutions I found will help you.

Thanks @Johnchrys ! I am also trying to create a habitat/.outer file with holes (i.e., to include only the water surrounding New Zealand) any chance you did this as well? The problem I have is that the birdtheme (kml creator webpage) only creates polygons with holes by listing the points in a counterclockwise and then clockwise fashion, and eems only wants these to be listed in one direction. If you have experience with this I would appreciate your recommendations, if not, dipetkov very nicely said he would offer me some help.

No, for a habitat with holes, the problem is not listing the points counterclockwise in the outer file. The issue is that outer file is expected to contain a "ring" geometry and, in the implementation, the habitat is represented by a ring data structure. That is, a simple closed polygon, without holes.

But you can use a pre-constructed graph (with the --gridpath argument) which doesn't contain demes inside holes. [The habitat itself will still be a ring though...]

@dipetkov
Copy link
Owner

It turns out this line of code does avoid the wrapping issues caused by points crossing the dateline however, it doesn't display the part of the map after lon +180.
Many people have tried to solve this, as I understand (1, Geoserver, Mapserver, 4, gdalwarp ),
but, no one can manage to do that with a PROJ4 line of code, so I think it's not possible right now to obtain the above image without the horizontal lines.

Although, I think it's possible to transform from WGS84 to EPSG4326 in a way similar to what they do here unless something more complicated is needed for EEMS.

It seems that something more complicated might be needed. The migration surface is plotted by the raster package. If we want to use a projection, then the surface is transformed with the raster::projectRaster function. Other elements like the grid and the map are transformed with sp::spTransform.

Take this New Zealand example: the grid and the map are both projected correctly but the migration surface is cut off just like you describe even though the same proj4 string is used to specify the projection. So maybe raster cannot handle this (general) case.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

3 participants