Dear R-sig-geoers,
I am stuck with a problem and i was hoping if you could help me.
Suppose i am given a coordinate in lat/long (say 25 lat 35 long) and a
great circle distance (say 10km). Also assume that i am given an angle
with respect to latitude (an angle from the x-axis say 30 degrees).
What would be the best way to determine a coordinate 10kms away from
25 lat 35 long with a given angle (30 degrees)?
I think in euclidean space, I can use tangent of the given angle
and combine with the Pythagorean theorem to find the target
coordinate.
Thank you very much in advance,
Harry
Determining a coordinate given another coordinate, a great circle distance, and an angle
5 messages · Tomislav Hengl, Robert J. Hijmans, Harry Kim
Hi Harry, The key issue is that you need to select the (Euclidean) coordinate system. For example, to represent the whole world, there are many possiblities (http://www.radicalcartography.net/?projectionref). If you are working with the great circle distance, make sure that the geographic coordinate system uses a sphere. Say if you want to use some equidistant coordinate system with center at 0,0 longlat (Lambert Azimuthal Equal-Area), you could then get the coordinates by using:
P1 <- data.frame(E=35, N=25)
lambda <- 60 # azimuth!
dist12 <- 300000
coordinates(P1) <- ~E+N
proj4string(P1) <- CRS("+proj=longlat +datum=WGS84")
P1.xy <- spTransform(P1, CRS("+proj=laea +lat_0 +lon_0 +x_0 +y_0"))
P2.xy <- data.frame(E=P1.xy at coords[1,1]+dist12*sin(lambda),
N=P1.xy at coords[1,2]+dist12*cos(lambda))
coordinates(P2.xy) <- ~E+N
proj4string(P2.xy) <- P1.xy at proj4string
P2 <- spTransform(P2.xy, CRS("+proj=longlat +datum=WGS84"))
P12 <- rbind(P1, P2)
P12
SpatialPoints:
E N
[1,] 35.00000 25.000
[2,] 33.55488 22.554
Coordinate Reference System (CRS) arguments:
+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
load(file=url("http://spatial-analyst.net/DATA/worldborders.RData"))
spplot(worldborders["NAME"], colorkey=FALSE, col.regions=rep(grey(0.5),
length(levels(worldborders$NAME))), sp.layout=list("sp.points", P12, pch="+", col="red", cex=3))
But I would consider using the UTM or which ever coordinate system is of your interest; then you
also need to convert the GCD to distance in the local system. For example, to represent the whole of
USA (contiguous 48-state area) I typically use:
+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83
+units=m +no_defs
See also:
http://spatial-analyst.net/DATA/usgrids5km.zip
cheers,
T. Hengl
http://spatial-accuracy.org/FromGEOSTAT2009
PS: I never got any photos from your trip to the Pyramids!
-----Original Message-----
From: r-sig-geo-bounces at stat.math.ethz.ch [mailto:r-sig-geo-bounces at stat.math.ethz.ch] On Behalf
Of Harry Kim
Sent: Tuesday, October 06, 2009 7:23 PM
To: r-sig-geo at stat.math.ethz.ch
Subject: [R-sig-Geo] Determining a coordinate given another coordinate,a great circle distance,
and an angle
Dear R-sig-geoers,
I am stuck with a problem and i was hoping if you could help me.
Suppose i am given a coordinate in lat/long (say 25 lat 35 long) and a
great circle distance (say 10km). Also assume that i am given an angle
with respect to latitude (an angle from the x-axis say 30 degrees).
What would be the best way to determine a coordinate 10kms away from
25 lat 35 long with a given angle (30 degrees)?
I think in euclidean space, I can use tangent of the given angle
and combine with the Pythagorean theorem to find the target
coordinate.
Thank you very much in advance,
Harry
_______________________________________________ R-sig-Geo mailing list R-sig-Geo at stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
An embedded and charset-unspecified text was scrubbed... Name: not available URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20091007/f5276aff/attachment.pl>
Hey Tom!
It's great to hear from you :). I did think about using a UTM
projection but the region of interest (Turkey) is too big. It requires
6 UTM zones to cover it and if the point is around the boundary i
thought that it would create some problems. I wasn't aware of other
distance/area preserving projections. I will take a look at the zip
file on your website. For some reason radicalcartography website isn't
working for me.
Thanks again for the detailed code tom!
Hope things are going well for you.
P.S: The pyramids were awesome but i thought Petra was better :-D.
http://www.facebook.com/cardboy?ref=profile#/album.php?aid=2477670&id=1204683&op=6
http://www.facebook.com/cardboy?ref=profile#/album.php?aid=2481932&id=1204683&op=6
On Wed, Oct 7, 2009 at 1:36 AM, Tomislav Hengl
<hengl at spatial-analyst.net> wrote:
Hi Harry, The key issue is that you need to select the (Euclidean) coordinate system. For example, to represent the whole world, there are many possiblities (http://www.radicalcartography.net/?projectionref). If you are working with the great circle distance, make sure that the geographic coordinate system uses a sphere. Say if you want to use some equidistant coordinate system with center at 0,0 longlat (Lambert Azimuthal Equal-Area), you could then get the coordinates by using:
P1 <- data.frame(E=35, N=25)
lambda <- 60 ?# azimuth!
dist12 <- 300000
coordinates(P1) <- ~E+N
proj4string(P1) <- CRS("+proj=longlat +datum=WGS84")
P1.xy <- spTransform(P1, CRS("+proj=laea +lat_0 +lon_0 +x_0 +y_0"))
P2.xy <- data.frame(E=P1.xy at coords[1,1]+dist12*sin(lambda),
N=P1.xy at coords[1,2]+dist12*cos(lambda))
coordinates(P2.xy) <- ~E+N
proj4string(P2.xy) <- P1.xy at proj4string
P2 <- spTransform(P2.xy, CRS("+proj=longlat +datum=WGS84"))
P12 <- rbind(P1, P2)
P12
SpatialPoints: ? ? ? ? ? ?E ? ? ?N [1,] 35.00000 25.000 [2,] 33.55488 22.554 Coordinate Reference System (CRS) arguments: +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
load(file=url("http://spatial-analyst.net/DATA/worldborders.RData"))
spplot(worldborders["NAME"], colorkey=FALSE, col.regions=rep(grey(0.5),
length(levels(worldborders$NAME))), sp.layout=list("sp.points", P12, pch="+", col="red", cex=3))
But I would consider using the UTM or which ever coordinate system is of your interest; then you
also need to convert the GCD to distance in the local system. For example, to represent the whole of
USA (contiguous 48-state area) I typically use:
+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83
+units=m +no_defs
See also:
http://spatial-analyst.net/DATA/usgrids5km.zip
cheers,
T. Hengl
http://spatial-accuracy.org/FromGEOSTAT2009
PS: I never got any photos from your trip to the Pyramids!
-----Original Message----- From: r-sig-geo-bounces at stat.math.ethz.ch [mailto:r-sig-geo-bounces at stat.math.ethz.ch] On Behalf Of Harry Kim Sent: Tuesday, October 06, 2009 7:23 PM To: r-sig-geo at stat.math.ethz.ch Subject: [R-sig-Geo] Determining a coordinate given another coordinate,a great circle distance, and an angle Dear R-sig-geoers, ? ? ? I am stuck with a problem and i was hoping if you could help me. Suppose i am given a coordinate in lat/long (say 25 lat 35 long) and a great circle distance (say 10km). Also assume that i am given an angle with respect to latitude (an angle from the x-axis say 30 degrees). What would be the best way to determine a coordinate 10kms away from 25 lat 35 long with a given angle (30 degrees)? ? ? ?I think in euclidean space, I can use tangent of the given angle and combine with the Pythagorean theorem to find the target coordinate. Thank you very much in advance, Harry
_______________________________________________ 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
Hi Robert,
Thank you so much! I just tested the code and it works
beautifully. This is exactly what I wanted to do. I had a feeling
that it could be calculated analytically but I kept messing up my
calculations.
Thanks again!
Harry
On Wed, Oct 7, 2009 at 10:47 AM, Robert J. Hijmans <r.hijmans at gmail.com> wrote:
Harry, I think you are after something like the function below (I did not test it much). Robert # based on code at http://www.movable-type.co.uk/scripts/latlong.html # ? 2002-2009 Chris Veness # license: LPGL; without any warranty express or implied # R version by R Hijmans # calculate destination point given start point, initial bearing (deg) and distance (km) # ? see http://williams.best.vwh.net/avform.htm#LL destPoint <- function(lon, lat, bearing, d, dms=FALSE, R=6378137) { # bearing in degrees # d in meters # R=earth's mean radius in m RadDeg <- pi / 180 lat1 = lat * RadDeg lon1 = lon * RadDeg brng = bearing * RadDeg lat2 <- asin( sin(lat1) * cos(d/R) + cos(lat1) * sin(d/R) * cos(brng) ) lon2 <- lon1 + atan2(sin(brng) * sin(d/R) * cos(lat1), cos(d/R) - sin(lat1) * sin(lat2)) lon2 <- (lon2 + pi) %% ( 2 * pi) - pi; # normalise to -180...+180 if (is.nan(lat2) | is.nan(lon2)) return(NULL) return ( (cbind(lon2, lat2)) / RadDeg ) } destPoint(35, 25, 30, 100000) On Tue, Oct 6, 2009 at 10:22 AM, Harry Kim <harryk at cal.berkeley.edu> wrote:
Dear R-sig-geoers, ? ? ?I am stuck with a problem and i was hoping if you could help me. Suppose i am given a coordinate in lat/long (say 25 lat 35 long) and a great circle distance (say 10km). Also assume that i am given an angle with respect to latitude (an angle from the x-axis say 30 degrees). What would be the best way to determine a coordinate 10kms away from 25 lat 35 long with a given angle (30 degrees)? ? ? I think in euclidean space, I can use tangent of the given angle and combine with the Pythagorean theorem to find the target coordinate. Thank you very much in advance, Harry
_______________________________________________ R-sig-Geo mailing list R-sig-Geo at stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo