root/Chameleon/trunk/Chameleon/geomath.php

<?php
/**
 * Math utils function
 *
 * @project     CWC2
 * @revision    $Id: 
 * @purpose     FindDealers Widget class
 * @author      DM Solutions Group (lacroix@dmsolutions.ca)
 * @copyright
 * <b>Copyright (c) 2003, DM Solutions Group Inc.</b>
 * Permission is hereby granted, free of charge, to any person obtaining a
 * copy of this software and associated documentation files (the "Software"),
 * to deal in the Software without restriction, including without limitation
 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
 * and/or sell copies of the Software, and to permit persons to whom the
 * Software is furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included
 * in all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL
 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 * DEALINGS IN THE SOFTWARE.
 */

/**
 * Number of KM in a MILE
 * there are 5280 feet in a mile        : start with 5280
 * there are 12 inches in a foot        : multiple by 12
 * there are 0.0254 meters in an inch   : multiple by 0.0254
 * there are 1000 meters in a kilometer : divide by 1000
 */
define( "KM_PER_MILE", (5280 * 12 * 0.0254 / 1000) );

/**
 * Constant returns miles per degree of latitude which is 69.
 */
define( "MilesPerDegreeLat", 69 );
define( "DegreesLatPerMile", (1 / 69.0) );


/**
 * function LonToLatRatio
 *
 * Given a latitude in decimal degrees, it will return the ratio of
 * longitude degress to latitude degrees.
 */
function LonToLatRatio ($latitude) 
{
    return cos($latitude * M_PI / 180);
}

/**
 * MilesPerDegreeLon($latitude)
 *
 * Given a latitude in decimal degrees, it will return the
 * miles per degree of longitude.
 */
function MilesPerDegreeLon ($latitude)
{
    return 69 * LonToLatRatio($latitude);
}

/**
 * DegreesLonPerMile($latitude)
 *
 * Given a latitude in decimal degrees, it will return the
 * degrees of longitude per mile.
 */
function DegreesLonPerMile ($latitude)
{
    return 1 / MilesPerDegreeLon( $latitude );
}

/**
 * decimal_degrees($degrees, $minutes, $seconds)
 *
 * Given whole degrees, minutes, and seconds, B<decimal_degrees> will
 * return the equivalent number of decimal degrees. B<dd> function
 * takes the same inputs as B<decimal_degrees> but returns a result
 * limited to 5 decimal places.
 */
function decimal_degrees ($degrees, $minutes, $seconds)
{
    // if degrees are negative, subtrack minutes and seconds.
    if ($degrees < 0) {
        return $degrees - ( $minutes / 60 ) - ( $seconds / 3600);
    }
    // if degrees are positve, add minutes and seconds.
    return $degrees + ( $minutes / 60 ) + ( $seconds / 3600);
}

function dd ($degrees, $minutes, $seconds)
{
    return sprintf("%.5f", decimal_degrees($degrees, $minutes, $seconds));
}

/**
 * Lon2DMS(I<longitude>)
 *
 * Given a longitude in decimal degrees, returns a string in the form of:
 *
 * I<degrees>B<d>I<minutes>B<'>I<seconds>B<">I<cardinal>
 *
 * where I<cardinal> is C<W> for negative longitudes and C<E> for
 * positive longitudes.
 */
function Lon2DMS ($longitude)
{
    $degrees = intval( $longitude );
    $minutes = intval( ( $longitude - $degrees ) * 60 );
    $seconds = intval( ( $longitude - $degrees ) * 3600 ) - ( $minutes * 60 );

      // west if longitude is negative
    if ($longitude < 0) {
        return abs($degrees) . 'd' .
               abs($minutes) . "'" .
               sprintf('%.3f', abs($seconds)) . '"W';
    }
      // east if longitude is positive
    return abs($degrees) . 'd' .
           abs($minutes) . "'" .
           sprintf('%.3f', abs($seconds)) . '"E';
}

/**
 * Lat2DMS(I<latitude>)
 *
 * Given a latitude in decimal degrees, returns a string in the form of:
 *
 * I<degrees>B<d>I<minutes>B<'>I<seconds>B<">I<cardinal>
 *
 * where I<cardinal> is C<S> for negative longitudes and C<N> for
 * positive longitudes.
 */
