[OSRS-PROJ] Re: [GRASS5] adding a new datum/projection?

Wed Jul 23 09:06:16 EDT 2003

Charlie,

I checked the great circle stuff, you are right that this was not the most accurate although difference generally within 1%

I replaced it with a shortened version of what you were using:

Private Function GeoDistance(Lon1 As Double, Lat1 As Double, Lon2 As Double, Lat2 As Double) As Double

'The great circle distance d between two points with coordinates {lat1,lon1} and {lat2,lon2} is given by:
'
'd = acos(Sin(Lat1) * Sin(Lat2) + Cos(Lat1) * Cos(Lat2) * Cos(Lon1 - Lon2))
'A mathematically equivalent formula, which is less subject to rounding error for short distances is:
'
'd=2*asin(sqrt((sin((lat1-lat2)/2))^2 +
'                 cos(lat1)*cos(lat2)*(sin((lon1-lon2)/2))^2))
'
' an additional check is put in to make sure the argument to asin is within the range 0-1

Const deg2rad = 1.74532925199433E-02                ' radians per degree
Const Pi = 3.14159265358979                         ' Pi
Const Radius = 6377000
Dim dLat As Double, dLon As Double

dLat = Lat2 - Lat1
dLon = Lon2 - Lon1

GeoDistance = Radius * 2 * asin(min(1, Sqr(Sin(deg2rad * dLat / 2) ^ 2 + Cos(deg2rad * Lat1) * Cos(deg2rad * Lat2) * Sin(deg2rad * dLon / 2) ^ 2)))

End Function

Private Function asin(ByVal x As Double) As Double
  asin = Atn(x / Sqr(-x * x + 1))
End Function

Private Function min(ByVal x As Double, ByVal y As Double)
  If x < y Then min = x Else min = y
End Function

The differences will be neglectable, i chose not to make my own sin2 function like you have, for simplicity's sake...
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:

  With mvarDisplay.Frame
    mapwidth = .xMax - .xMin
    mapheight = .yMax - .yMin
  End With

  With mvarDisplay.ZoomExtent
    ' extent of central pixel
    x1 = (.xMin + .xMax) / 2
    x2 = x1 + (.xMax - .xMin) / mapwidth
    y1 = (.yMin = .yMax) / 2
    y2 = y1 + (.yMax - .yMin) / mapheight
    If mvarUnitMeter = 0 Then
      ' great circle distance for lat-lon
      pixelmeters = GeoDistance(x1, y1, x2, y2)
    Else
      ' pythagoras for projected maps in meters
      pixelmeters = Sqr((mvarUnitMeter * (.xMax - .xMin)) ^ 2 + (mvarUnitMeter * (.yMax - .yMin)) ^ 2)
    End If

  End With
  DoGetZoomFactor = pixelmeters / 0.00025

I'll try to get a server running here today, but first i have to add the scalehints to the XML output.

Best regards,

Bart Adriaanse
Demis bv

----- Original Message ----- 
From: "Paul Kelly" <paul-grass at stjohnspoint.co.uk>
To: "H Bowman" <hamish_nospam at yahoo.com>
Cc: <grass5 at grass.itc.it>; <osrs-proj at remotesensing.org>
Sent: Wednesday, July 23, 2003 13:09
Subject: [OSRS-PROJ] Re: [GRASS5] adding a new datum/projection?

