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