{-# LANGUAGE CPP #-}moduleData.GPS.Core(-- * TypesDistance,Heading,Speed,Vector,Trail-- * Constants,north,south,east,west,radiusOfEarth,circumferenceOfEarth-- * Coordinate Functions,heading,distance,speed,getVector,addVector,getRadianPair,getDMSPair,divideArea,interpolate-- * IO helpers,writeGPX,readGPX,readGPXSegments-- * Utility,getUTCTime,moduleData.Geo.GPX)whereimportData.TimeimportData.MaybeimportControl.MonadimportText.XML.HXT.CoreimportText.XML.XSD.DateTime(DateTime,toUTCTime)importData.Geo.GPX-- |Distances are expressed in meterstypeDistance=Double-- |Angles are expressed in radians from North.-- 0 == North-- pi/2 == West-- pi == South-- (3/2)pi == East == - (pi / 2)typeHeading=DoubletypeAngle=Double-- |Speed is hard coded as meters per secondtypeSpeed=DoubletypeVector=(Distance,Heading)typeTraila=[a]getUTCTime::(Timea)=>a->MaybeUTCTimegetUTCTime=fmaptoUTCTime.timedistance::(Lata,Lona,Latb,Lonb)=>a->b->Distancedistancexy=let(lat1,lon1)=getRadianPairDx(lat2,lon2)=getRadianPairDydeltaLat=lat2-lat1deltaLon=lon2-lon1a=(sin(deltaLat/2))^2+coslat1*coslat2*(sin(deltaLon/2))^2c=2*atan2(a**0.5)((1-a)**0.5)inradiusOfEarth*c-- | Direction two points aim toward (0 = North, pi/2 = West, pi = South, 3pi/2 = East)heading::(Lata,Lona,Latb,Lonb)=>a->b->Headingheadingab=atan2(sin(diffLon)*cos(lat2))(cos(lat1)*sin(lat2)-sin(lat1)*coslat2*cos(diffLon))where(lat1,lon1)=getRadianPairDa(lat2,lon2)=getRadianPairDbdiffLon=lon2-lon1getVector::(Lata,Lona,Latb,Lonb)=>a->b->VectorgetVectorab=(distanceab,headingab)-- |Given a vector and coordinate, computes a new coordinate.-- Within some epsilon it should hold that if---- @dest = addVector (dist,heading) start@---- then---- @heading == heading start dest@-- -- @dist == distance start dest@addVector::(Latc,Lonc)=>Vector->c->caddVector(d,h)p=setLon(longitudeType$toDegreeslon2).setLat(latitudeType$toDegreeslat2)$pwhere(lat,lon)=getRadianPairDplat2=asin(sin(lat)*cos(d/radiusOfEarth)+cos(lat)*sin(d/radiusOfEarth)*cosh)lon2=lon+atan2(sinh*sin(d/radiusOfEarth)*coslat)(cos(d/radiusOfEarth)-sinlat*sinlat2)-- | Speed in meters per second, only if a 'Time' was recorded for each waypoint.speed::(Latloc,Lonloc,Timeloc,Latb,Lonb,Timeb)=>loc->b->MaybeSpeedspeedab=case(getUTCTimeb,getUTCTimea)of(Justx,Justy)->lettimeDiff=realToFrac(diffUTCTimexy)iniftimeDiff==0thenNothingelseJust$(distanceab)/timeDiff_->Nothing-- |radius of the earth in metersradiusOfEarth::DoubleradiusOfEarth=6378700-- |Circumference of earht (meters)circumferenceOfEarth::DoublecircumferenceOfEarth=radiusOfEarth*2*pi-- |North is 0 radiansnorth::Headingnorth=0-- |South, being 180 degrees from North, is pi.south::Headingsouth=pi-- |East is 270 degrees (3 pi / 2)east::Headingeast=(3/2)*pi-- |West is 90 degrees (pi/2)west::Headingwest=pi/2toDegrees=(*)(180/pi)getRadianPairD::(Latc,Lonc)=>c->(Double,Double)getRadianPairD=(\(a,b)->(realToFraca,realToFracb)).getRadianPairgetDMSPair::(Latc,Lonc)=>c->(LatitudeType,LongitudeType)getDMSPairc=(latc,lonc)-- |Provides a lat/lon pair of doubles in radiansgetRadianPair::(Latp,Lonp)=>p->(LatitudeType,LongitudeType)getRadianPairp=(toRadians(latp),toRadians(lonp))toRadians::Floatingf=>f->ftoRadians=(*)(pi/180)-- | @interpolate c1 c2 w@ where 0 <= w <= 1 Gives a point on the line-- between c1 and c2 equal to @c1 when @w == 0@ (weighted linearly-- toward c2).interpolate::(Lata,Lona)=>a->a->Double->ainterpolatec1c2w|w<0||w>1=error"Interpolate only works with a weight between zero and one"|otherwise=let(h,d)=(headingc1c2,distancec1c2)v=(d*w,h)inaddVectorvc1-- |@divideArea vDist hDist nw se@ divides an area into a grid of equally-- spaced coordinates within the box drawn by the northwest point (nw) and-- southeast point (se). Because this uses floating point there might be a-- different number of points in some rows (the last might be too far east based-- on a heading from the se point).divideArea::(Latc,Lonc)=>Distance->Distance->c->c->[[c]]divideAreavDisthDistnwse=let(top,left)=(latnw,lonnw)(btm,right)=(latse,lonse)columnOne=takeWhile((<=west).headingse).iterate(addVector(vDist,south))$nwbuildRow=takeWhile((>=north).headingse).iterate(addVector(hDist,east))inmapbuildRowcolumnOne-- |Reads a GPX file (using the GPX library) by simply concatenating all the-- tracks, segments, and points ('trkpts', 'trksegs', 'trks') into a single 'Trail'.readGPX::FilePath->IO(TrailWptType)readGPX=liftM(concatMaptrkpts.concatMaptrksegs.concatMaptrks).readGpxFilewriteGPX::FilePath->TrailWptType->IO()writeGPXfpps=writeGpxFilefp$gpx$gpxType"1.0""Haskell GPS Package (via the GPX package)"Nothing[][][trkTypeNothingNothingNothingNothing[]NothingNothingNothing[trksegTypepsNothing]]Nothing-- writeGpxFile should go in the GPX packagewriteGpxFile::FilePath->Gpx->IO()writeGpxFilefpgpx=runX_(constAgpx>>>xpickleDocument(xpickle::PUGpx)[]fp)runX_t=runXt>>return()readGPXSegments::FilePath->IO[TrailWptType]readGPXSegments=liftM(map(concatMaptrkpts).maptrksegs.concatMaptrks).readGpxFile