[R-sig-Geo] Good projection for N/S America?

[Roger Bivand]

On Tue, 10 Apr 2007, Roger Bivand wrote:

> On Mon, 9 Apr 2007 White.Denis at epamail.epa.gov wrote:
> 
> > Roger Bivand <Roger.Bivand at nhh.no> wrote on 2007-04-09 13:52:45:
> > 
> > > On Mon, 9 Apr 2007 White.Denis at epamail.epa.gov wrote:
> > >
> > > > Roger's solution makes sense to me.
> > > >
> > > > The sinusoidal does have the appearance of pinching poleward, a
> > > > consequence of allocating equal area by spacing equally in both x
> > and y.
> > >
> > > Both Canada and Norway - especially including Spitzbergen - suffer
> > from
> > > this, but good compromises are hard to find if one needs Patagonia
> > too.
> > 
> > Goode's Homolosine is one.  USGS distributes (or did) some of their
> > global land cover data in that projection, see
> > http://landcover.usgs.gov/globallandcover.php.  Not only is it
> > interrupted, but it combines the sinusoidal from -40 to +40 and the
> > Mollweide poleward of +/- 40.  (I don't see it on the geotiff list.)
> 
> projInfo() in rgdal gives a list of supported projections, here:
> 
> +proj=goode +lon_0=-80
> 
> and a PNG is attached. However, in rgdal on one platform I'm seeing an
> EDOM (trig/math function out of domain) that I don't see running proj
> through system(). The grid is OK, so there may be a weakness in the
> project/transform interface in this case - I'm looking into it. Like the
> PROJ.4 sinusoidal, it disregards the ellipsoid anyway, just spherical.

Following this problem up, I have improved error handling in project() and
spTransform() in rgdal, but there is a problem with "+proj=goode" on some
platforms (not the current Windows binary), which seems to get muddled
when switching across +/- 40d44m (in proj source file pj_goode.c) or
40d44m12s (Snyder p. 249). Running "+proj=moll" or "+proj=sinu" on the
same data is not a problem on the affected platform(s). My suspicions are 
that there are differences in compiler headers between the proj library 
and the rgdal/R binaries, which may be platform dependent - my affected 
platforms are RHEL 4.

Roger

> 
> I haven't checked whether the same offset approach can be used to split
> the Northern and Southern hemispheres, but I don't see why not, given 
> that this is mostly sinusoidal.
> 
> Roger
> 
> > 
> > > I can see that getting rgdal binaries for popular platforms is a real
> > > issue, any suggestions?
> > >
> > > Roger
> > >
> > > > Also the Lambert cylindrical sent in before should have had standard
> > > > parallels set to +/- 30, i.e.,
> > > >
> > > > (See attached file: whemi.projs.png)
> > > >
> > > > Tim Keitt <tkeitt at gmail.com> wrote on 2007-04-09 13:38:33:
> > > >
> > > > > Canada looks pinched in this projection. S. Am is perfect.
> > > > >
> > > > > THK
> > > > >
> > > > > On 4/9/07, Roger Bivand <Roger.Bivand at nhh.no> wrote:
> > > > > > On Mon, 9 Apr 2007 White.Denis at epamail.epa.gov wrote:
> > > > > >
> > > > > > > For preview graphics and for large areas such as continents,
> > large
> > > > > > > countries, hemispheres, or the whole earth, spherical
> > projections
> > > > are
> > > > > > > often adequate.  I can provide some of the ones I have used.
> > For
> > > > > > > detailed work at sites and small areas, ellipsoidal
> > projections
> > > > such as
> > > > > > > UTM are usually used, and then the coding gets more
> > complicated
> > > > with
> > > > > > > choices of datums and so forth.
> > > > > > >
> > > > > >
> > > > > > The attached script shows how to do the interrupted sinusoidal
> > > > projection
> > > > > > using spTransform in rgdal, for the whemi.lin data posted with
> > the
> > > > > > free-standing functions by Denis White a couple of days ago.
> > Once
> > > > the
> > > > > > lines are converted into SpatialLines objects, the rest is
> > robust
> > > > and
> > > > > > simple, as is the use of gridlines() in sp. The one catch is
> > > > calculating
> > > > > > the offset, here in an x_0= offset along the Equator in metres
> > > > between the
> > > > > > two central longitude values. The output is attached as a PNG
> > image.
> > > > The
> > > > > > point about the sp objects is that they contain enough metadata
> > > > (here a
> > > > > > PROJ.4 projection description) to let them be moved to other R
> > > > packages or
> > > > > > external software.
> > > > > >
> > > > > > The half-dozen basic projections are easy to specify in PROJ.4,
> > for
> > > > > > example from the geotiff list:
> > > > > >
> > > > > > http://www.remotesensing.org/geotiff/proj_list/
> > > > > >
> > > > > > which is what I used here. The other projections mentioned are:
> > > > > >
> > > > > > Lambert Cylindrical Equal Area "+proj=cea +lon_0=-80"
> > > > > >
> > > > > > Lambert Azimuthal Equal Area "+proj=laea +lat_0=0 +lon_0=-80"
> > > > > >
> > > > > > while the Northern hemisphere sinusoidal is:
> > > > > >
> > > > > > Sinusoidal "+proj=sinu +lon_0=-100"
> > > > > >
> > > > > > So I'd argue that PROJ.4 projection descriptions are not
> > difficult
> > > > to use,
> > > > > > and with sp objects, do stay stuck to the data (has anyone else
> > ever
> > > > > > forgotten what projection was used when revisiting data, not
> > just
> > > > me?).
> > > > > >
> > > > > > Using the maptools map2SpatialLines() interface function, or the
> > > > Rgshhs()
> > > > > > interface to GSHHS shorelines, even getting the lines is quite
> > easy,
> > > > > > qualified by clipping and bounding box issues in extremities for
> > > > > > projection from geographical coordinates.
> > > > > >
> > > > > > Of course, it would help to have MacOS X and selected Linux
> > binaries
> > > > of
> > > > > > rgdal, we're very lucky that Uwe Ligges is so helpful with the
> > > > Windows
> > > > > > binaries.
> > > > > >
> > > > > > Roger
> > > > > >
> > > > > > > r-sig-geo-bounces at stat.math.ethz.ch wrote on 2007-04-08
> > 07:56:03:
> > > > > > >
> > > > > > > > Denis,
> > > > > > > >
> > > > > > > > That's really useful. It occurs to me that we only really
> > need a
> > > > > > > > half-dozen basic projections to cover 90% of user cases.
> > Perhaps
> > > > these
> > > > > > > > could be incorporated into the 'sp' group somewhere and
> > relieve
> > > > the
> > > > > > > > dependence on proj4. (It could be packaged separately for R
> > for
> > > > the
> > > > > > > > other 10% of cases where its needed.)
> > > > > > > >
> > > > > > > > THK
> > > > > > > >
> > > > > > > > On 4/6/07, White.Denis at epamail.epa.gov
> > > > <White.Denis at epamail.epa.gov>
> > > > > > > wrote:
> > > > > > > > > Thanks, Roger.  There was a request to see the R code for
> > > > these
> > > > > > > figures.
> > > > > > > > > Attached is the script for the second PDF file plus the
> > input
> > > > > > > boundary
> > > > > > > > > file I used for the hemisphere.  The three projection
> > > > functions are
> > > > > > > for
> > > > > > > > > simple spherical, rather than ellipsoidal, models of the
> > > > earth.  The
> > > > > > > > > graticule generating function could be more elegant.  I'm
> > not
> > > > yet up
> > > > > > > to
> > > > > > > > > speed with sp and the many new spatial capabilities in R
> > so
> > > > please
> > > > > > > > > excuse the old style "lines()" format encoding and
> > graphics.
> > > > > > > > >
> > > > > > > > > Tim, I don't know whether proj4 could do the interrupted
> > > > sinusoidal.
> > > > > > > > >
> > > > > > > > > (See attached file: whemi.projs.r)(See attached file:
> > > > whemi.lin)
> > > > > > > > >
> > > > > > > > > r-sig-geo-bounces at stat.math.ethz.ch wrote on 2007-04-06
> > > > 04:51:53:
> > > > > > > > >
> > > > > > > > > > Since this topic is of general interest, I've made an
> > > > exception
> > > > > > > and
> > > > > > > > > > allowed (this once!) a posting of more than 200K. In
> > > > general, if
> > > > > > > > > graphics
> > > > > > > > > > are big, please consider either an alternative device
> > (png
> > > > is
> > > > > > > often
> > > > > > > > > OK),
> > > > > > > > > > or posting just a URL to the real file.
> > > > > > > > > >
> > > > > > > > > > With apologies to list members on dial-up connections in
> > the
> > > > > > > field,
> > > > > > > > > >
> > > > > > > > > > Roger
> > > > > > > > > >
> > > > > > > > > >
> > > > > > > > > > On Thu, 5 Apr 2007 White.Denis at epamail.epa.gov wrote:
> > > > > > > > > >
> > > > > > > > > > > Yes, for many uses that is my choice also.  For the
> > > > conterminous
> > > > > > > US
> > > > > > > > > for
> > > > > > > > > > > example, the Lambert azimuthal has lower mean
> > distortion
> > > > than
> > > > > > > the
> > > > > > > > > > > commonly used standard projection, the Albers conical
> > > > equal
> > > > > > > area,
> > > > > > > > > > > although Albers was chosen by USGS as a standard
> > because
> > > > of
> > > > > > > lower
> > > > > > > > > > > extreme distortion than many other possible
> > projections.
> > > > > > > > > > >
> > > > > > > > > > > For our hemispherical application, because we were
> > > > gridding the
> > > > > > > > > data, we
> > > > > > > > > > > wanted parallels of latitude to be parallel in the
> > > > projected
> > > > > > > > > coordinate
> > > > > > > > > > > space, which we wouldn't get with the Lambert
> > azimuthal.
> > > > > > > > > > >
> > > > > > > > > > > (See attached file: whemi.projs.pdf)
> > > > > > > > > > >
> > > > > > > > > > > Tim Keitt <tkeitt at gmail.com> wrote on 2007-04-05
> > 10:56:09:
> > > > > > > > > > >
> > > > > > > > > > > > Thanks. My application is not that demanding.
> > Really, I
> > > > just
> > > > > > > want
> > > > > > > > > it
> > > > > > > > > > > > to look reasonable. My plan is to lay out the
> > postings
> > > > in the
> > > > > > > > > > > > projected coordinates and then back transform into
> > > > geographic
> > > > > > > > > > > > coordinates for analysis. I tried lots of
> > projections
> > > > and
> > > > > > > found
> > > > > > > > > > > > Lamberts Azimuthal Equal Area to be quite good. I
> > like
> > > > the
> > > > > > > look of
> > > > > > > > > the
> > > > > > > > > > > > Azimuthal Equidistant better, but figured equal area
> > was
> > > > a
> > > > > > > good
> > > > > > > > > > > > choice.
> > > > > > > > > > > >
> > > > > > > > > > > > THK
> > > > > > > > > > > >
> > > > > > > > > > > > On 4/4/07, White.Denis at epamail.epa.gov
> > > > > > > > > <White.Denis at epamail.epa.gov>
> > > > > > > > > > > wrote:
> > > > > > > > > > > > > Tim,
> > > > > > > > > > > > >
> > > > > > > > > > > > > It depends on which kind of distortion is of most
> > > > concern.
> > > > > > > For
> > > > > > > > > many
> > > > > > > > > > > > > types of extensive data, especially counts, for
> > > > example, the
> > > > > > > > > equal
> > > > > > > > > > > area
> > > > > > > > > > > > > property is desirable.  We used the Lambert
> > > > cylindrical
> > > > > > > equal
> > > > > > > > > area
> > > > > > > > > > > > > projection with standard parallels of +/- 30
> > degrees
> > > > for
> > > > > > > some
> > > > > > > > > > > western
> > > > > > > > > > > > > hemispherical work, see reference below.  (The
> > center
> > > > > > > longitude
> > > > > > > > > > > could be
> > > > > > > > > > > > > -80 west, but that is less important than the
> > choice
> > > > of
> > > > > > > > > parallels.)
> > > > > > > > > > > > >
> > > > > > > > > > > > > Before falling back on the Lambert as an easy to
> > use
> > > > > > > projection,
> > > > > > > > > I
> > > > > > > > > > > tried
> > > > > > > > > > > > > to get several ESRI products to implement an
> > > > interrupted
> > > > > > > > > projection
> > > > > > > > > > > > > using the sinusoidal projection, in part for
> > reasons
> > > > given
> > > > > > > in
> > > > > > > > > the
> > > > > > > > > > > second
> > > > > > > > > > > > > reference.  I used a separate center longitude for
> > > > north and
> > > > > > > > > south
> > > > > > > > > > > of
> > > > > > > > > > > > > the equator and the appearance is certainly more
> > > > > > > satisfactory
> > > > > > > > > than
> > > > > > > > > > > the
> > > > > > > > > > > > > Lambert in my opinion.  I'll attach a PDF of an
> > > > illustration
> > > > > > > of
> > > > > > > > > this
> > > > > > > > > > > > > approach generated in R that I hope you will get
> > but
> > > > not the
> > > > > > > > > rest of
> > > > > > > > > > > the
> > > > > > > > > > > > > list unfortunately.  I can send PDFs of the
> > references
> > > > also
> > > > > > > if
> > > > > > > > > > > needed.
> > > > > > > > > > > > >
> > > > > > > > > > > > > Denis
> > > > > > > > > > > > >
> > > > > > > > > > > > > Lawler JJ, White D, Neilson RP, Blaustein AR.
> > 2006.
> > > > > > > Predicting
> > > > > > > > > > > > > climate-induced range shifts: model differences
> > and
> > > > model
> > > > > > > > > > > reliability.
> > > > > > > > > > > > > Global Change Biology 12:1568-1584.
> > > > > > > > > > > > >
> > > > > > > > > > > > > White D.  2006.  Display of pixel loss and
> > replication
> > > > in
> > > > > > > > > > > reprojecting
> > > > > > > > > > > > > raster data from the sinusoidal projection.
> > Geocarto
> > > > > > > > > International
> > > > > > > > > > > > > 21(2):19-22.
> > > > > > > > > > > > >
> > > > > > > > > > > > > (See attached file: whemi.sinus.pdf)
> > > > > > > > > > > > >
> > > > > > > > > > > > > r-sig-geo-bounces at stat.math.ethz.ch wrote on
> > > > 2007-04-04
> > > > > > > > > 12:17:39:
> > > > > > > > > > > > >
> > > > > > > > > > > > > > Anyone know of a particularly good map
> > projection
> > > > for
> > > > > > > showing
> > > > > > > > > all
> > > > > > > > > > > of
> > > > > > > > > > > > > > North and South America without too much
> > distortion?
> > > > > > > > > > > > > >
> > > > > > > > > > > > > > THK
> > > > > > > > > > > > > >
> > > > > > > > > > > > > > --
> > > > > > > > > > > > > > Timothy H. Keitt, University of Texas at Austin
> > > > > > > > > > > > > > Contact info and schedule at
> > > > > > > http://www.keittlab.org/tkeitt/
> > > > > > > > > > > > > > Reprints at
> > http://www.keittlab.org/tkeitt/papers/
> > > > > > > > > > > > > > ODF attachment? See http://www.openoffice.org/
> > > > > > > > > > > > > >
> > > > > > > > > > > > > > _______________________________________________
> > > > > > > > > > > > > > R-sig-Geo mailing list
> > > > > > > > > > > > > > R-sig-Geo at stat.math.ethz.ch
> > > > > > > > > > > > > > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> > > > > > > > > > > > >
> > > > > > > > > > > >
> > > > > > > > > > > >
> > > > > > > > > > > > --
> > > > > > > > > > > > Timothy H. Keitt, University of Texas at Austin
> > > > > > > > > > > > Contact info and schedule at
> > > > http://www.keittlab.org/tkeitt/
> > > > > > > > > > > > Reprints at http://www.keittlab.org/tkeitt/papers/
> > > > > > > > > > > > ODF attachment? See http://www.openoffice.org/
> > > > > > > > > >
> > > > > > > > > > --
> > > > > > > > > > Roger Bivand
> > > > > > > > > > Economic Geography Section, Department of Economics,
> > > > Norwegian
> > > > > > > School
> > > > > > > > > of
> > > > > > > > > > Economics and Business Administration, 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
> > > > > > > > > >
> > > > > > > > > > _______________________________________________
> > > > > > > > > > R-sig-Geo mailing list
> > > > > > > > > > R-sig-Geo at stat.math.ethz.ch
> > > > > > > > > > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> > > > > > > > >
> > > > > > > >
> > > > > > > >
> > > > > > > > --
> > > > > > > > Timothy H. Keitt, University of Texas at Austin
> > > > > > > > Contact info and schedule at http://www.keittlab.org/tkeitt/
> > > > > > > > Reprints at http://www.keittlab.org/tkeitt/papers/
> > > > > > > > ODF attachment? See http://www.openoffice.org/
> > > > > > > >
> > > > > > > > _______________________________________________
> > > > > > > > R-sig-Geo mailing list
> > > > > > > > R-sig-Geo at stat.math.ethz.ch
> > > > > > > > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> > > > > > >
> > > > > > > _______________________________________________
> > > > > > > R-sig-Geo mailing list
> > > > > > > R-sig-Geo at stat.math.ethz.ch
> > > > > > > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> > > > > > >
> > > > > >
> > > > > > --
> > > > > > Roger Bivand
> > > > > > Economic Geography Section, Department of Economics, Norwegian
> > > > School of
> > > > > > Economics and Business Administration, 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
> > > > > >
> > > > > >
> > > > > >
> > > > >
> > > > >
> > > > > --
> > > > > Timothy H. Keitt, University of Texas at Austin
> > > > > Contact info and schedule at http://www.keittlab.org/tkeitt/
> > > > > Reprints at http://www.keittlab.org/tkeitt/papers/
> > > > > ODF attachment? See http://www.openoffice.org/
> > >
> > > --
> > > Roger Bivand
> > > Economic Geography Section, Department of Economics, Norwegian School
> > of
> > > Economics and Business Administration, 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
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, 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