> Hello
> 
> On Wed, 23 Jul 2003, H Bowman wrote:
> 
> > Hi again --
> >
> >
> > I've got some new data in yet another new projection, which uses its own
> > new datum.
> >
> >
> > The datum is New Zealand Geodetic Datum 2000 (NZGD2000), which is based
> > on the International Terrestrial Reference Frame 1996 (ITRF96). This is
> > pretty much the same as WGS84 but uses the GRS80 ellipsoid (slightly
> > different flattening parameter).
> >
> > So as far as I understand: NZGD2k (==ITRF96) is the WGS84 datum combined
> > with the GRS80 ellipsoid.
> 
> We are pushing at the limits of the accuracy of the current GRASS/PROJ
> system here. I looked at the document and to me it seems that NZGD2000 is
> based on the ITRS datum. It is 'realised' (which I think means a
> description of the current location of the earth's plates based on
> measurements from a number of stations around the world) by several ITRF
> defintions, and each can be realised at any epoch (point in time).
> 
> For high accuracy you would define a shift that related each epoch to WGS84
> (would change because of plate tectonics) and each would have a separate
> datum entry in GRASS. It would get out of hand very quickly if these were
> all included but it may be possible to use GRASS in this way if you could
> get figures for the shift (for New Zealand you need ITRF96 at 2000.0
> epoch). http://www.iers.org/iers/products/itrf/ looks relevant but I don't
> really know.
> 
> > Except g.setproj won't let me do that, either with a new location or
> > retroactively. The old 5.0.2 version asks for ellipsoid first, then
> > datum, but complains when you try to mix wgs84 and grs80 and quits
> > without writing anything. The latest g.setproj never gives you a chance
> > to set the ellipsoid, it is chosen from your datum.
> 
> No it doesn't seem logical to do that. GRASS uses a very limited definition
> of a datum which really just comes down to the ellipsoid used, so by
> changing the ellipsoid you are really taking away everything the datum
> stands for.
> 
> > I'm not sure if I get the correct result if I just change the ellps:
> > line to grs80 in the PROJ_INFO file?
> 
> ellps: grs80
> towgs84: 0,0,0
> 
> and make sure there is no 'datum' line.
> 
> > Is anything else based on ITRF96? Would it be worth adding that as a
> > common denominator datum? I think the Map Grid of Australia might be in
> > a similar situation (??).
> 
> It might be worth adding ITRS (International Terrestrial Reference
> *System*), but to be compatible with other GIS we should probably
> explicitly add New_Zealand_Geodetic_Datum_2000 and the others. But as
> they are the same as far as the level of accuracy GRASS/PROJ currently
> provides goes, it seems like we could end up with a lot of clutter.
> However I have already done that for several countries in Europe which use
> the international ellipsoid but have different datum names, so there is no
> point in discriminating....
> 
> 
> >
> > see:  (102k)
> > http://www.linz.govt.nz/rcs/linz/17880/difference_wgs84_nzgd2000.pdf
> >
> >
> >
> > As for the new projection, it is New Zealand Transverse Mercator (NZTM),
> > which is in terms of that new NZGD2k datum.
> > Details:
> > Datum: NZGD2k
> > Origin Latitude: 0° South
> > Origin Longitude: 173° East
> > False Northing: 10 000 000m N
> > False Easting: 1 600 000m E
> > Scale Factor: 0.9996
> >
> > Setting projection as tmerc and entering those terms goes smoothly.
> > Should NZTM get its own projection entry or is entering by hand with the
> > tmerc projection every time the way to go?
> 
> Yes that is the way to do it until GRASS has co-ordinate system support,
> which I don't see happening any time soon.. you have it very easy with
> New Zealand Map Grid which puts in all the parameters for you.
> 
> >
> > for more details see:  (44k)
> > http://www.linz.govt.nz/rcs/linz/5684/nztransverse_mercator.pdf
> >
> > That PDF also has a page of useful ellipsoid to grid formulae which
> > might be useful to someone.
> >
> >
> > It doesn't help much that half the country is moving 5cm/year in the
> > opposite direction to the other half..
> 
> It would be useful to keep an eye out for any Free Software that handles
> that level of accuracy
> http://www-gpsg.mit.edu/~simon/gtgk/ has some interesting links
> 
> >
> >
> > thanks for any insight,
> 
> Well hopefully someone will correct the inevitable errors and
> misunderstandings in what I've written above
> 
> Paul
> ----------------------------------------
> PROJ.4 Discussion List
> See http://www.remotesensing.org/proj for subscription, unsubscription
> and other information.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.maptools.org/pipermail/proj/attachments/20030723/ccbe7ebe/attachment.html