function Lat2DMS ($latitude)
{
    $degrees = intval( $latitude );
    $minutes = intval( ( $latitude - $degrees ) * 60 );
    $seconds = intval( ( $latitude - $degrees ) * 3600 ) - ( $minutes * 60 );

    // south if latitude is negative
    if ($latitude < 0) {
        return abs($degrees) . 'd' .
            abs($minutes) . "'" .
            sprintf('%.3f', abs($seconds)) . '"S';
    }
    // north if latitude is positive
    return abs($degrees) . 'd' .
        abs($minutes) . "'" .
        sprintf('%.3f', abs($seconds)) . '"N';
}

/*
 * deprecated - only used by mapimagesharedresource widget, changed to Pix2Geo (below)
function Pix2Georef ( $nX,   $nY,
             $dfGeoMiddleX, $dfGeoMiddleY,
             $dfGeoWidth,   $dfGeoHeight,
             $nPixelWidth,  $nPixelHeight  )
{

    // --------------------------------------------------------------------
    //      Calculate first the extent in georef as does mapserver
    //      (fcution msAdjustScale in maputil.c).
    // --------------------------------------------------------------------
    
    $dfGeoMinX = $dfGeoMiddleX - ($dfGeoWidth/2);
    $dfGeoMaxX = $dfGeoMiddleX + ($dfGeoWidth/2);
    $dfGeoMinY = $dfGeoMiddleY - ($dfGeoHeight/2);
    $dfGeoMaxY = $dfGeoMiddleY + ($dfGeoHeight/2);
    
    $dfTmp1 = ($dfGeoMaxX - $dfGeoMinX)/$nPixelWidth;
    $dfTmp2 = ($dfGeoMaxY - $dfGeoMinY)/$nPixelHeight;
    
    $cellsize = max($dfTmp1, $dfTmp2);
    
    if ($cellsize > 0) 
    {
    
        $dfTmp1 = max(($nPixelWidth - ($dfGeoMaxX - $dfGeoMinX)/$cellsize)/2 , 0);
        $dfTmp2 = max(($nPixelHeight - ($dfGeoMaxY - $dfGeoMinY)/$cellsize)/2 , 0);
    
        $ox = $oy = 0;
    
        if ($dfTmp1 >= 0.0) 
        {
            $ox =   intval($dfTmp1 + 0.5);
        }
        else 
        {
            $ox =   intval($dfTmp1 - 0.5);
        }
    
        if ($dfTmp2 >= 0.0) {
            $oy =   intval($dfTmp2 + 0.5);
        }
        else 
        {
            $oy =   intval($dfTmp2 - 0.5);
        }
    
    
        $dfGeoMinX = $dfGeoMinX - $ox*$cellsize;
        $dfGeoMaxX = $dfGeoMaxX + $ox*$cellsize;
        $dfGeoMinY = $dfGeoMinY - $oy*$cellsize;
        $dfGeoMaxY = $dfGeoMaxY + $oy*$cellsize;
    
    
        $dfPix2GeoX = ($dfGeoMaxX - $dfGeoMinX) / $nPixelWidth;
        $dfPix2GeoY = ($dfGeoMaxY - $dfGeoMinY) / $nPixelHeight;
    
        $nDeltaPixX =  $nX;
        $nDeltaPixY =  $nPixelHeight - $nY;# - $nPixelHeight;
    
    
        $dfDeltaGeoX = $nDeltaPixX * $dfPix2GeoX;
        $dfDeltaGeoY = $nDeltaPixY * $dfPix2GeoY;
    
        $dfPosGeoX = $dfGeoMinX + $dfDeltaGeoX;
        $dfPosGeoY = $dfGeoMinY + $dfDeltaGeoY;
    
        $aLatLong = array();
        $aLatLong[0] = $dfPosGeoX;
        $aLatLong[1] = $dfPosGeoY;
    
        return $aLatLong;
    }
}
*/
    
