<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0 Transitional//EN">
<HTML><HEAD>
<META http-equiv=Content-Type content="text/html; charset=iso-8859-1">
<META content="MSHTML 6.00.2800.1170" name=GENERATOR>
<STYLE></STYLE>
</HEAD>
<BODY>
<DIV><FONT face=Arial size=2>Charlie,</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2>I checked the great circle stuff, you are right
that this was not the most accurate although difference generally within
1%<BR></FONT></DIV>
<DIV><FONT face=Arial size=2>I replaced it with a shortened version of what you
were using:</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2><EM>Private Function GeoDistance(Lon1 As Double,
Lat1 As Double, Lon2 As Double, Lat2 As Double) As Double</EM></FONT></DIV>
<DIV><FONT face=Arial size=2><EM></EM></FONT> </DIV>
<DIV><FONT face=Arial size=2><EM>'The great circle distance d between two points
with coordinates {lat1,lon1} and {lat2,lon2} is given by:<BR>'<BR>'d =
acos(Sin(Lat1) * Sin(Lat2) + Cos(Lat1) * Cos(Lat2) * Cos(Lon1 - Lon2))<BR>'A
mathematically equivalent formula, which is less subject to rounding error for
short distances is:<BR>'<BR>'d=2*asin(sqrt((sin((lat1-lat2)/2))^2
+<BR>'
cos(lat1)*cos(lat2)*(sin((lon1-lon2)/2))^2))<BR>'<BR>' an additional check is
put in to make sure the argument to asin is within the range
0-1</EM></FONT></DIV>
<DIV><FONT face=Arial size=2><EM></EM></FONT> </DIV>
<DIV><FONT face=Arial size=2><EM>Const deg2rad =
1.74532925199433E-02
' radians per degree<BR>Const Pi =
3.14159265358979
' Pi<BR>Const Radius = 6377000<BR>Dim dLat As Double, dLon As
Double</EM></FONT></DIV>
<DIV><FONT face=Arial size=2><EM></EM></FONT> </DIV>
<DIV><FONT face=Arial size=2><EM>dLat = Lat2 - Lat1<BR>dLon = Lon2 -
Lon1</EM></FONT></DIV>
<DIV><FONT face=Arial size=2><EM></EM></FONT> </DIV>
<DIV><FONT face=Arial size=2><EM>GeoDistance = Radius * 2 * asin(min(1,
Sqr(Sin(deg2rad * dLat / 2) ^ 2 + Cos(deg2rad * Lat1) * Cos(deg2rad * Lat2) *
Sin(deg2rad * dLon / 2) ^ 2)))</EM></FONT></DIV>
<DIV><FONT face=Arial size=2><EM></EM></FONT> </DIV>
<DIV><FONT face=Arial size=2><EM>End Function</EM></FONT></DIV>
<DIV><FONT face=Arial size=2><EM></EM></FONT> </DIV>
<DIV><FONT face=Arial size=2><EM>Private Function asin(ByVal x As Double) As
Double<BR> asin = Atn(x / Sqr(-x * x + 1))<BR>End
Function</EM></FONT></DIV>
<DIV><FONT face=Arial size=2><EM></EM></FONT> </DIV>
<DIV><FONT face=Arial size=2><EM>Private Function min(ByVal x As Double, ByVal y
As Double)<BR> If x < y Then min = x Else min = y<BR>End
Function<BR></EM></FONT></DIV>
<DIV><FONT face=Arial size=2>The differences will be neglectable, i chose not to
make my own sin2 function like you have, for simplicity's sake...</FONT></DIV>
<DIV><FONT face=Arial size=2>Now if we want do do this completely right the
first time, i adapted the rest of it to exactly compute one pixel in the middle
rather then 10% of the window, just to stick to the openGIS definition
exactly:</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2><EM> With
mvarDisplay.Frame<BR> mapwidth = .xMax -
.xMin<BR> mapheight = .yMax - .yMin<BR> End
With<BR> <BR> With mvarDisplay.ZoomExtent<BR> '
extent of central pixel<BR> x1 = (.xMin + .xMax) /
2<BR> x2 = x1 + (.xMax - .xMin) /
mapwidth<BR> y1 = (.yMin = .yMax) / 2<BR> y2
= y1 + (.yMax - .yMin) / mapheight<BR> If mvarUnitMeter = 0
Then<BR> ' great circle distance for
lat-lon<BR> pixelmeters = GeoDistance(x1, y1, x2,
y2)<BR> Else<BR> ' pythagoras
for projected maps in meters<BR> pixelmeters =
Sqr((mvarUnitMeter * (.xMax - .xMin)) ^ 2 + (mvarUnitMeter * (.yMax - .yMin)) ^
2)<BR> End If<BR> <BR> End
With<BR> DoGetZoomFactor = pixelmeters / 0.00025<BR></EM></FONT></DIV>
<DIV><FONT face=Arial size=2>I'll try to get a server running here today, but
first i have to add the scalehints to the XML output.</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2>Best regards,</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2>Bart Adriaanse</FONT></DIV>
<DIV><FONT face=Arial size=2>Demis bv</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2>----- Original Message ----- </FONT></DIV>
<DIV><FONT face=Arial size=2>From: "Paul Kelly" <</FONT><A
href="mailto:paul-grass@stjohnspoint.co.uk"><FONT face=Arial
size=2>paul-grass@stjohnspoint.co.uk</FONT></A><FONT face=Arial
size=2>></FONT></DIV>
<DIV><FONT face=Arial size=2>To: "H Bowman" <</FONT><A
href="mailto:hamish_nospam@yahoo.com"><FONT face=Arial
size=2>hamish_nospam@yahoo.com</FONT></A><FONT face=Arial
size=2>></FONT></DIV>
<DIV><FONT face=Arial size=2>Cc: <</FONT><A
href="mailto:grass5@grass.itc.it"><FONT face=Arial
size=2>grass5@grass.itc.it</FONT></A><FONT face=Arial size=2>>; <</FONT><A
href="mailto:osrs-proj@remotesensing.org"><FONT face=Arial
size=2>osrs-proj@remotesensing.org</FONT></A><FONT face=Arial
size=2>></FONT></DIV>
<DIV><FONT face=Arial size=2>Sent: Wednesday, July 23, 2003 13:09</FONT></DIV>
<DIV><FONT face=Arial size=2>Subject: [OSRS-PROJ] Re: [GRASS5] adding a new
datum/projection?</FONT></DIV>
<DIV><FONT face=Arial><BR><FONT size=2></FONT></FONT></DIV><FONT face=Arial
size=2>> Hello<BR>> <BR>> On Wed, 23 Jul 2003, H Bowman wrote:<BR>>
<BR>> > Hi again --<BR>> ><BR>> ><BR>> > I've got some
new data in yet another new projection, which uses its own<BR>> > new
datum.<BR>> ><BR>> ><BR>> > The datum is New Zealand Geodetic
Datum 2000 (NZGD2000), which is based<BR>> > on the International
Terrestrial Reference Frame 1996 (ITRF96). This is<BR>> > pretty much the
same as WGS84 but uses the GRS80 ellipsoid (slightly<BR>> > different
flattening parameter).<BR>> ><BR>> > So as far as I understand:
NZGD2k (==ITRF96) is the WGS84 datum combined<BR>> > with the GRS80
ellipsoid.<BR>> <BR>> We are pushing at the limits of the accuracy of the
current GRASS/PROJ<BR>> system here. I looked at the document and to me it
seems that NZGD2000 is<BR>> based on the ITRS datum. It is 'realised' (which
I think means a<BR>> description of the current location of the earth's
plates based on<BR>> measurements from a number of stations around the world)
by several ITRF<BR>> defintions, and each can be realised at any epoch (point
in time).<BR>> <BR>> For high accuracy you would define a shift that
related each epoch to WGS84<BR>> (would change because of plate tectonics)
and each would have a separate<BR>> datum entry in GRASS. It would get out of
hand very quickly if these were<BR>> all included but it may be possible to
use GRASS in this way if you could<BR>> get figures for the shift (for New
Zealand you need ITRF96 at 2000.0<BR>> epoch). </FONT><A
href="http://www.iers.org/iers/products/itrf/"><FONT face=Arial
size=2>http://www.iers.org/iers/products/itrf/</FONT></A><FONT face=Arial
size=2> looks relevant but I don't<BR>> really know.<BR>> <BR>> >
Except g.setproj won't let me do that, either with a new location or<BR>>
> retroactively. The old 5.0.2 version asks for ellipsoid first, then<BR>>
> datum, but complains when you try to mix wgs84 and grs80 and quits<BR>>
> without writing anything. The latest g.setproj never gives you a
chance<BR>> > to set the ellipsoid, it is chosen from your datum.<BR>>
<BR>> No it doesn't seem logical to do that. GRASS uses a very limited
definition<BR>> of a datum which really just comes down to the ellipsoid
used, so by<BR>> changing the ellipsoid you are really taking away everything
the datum<BR>> stands for.<BR>> <BR>> > I'm not sure if I get the
correct result if I just change the ellps:<BR>> > line to grs80 in the
PROJ_INFO file?<BR>> <BR>> ellps: grs80<BR>> towgs84: 0,0,0<BR>>
<BR>> and make sure there is no 'datum' line.<BR>> <BR>> > Is
anything else based on ITRF96? Would it be worth adding that as a<BR>> >
common denominator datum? I think the Map Grid of Australia might be in<BR>>
> a similar situation (??).<BR>> <BR>> It might be worth adding ITRS
(International Terrestrial Reference<BR>> *System*), but to be compatible
with other GIS we should probably<BR>> explicitly add
New_Zealand_Geodetic_Datum_2000 and the others. But as<BR>> they are the same
as far as the level of accuracy GRASS/PROJ currently<BR>> provides goes, it
seems like we could end up with a lot of clutter.<BR>> However I have already
done that for several countries in Europe which use<BR>> the international
ellipsoid but have different datum names, so there is no<BR>> point in
discriminating....<BR>> <BR>> <BR>> ><BR>> > see:
(102k)<BR>> > </FONT><A
href="http://www.linz.govt.nz/rcs/linz/17880/difference_wgs84_nzgd2000.pdf"><FONT
face=Arial
size=2>http://www.linz.govt.nz/rcs/linz/17880/difference_wgs84_nzgd2000.pdf</FONT></A><BR><FONT
face=Arial size=2>> ><BR>> ><BR>> ><BR>> > As for the
new projection, it is New Zealand Transverse Mercator (NZTM),<BR>> > which
is in terms of that new NZGD2k datum.<BR>> > Details:<BR>> > Datum:
NZGD2k<BR>> > Origin Latitude: 0° South<BR>> > Origin Longitude:
173° East<BR>> > False Northing: 10 000 000m N<BR>> > False Easting:
1 600 000m E<BR>> > Scale Factor: 0.9996<BR>> ><BR>> > Setting
projection as tmerc and entering those terms goes smoothly.<BR>> > Should
NZTM get its own projection entry or is entering by hand with the<BR>> >
tmerc projection every time the way to go?<BR>> <BR>> Yes that is the way
to do it until GRASS has co-ordinate system support,<BR>> which I don't see
happening any time soon.. you have it very easy with<BR>> New Zealand Map
Grid which puts in all the parameters for you.<BR>> <BR>> ><BR>>
> for more details see: (44k)<BR>> > </FONT><A
href="http://www.linz.govt.nz/rcs/linz/5684/nztransverse_mercator.pdf"><FONT
face=Arial
size=2>http://www.linz.govt.nz/rcs/linz/5684/nztransverse_mercator.pdf</FONT></A><BR><FONT
face=Arial size=2>> ><BR>> > That PDF also has a page of useful
ellipsoid to grid formulae which<BR>> > might be useful to
someone.<BR>> ><BR>> ><BR>> > It doesn't help much that half
the country is moving 5cm/year in the<BR>> > opposite direction to the
other half..<BR>> <BR>> It would be useful to keep an eye out for any Free
Software that handles<BR>> that level of accuracy<BR>> </FONT><A
href="http://www-gpsg.mit.edu/~simon/gtgk/"><FONT face=Arial
size=2>http://www-gpsg.mit.edu/~simon/gtgk/</FONT></A><FONT face=Arial size=2>
has some interesting links<BR>> <BR>> ><BR>> ><BR>> >
thanks for any insight,<BR>> <BR>> Well hopefully someone will correct the
inevitable errors and<BR>> misunderstandings in what I've written
above<BR>> <BR>> Paul<BR>>
----------------------------------------<BR>> PROJ.4 Discussion List<BR>>
See </FONT><A href="http://www.remotesensing.org/proj"><FONT face=Arial
size=2>http://www.remotesensing.org/proj</FONT></A><FONT face=Arial size=2> for
subscription, unsubscription<BR>> and other information.</FONT></BODY></HTML>