[R] intersecting polygons and conversion from decimal degree to km

[Dr. M.Sawada]

Renaud:
Here is a simple Albers-Equal Area Conic program in R that you can use
(actually I wrote it for S+ but should be the same - you may have to
modify the dimnames() function for R).  This is preferred over an
equidistant projection.  Long/Lat, as you are probably aware, are a
spherical system and so distances are not isotropic (1 degree of long at
the equator is 116 miles vs 0 miles at the pole) and euclidean geometry
does not apply strictly (length x width <> area).  Thus it would be more
difficult to calculate area and you would need to resort to spherical
trig.  It is easier to use an equal area projection like this one with
standard lines in the bottom 1/3 and top 2/3 of your study area.  Modify
R in the code which represents the radius of the earth in km - here
multiplied by 1000 to provide a coordinate system in meters.  If you
want to calculate areas on an ellipsoidal earth (e.g., WGS84) you will
need to use a GIS.

albers.mat.f<-function(xy, origin = c(0, 0), stdlns = c(0, 0)) { 
#inputs############# 
#xy - a two column dataframe or matrix with longitudes in the first
column (x) and latitudes in the second column 
#(y) 
#origin = c(0,0): is the Longitude and Latitude for the center of the
output coordinate system e.g., c(-100,50) 
#	for -100 West and 50 N.
#stdlns = c(0,0) are the standard lines (parallels) where true scale is
found e.g., for North America perhaps
#c(35,80) for 35 N and 80 N since this distributes distortion throughout
the map. #R - Earth Radius #phi1 - Lat; Stdln 1 #phi2 - Lat; Stdln 2
#phi0 - Lat; Origin #lambda0 - Lon; Origin #xy - c(phi,lambda); point to
be projected #C = (cos(phi1))^2+2nsin(phi1) #n=(sin(phi1)+sin(phi2))/2
#theta = n(lambda-lambda0) #p = (R(C-2n sin(phi))^.5)/n #pnot = (R(C-2n
sin (phi0))^.5)/n #x = p sin(theta) #y = pnot - p cos(theta)
	xy <- as.matrix(xy)
	dimnames(xy) <- list(NULL, c("x", "y"))	
	xy <- (xy * pi)/180
	R <- 6378 * 1000
	phi1 <- (stdlns[1] * pi)/180
	phi2 <- (stdlns[2] * pi)/180
	phi0 <- (origin[2] * pi)/180
	lambda0 <- (origin[1] * pi)/180	
	n <- (sin(phi1) + sin(phi2))/2
	twon <- 2 * n
	C <- ((cos(phi1))^2 + 2 * n * sin(phi1))
	theta <- n * (xy[, 1] - lambda0)
	xy[, 1] <- ((R * (C - twon * sin(xy[, 2]))^0.5)/n) * sin(theta)
	xy[, 2] <- ((R * (C - twon * sin(phi0))^0.5)/n) - ((R * (C - 2 *
n * sin(xy[, 2]))^0.5)/n) * cos(theta)
	xy
}

I've implemented the Hodgman-Sutherland Algorithm to clip polygons to
graph boxes (S+ won't clip filled polygons) as a DLL interfaced with S+
and you can see it at:

http://www.uottawa.ca/academic/arts/geographie/lpcweb/newlook/data_and_d
ownloads/download/pclip/pclip.htm

You might be able to take this approach and re-implement it for
intersecting two convex polygons but it will be some work.  I don't know
of any intersection algorithms in R either - use ArcView or ArcGIS or
Mapinfo.  You could also call the WinAPI within a DLL to do intersection
areas.

Cheers,
Mike

______________________________________________
Dr. M.Sawada
Assistant Professor, GIS
University of Ottawa
Dept. Geography
P.O. Box 450 Stn. A
Ottawa, ON K1N 6N9
Tel: 613-562-5800 x1040
Fax: 613-562-5145
Eml: msawada at uottawa.ca 

-----Original Message-----
From: owner-r-help at stat.math.ethz.ch
[mailto:owner-r-help at stat.math.ethz.ch] On Behalf Of Renaud Lancelot
Sent: Friday, May 24, 2002 2:20 PM
To: R-Help
Subject: [R] intersecting polygons and conversion from decimal degree to
km


Dear all,

1. How can I compute the intersecting area between 2 polygons ?

2. I have polygons with coordinates in decimal degrees (i.e. 13 deg 30
min = 13.5 decimal degrees). I want to compute their area and get the
results in square meters or square kiometers. Can anyone give me a
conversion coefficient or a pointer where I can find this information
(sorry for this off topic question) ?

Thanks in advance,

Renaud

-- 
Dr Renaud Lancelot, vétérinaire
CIRAD, Département Elevage et Médecine Vétérinaire (CIRAD-Emvt)
Programme Productions Animales
http://www.cirad.fr/presentation/programmes/prod-ani.shtml (Français)
http://www.cirad.fr/presentation/en/program-eng/prod-ani.shtml (English)

ISRA-LNERV                      tel    (221) 832 49 02
BP 2057 Dakar-Hann              fax    (221) 821 18 79 (CIRAD)
Senegal                         e-mail renaud.lancelot at cirad.fr
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.
-.-.-.-
r-help mailing list -- Read
http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._.
_._._._

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._