[R-sig-Geo] Strange lines when projecting a world map

[Roger Bivand]

On Fri, 6 Dec 2013, Francesco Carotenuto wrote:

> Roger you are right, my explanation is still a bit confusing. Keeping 
> apart the automatic procedure, that is not so much important now, the 
> real problem, as you said, is only aesthetic. When the dataset is 
> worldwide distributed, I don't need to change the default longitude, but 
> when, for instance, I have to consider a dataset distributed in Africa 
> and the whole Eurasia, I need to focus my study on the particular 
> distribution of the data. The grid itself is used to sample values in 
> each cell from a points distribution and, to keep constant the sampling 
> power, I need an equal area projection that can fit well this so huge 
> geographical extent. But in this procedure I don't need  a projected map 
> of the word. The process of intersecting the grid by using the world map 
> is only an aesthetic need to better interpret (and submit) the results 
> of the interpolation.

If you don't need the whole world, then it may make sense to project just 
the continents you do need, avoiding the lines running across the "back" 
of the map. Using wrld_simpl:

library(maptools)
data(wrld_simpl)
EAAw <- wrld_simpl[wrld_simpl$REGION==2 | wrld_simpl$REGION==142 |
  wrld_simpl$REGION==150,]
library(rgdal)
EAAwmoll <- spTransform(EAAw, CRS("+proj=moll +lon_0=66"))
plot(EAAwmoll)

or:

library(rworldmap)
data(countriesCoarse)
EAA <- countriesCoarse[grep("Asia|Europe|Africa", countriesCoarse$GBD),]
library(rgdal)
EAAmoll <- spTransform(EAA, CRS("+proj=moll +lon_0=66"))
plot(EAAmoll)

Both of these give Eurasia+Africa maps that can be projected without the 
stripe artefacts. I suggest doing this, and staying global on +lon_0=0 for 
the global datasets. The rworldmap package is better supported than the 
wrld_simpl dataset in maptools. All of them have geometry issues probably 
introduced under line generalisation from detailed coastlines. It may be 
possible to use maptools::Rgshhs() as an alternative, but here there are 
no country polygons.

Hope this helps,

Roger

> As you can see in the attached figure, the

PS. Your reply didn't reach the list, possibly because of a large 
attachment. On the list, it is best not to attach anything, but to provide 
a link to where the material can be accessed (3000 subscribers are not 
equally interested in all topics).

