Renaud Lancelot
2002-May-24 18:20 UTC
[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 _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
White.Denis@epamail.epa.gov
2002-May-24 19:56 UTC
[R] intersecting polygons and conversion from decimal degree to km
> 1. How can I compute the intersecting area between 2 polygons ?The general case (with interior exclusions and/or concave borders) is a little work. I don't know whether there are any solutions in R, but many GIS have this feature.> 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) ?Here is a very simple map projection that you could use. It ignores the slightly non-spherical shape of the earth, and is neither equal-area nor conformal. After converting your geographic coordinates to kilometers with this function, you can use function "areapl" in package "splancs" to calculate the area. marinus <- function (lat, lon) { # Arguments are vectors of latitude and longitude, # that is, geographic coordinates in decimal degrees. # Projection center defined as midpoint in lat-lon space. # Returns a two column matrix with column names "x" and "y", # in units of kilometers. This is the equidistant cylindric # map projection, here named after Marinus of Tyre # (see JP Snyder. USGS Prof Paper 1395, p 90). R <- 6371 rlat <- range (lat, na.rm=TRUE) rlon <- range (lon, na.rm=TRUE) midlat <- diff (rlat) / 2 + rlat[1] midlon <- diff (rlon) / 2 + rlon[1] x <- R * (lon-midlon)*pi/180 * cos (midlat*pi/180) y <- R * (lat-midlat)*pi/180 xy <- cbind (x=x, y=y) xy } -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- 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 _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
Dr. M.Sawada
2002-May-24 22:51 UTC
[R] intersecting polygons and conversion from decimal degree to km
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
White.Denis@epamail.epa.gov
2002-May-24 23:51 UTC
[R] intersecting polygons and conversion from decimal degree to km
Dr. M.Sawada wrote:> 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) ...At the equator, 1 degree of longitude is 111+ km on the Clarke 1866 ellipsoid. Choosing a map projection depends on the problem. Equal-area may not be preferable to conformal or equidistant. And simplicity may be fine depending on precision requirements. For background on map projections try a web search, e.g. http://www.colorado.edu/geography/gcraft/notes/mapproj/mapproj_f.html. There have been a number of discussions on r-help on map projections. See threads about 20 April 2001, and 22 June 2001. Denis White US EPA, 200 SW 35th St, Corvallis, Oregon, 97333 USA voice: 541.754.4476, email: white.denis at epa.gov web: www.epa.gov/wed/pages/staff/white/ -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- 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 _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
Roger Peng
2002-May-25 00:58 UTC
[R] intersecting polygons and conversion from decimal degree to km
Intersecting polygons is a fairly non-trivial task, one which I was confronted with a few years ago. Luckily, Alan Murta has written a very nice C library which can be downloaded at http://www.cs.man.ac.uk/aig/staff/alan/software/ A little while back I wrote an R interface to this library which is fairly straightforward to use but it's a little rough and ready. The good thing is that Murta's library is *very* fast. It can intersect/difference/union polygons with over a thousand points very quickly. If you're interested in my R package please email me privately. -roger _______________________________ UCLA Department of Statistics rpeng at stat.ucla.edu http://www.stat.ucla.edu/~rpeng On Fri, 24 May 2002, Renaud Lancelot wrote:> 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 _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
I have a very much simpler problem - radiotracking data, data limited and grubby so anything more than simple convex polygons is unwarranted (obtained using chull()). Distances so small we can regard the Earth as flat and Euclidean. Has anyone a smart idea to detect the (possible) intersection points of two such bounding polygons? (once the intersection points are detected then the overlap area is another simple convex polygon), Richard Rowe Richard Rowe Senior Lecturer Department of Zoology and Tropical Ecology, James Cook University Townsville, Queensland 4811, Australia fax (61)7 47 25 1570 phone (61)7 47 81 4851 e-mail: Richard.Rowe at jcu.edu.au http://www.jcu.edu.au/school/tbiol/zoology/homepage.html -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- 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 _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._