/**
 * convert a pixel position to a georef position
 *
 * @param nPixPos integer the pixel position
 * @param dfPixMin float the minimum pixel position
 * @param dfPixMax float the maximum pixel position
 * @param dfGeoMin float the minimum georeferenced position
 * @param dfGeoMax float the maximum georeferenced position
 * @param nInversePix boolean, flag to invert Y coordinates if UL > LR
 * @return double geocoded position.
 */
function Pix2Geo($nPixPos, $dfPixMin, $dfPixMax, $dfGeoMin, $dfGeoMax,
                 $nInversePix = false)
{
    $dfWidthGeo = $dfGeoMax - $dfGeoMin;
    $dfWidthPix = $dfPixMax - $dfPixMin;

    $dfPixToGeo = $dfWidthGeo / $dfWidthPix;

    if (!$nInversePix)
        $dfDeltaPix = $nPixPos - $dfPixMin;
    else
        $dfDeltaPix = $dfPixMax - $nPixPos;

    $dfDeltaGeo = $dfDeltaPix * $dfPixToGeo;


    $dfPosGeo = $dfGeoMin + $dfDeltaGeo;

    return ($dfPosGeo);
}
  
/**
 * convert a geocoded position to pixel coord
 *
 * @param nGeoPos double Geocoded position
 * @param dfPixMin double minimum map pixel value
 * @param dfPixMax double maximum map pixel value
 * @param dfGeoMin double minimum map geocoded value
 * @param dfGeoMax double maximum map geocoded value
 * @param nInverseGeo integer optional flag to inverse , set to 1 for
 *                            Y pixel coordinates where UL > LR
 * @return double geocoded position
 */
function Geo2Pix ($nGeoPos, $dfPixMin, $dfPixMax, $dfGeoMin,
                     $dfGeoMax, $nInverseGeo = "")
{
    // calculate the geocoded & pixel width
    $dfWidthGeo = abs($dfGeoMax - $dfGeoMin);
    $dfWidthPix = abs($dfPixMax - $dfPixMin);

    // get ratio
    if ( $dfWidthGeo <= 0 )
    {
        return 0;
    }
    $dfGeoToPix = $dfWidthPix / $dfWidthGeo;

    // get difference
    if (!$nInverseGeo)
        $dfDeltaGeo = $nGeoPos - $dfGeoMin;
    else
        $dfDeltaGeo = $dfGeoMax - $nGeoPos;

    // calculate
    $dfDeltaPix = $dfDeltaGeo * $dfGeoToPix;
    $dfPosPix = $dfPixMin + $dfDeltaPix;

    // return value
    return round ($dfPosPix);

    // end pixel_to_geo function
}
    
/**
 * reproject a point nX,nY
 * no reprojection happens if either projection is
 * empty or if they are the same. Will prepend
 * init= if either contains epsg (or EPSG)
 *
 * @param nX float reference the x coordinate to project
 * @param nY float reference to the y coordinate to project
 * @param szFrom string the projection of nX,nY
 * @param szTo string the projection to convert to
 * @return none.  Variables nX, nY are modified directly
 */
function reprojectPoint( &$nX, &$nY, $szFrom, $szTo )
{
    //echo "reproject $nX,$nY from $szFrom to $szTo<BR>\n";
    
    //check to see if reprojection is necessary
    if ($szFrom == '' || $szTo == '')
    {
        return;
    }
    
    if (stristr($szFrom, "epsg") !== false &&
        stristr( $szFrom, "init=") == false)
    {
        $szFrom = "init=".strtolower($szFrom);
    }
    if (stristr($szTo, "epsg") !== false &&
        stristr( $szTo, "init=") == false)
    {
        $szTo = "init=".strtolower($szTo);
    }
    
    if ($szFrom == $szTo)
    {
        return;
    }
    
    $oPoint = ms_newPointObj();
    $oPoint->setXY( $nX, $nY );
    $oPoint->project( ms_newProjectionObj( $szFrom ),
                      ms_newProjectionObj( $szTo )
                       );
    $nX = $oPoint->x;
    $nY = $oPoint->y;    
    
    //echo "reprojection result: $nX, $nY<BR>\n";
}
?>