> interpolation output is difficult to interpret because you can’t 
> understand where are African data and Eurasian ones (the example A in 
> the figure). In the panel C of the same figure I removed the sea portion 
> by intersecting the grid with the world map and this process produces 
> the problem we are talking about. In the panel B I can simply overlap 
> the grid and the world map but the “strange lines” are still there. I 
> could manage the output in the panel A of the figure externally, by 
> using image or gis softwares, but as I’m creating a script that allows 
> one to perform automatically analyses like this, I’d like to provide the 
> users a "finished result". Honestly, the polygons of the states are not 
> necessary because these are pure "macroecological" analyses and I don't 
> need political boundaries, but I really need the continental ones and to 
> exclude sea/ocean portions correctly.
>
>  
> Dr. Francesco Carotenuto (Ph.D.)
> Department of Earth Sciences
> Federico II University
> Largo S. Marcellino 10,
> 80138 NAPLES (ITALY)
> http://francesco-carotenuto.blogspot.it/
> http://3dpaleontology.blogspot.it/ 
>
>
>
> Il Giovedì 5 Dicembre 2013 23:58, Roger Bivand <Roger.Bivand at nhh.no> ha scritto:
> 
> On Thu, 5 Dec 2013, Francesco Carotenuto wrote:
>
>> Thank you for the answer. The object is the prediction location in an 
>> interpolation process. The function of "overlay" is involved in an 
>> automatic process and the longitude (and the kind of projection too) can 
>> change according to the dataset that is, in some cases, worldwide 
>> distributed.
>
> This doesn't help much - do you mean that the grid is of prediction 
> locations? Are the polygons used for anything other than visualization? If 
> you change the visualization projection (shifting the meridian (+lon_0) 
> for worldwide distributions) repeatedly, you will have to clip to suit. If 
> visualization, does it make sense to modify the meridian (the reader may 
> see the shifts as confusing)?
>
> You are trying to deal with multiple somethings that you are interpolating 
> (to a grid?) and displaying in possibly different windows possibly 
> depending on the extent of the prediction grid. Your workflow isn't 
> obvious.
>
> You haven't explained why you need the projection so far (other than 
> aesthetics, maybe?). You may be able to interpolate on the sphere, so 
> unless the target grid cells must be of fixed area - and possibly not even 
> then, you don't need to make life complicated. There are problems in 
> creating global visualizations in any case, as the lines/polygons need 
> clipping at +/- 180 degrees from the meridian.
>
> Roger
>
>>
>>
>>  
>> Dr. Francesco Carotenuto (Ph.D.)
>> Department of Earth Sciences
>> Federico II University
>> Largo S. Marcellino 10,
>> 80138 NAPLES (ITALY)
>> http://francesco-carotenuto.blogspot.it/
>> http://3dpaleontology.blogspot.it/ 
>>
>>
>>
>> Il Giovedì 5 Dicembre 2013 21:23, Roger Bivand <Roger.Bivand at nhh.no> ha scritto:
>> 
>> On Thu, 5 Dec 2013, Francesco Carotenuto wrote:
>>
>>> Dear list, as in
>>> the title, I found that when projecting a world map by mollweide projection and
>>> setting a?? longitude value, some strange
>>> parallels are drawn. The real problem is that these lines ???interact??? with a
>>> raster when performing the ???over??? function.
>>> How can I manage
>>> it?
>>
>> The underlying question is what you want to do with the objects after 
>> projection. The problem in your example is that some polygons/lines cross 
>> -180+66 degrees, so get reflected back across. A possible solution is to 
>> clip the polygons before projection at -180+66 degrees +/- a small number, 
>> but before looking at this in more detail, I need to know what the 
>> projected object is intended for. When +lon_0=0, there is no such effect, 
>> as this data set is clipped at -180 already.
>>
>> Roger
>>
>>
>>> Here I provide a
>>> simulated example to let you understand what I mean
>>> Thanks.??
>>>
>>>
>>> #Example of a
>>> world map mollweide projection with a specific longitude value
>>> library(sp)
>>> library(maptools)
>>> library(rgdal)
>>> data(wrld_simpl)
>>> w_pol<-wrld_simpl
>>> w_pol<-spTransform(w_pol,
>>> CRS("+proj=moll +lon_0=66"))
>>> plot(w_pol)
>>> ??
>>> # Intersecting a
>>> grid with the projected world map
>>> gr<-expand.grid(x
>>> = seq(-180, 180, length.out = 100), y = seq(-90, 90, length.out = 100))
>>> coordinates(gr)<-~x+y
>>> proj4string(gr)<-CRS("+proj=longlat
>>> +ellps=WGS84 +datum=WGS84 +no_defs")
>>> gr<-spTransform(gr,
>>> CRS("+proj=moll +lon_0=66"))
>>> plot(gr[w_pol,])
>>> ??
>>>
>>> ??
>>> Dr. Francesco Carotenuto (Ph.D.)
>>> Department of Earth Sciences
>>> Federico II University
>>> Largo S. Marcellino 10,
>>> 80138 NAPLES (ITALY)
>>> http://francesco-carotenuto.blogspot.it/
>>> http://3dpaleontology.blogspot.it/??
>>>     [[alternative HTML version deleted]]
>>>
>>>
>>
>> -- 
>> Roger Bivand
>> Department of Economics, Norwegian School of Economics,
>> 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
>
>>
>
> -- 
> Roger Bivand
> Department of Economics, Norwegian School of Economics,
> 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
>

-- 
Roger Bivand
Department of Economics, Norwegian School of Economics,
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