geo_extra.c

GIS系统支持库Geospatial Data Abstraction Library代码.GDAL is a translator library for raster geospatial dat

    PCS_NAD27_Texas_North,		Proj_Texas_CS27_North,
    PCS_NAD27_Texas_North_Cen,		Proj_Texas_CS27_North_Central,
    PCS_NAD27_Texas_Central,		Proj_Texas_CS27_Central,
    PCS_NAD27_Texas_South_Cen,		Proj_Texas_CS27_South_Central,
    PCS_NAD27_Texas_South,		Proj_Texas_CS27_South,

    PCS_NAD27_Utah_North,		Proj_Utah_CS27_North,
    PCS_NAD27_Utah_Central,		Proj_Utah_CS27_Central,
    PCS_NAD27_Utah_South,		Proj_Utah_CS27_South,

    PCS_NAD27_Vermont,			Proj_Vermont_CS27,

    PCS_NAD27_Virginia_North,		Proj_Virginia_CS27_North,
    PCS_NAD27_Virginia_South,		Proj_Virginia_CS27_South,

    PCS_NAD27_Washington_North,		Proj_Washington_CS27_North,
    PCS_NAD27_Washington_South,		Proj_Washington_CS27_South,

    PCS_NAD27_West_Virginia_N,		Proj_West_Virginia_CS27_North,
    PCS_NAD27_West_Virginia_S,		Proj_West_Virginia_CS27_South,

    PCS_NAD27_Wisconsin_North,		Proj_Wisconsin_CS27_North,
    PCS_NAD27_Wisconsin_Cen,		Proj_Wisconsin_CS27_Central,
    PCS_NAD27_Wisconsin_South,		Proj_Wisconsin_CS27_South,

    PCS_NAD27_Wyoming_East,		Proj_Wyoming_CS27_East,
    PCS_NAD27_Wyoming_E_Cen,		Proj_Wyoming_CS27_East_Central,
    PCS_NAD27_Wyoming_W_Cen,		Proj_Wyoming_CS27_West_Central,
    PCS_NAD27_Wyoming_West,		Proj_Wyoming_CS27_West,
    
    PCS_NAD27_Puerto_Rico,		Proj_Puerto_Rico_CS27,

    KvUserDefined
};

/************************************************************************/
/*                          GTIFMapSysToPCS()                           */
/*                                                                      */
/*      Given a Datum, MapSys and zone value generate the best PCS      */
/*      code possible.                                                  */
/************************************************************************/

int	GTIFMapSysToPCS( int MapSys, int Datum, int nZone )

{
    int		PCSCode = KvUserDefined;

    if( MapSys == MapSys_UTM_North )
    {
	if( Datum == GCS_NAD27 )
	    PCSCode = PCS_NAD27_UTM_zone_3N + nZone - 3;
	else if( Datum == GCS_NAD83 )
	    PCSCode = PCS_NAD83_UTM_zone_3N + nZone - 3;
	else if( Datum == GCS_WGS_72 )
	    PCSCode = PCS_WGS72_UTM_zone_1N + nZone - 1;
	else if( Datum == GCS_WGS_72BE )
	    PCSCode = PCS_WGS72BE_UTM_zone_1N + nZone - 1;
	else if( Datum == GCS_WGS_84 )
	    PCSCode = PCS_WGS84_UTM_zone_1N + nZone - 1;
    }
    else if( MapSys == MapSys_UTM_South )
    {
	if( Datum == GCS_WGS_72 )
	    PCSCode = PCS_WGS72_UTM_zone_1S + nZone - 1;
	else if( Datum == GCS_WGS_72BE )
	    PCSCode = PCS_WGS72BE_UTM_zone_1S + nZone - 1;
	else if( Datum == GCS_WGS_84 )
	    PCSCode = PCS_WGS84_UTM_zone_1S + nZone - 1;
    }
    else if( MapSys == MapSys_State_Plane_27 )
    {
	int		i;

        PCSCode = 10000 + nZone;
	for( i = 0; StatePlaneTable[i] != KvUserDefined; i += 2 )
	{
	    if( StatePlaneTable[i+1] == PCSCode )
	        PCSCode = StatePlaneTable[i];
	}

        /* Old EPSG code was in error for Tennesse CS27, override */
        if( nZone == 4100 )
            PCSCode = 2204;
    }
    else if( MapSys == MapSys_State_Plane_83 )
    {
	int		i;

        PCSCode = 10000 + nZone + 30;

	for( i = 0; StatePlaneTable[i] != KvUserDefined; i += 2 )
	{
	    if( StatePlaneTable[i+1] == PCSCode )
	        PCSCode = StatePlaneTable[i];
	}

        /* Old EPSG code was in error for Kentucky North CS83, override */
        if( nZone == 1601 )
            PCSCode = 2205;
    }

    return( PCSCode );
}

/************************************************************************/
/*                          GTIFMapSysToProj()                          */
/*                                                                      */
/*      Given a MapSys and zone value generate the best Proj_           */
/*      code possible.                                                  */
/************************************************************************/

int	GTIFMapSysToProj( int MapSys, int nZone )

{
    int		ProjCode = KvUserDefined;

    if( MapSys == MapSys_UTM_North )
    {
        ProjCode = Proj_UTM_zone_1N + nZone - 1;
    }
    else if( MapSys == MapSys_UTM_South )
    {
        ProjCode = Proj_UTM_zone_1S + nZone - 1;
    }
    else if( MapSys == MapSys_State_Plane_27 )
    {
        ProjCode = 10000 + nZone;

        /* Tennesse override */
        if( nZone == 4100 )
            ProjCode = 15302;
    }
    else if( MapSys == MapSys_State_Plane_83 )
    {
        ProjCode = 10000 + nZone + 30;

        /* Kentucky North override */
        if( nZone == 1601 )
            ProjCode = 15303;
    }

    return( ProjCode );
}

/************************************************************************/
/*                          GTIFPCSToMapSys()                           */
/************************************************************************/

/**
 * Translate a PCS_ code into a UTM or State Plane map system, a datum,
 * and a zone if possible.
 *
 * @param PCSCode The projection code (PCS_*) as would be stored in the
 * ProjectedCSTypeGeoKey of a GeoTIFF file.
 *
 * @param pDatum Pointer to an integer into which the datum code (GCS_*)
 * is put if the function succeeds.
 *
 * @param pZone Pointer to an integer into which the zone will be placed
 * if the function is successful.
 *
 * @return Returns either MapSys_UTM_North, MapSys_UTM_South,
 * MapSys_State_Plane_83, MapSys_State_Plane_27 or KvUserDefined.
 * KvUserDefined indicates that the
 * function failed to recognise the projection as UTM or State Plane.
 *
 * The zone value is only set if the return code is other than KvUserDefined.
 * For utm map system the returned zone will be between 1 and 60.  For
 * State Plane, the USGS state plane zone number is returned.  For instance,
 * Alabama East is zone 101.
 *
 * The datum (really this is the GCS) is set to a GCS_ value such as GCS_NAD27.
 *
 * This function is useful to recognise (most) UTM and State Plane coordinate
 * systems, even if CSV files aren't available to translate them automatically.
 * It is used as a fallback mechanism by GTIFGetDefn() for normalization when
 * CSV files aren't found. 
 */

int GTIFPCSToMapSys( int PCSCode, int * pDatum, int * pZone )

{
    int		Datum = KvUserDefined, Proj = KvUserDefined;
    int		nZone = KvUserDefined, i;

/* -------------------------------------------------------------------- */
/*      UTM with various datums.  Note there are lots of PCS UTM        */
/*      codes not done yet which use strange datums.                    */
/* -------------------------------------------------------------------- */
    if( PCSCode >= PCS_NAD27_UTM_zone_3N && PCSCode <= PCS_NAD27_UTM_zone_22N )
    {
	Datum = GCS_NAD27;
	Proj = MapSys_UTM_North;
	nZone = PCSCode - PCS_NAD27_UTM_zone_3N + 3;
    }
    else if( PCSCode >= PCS_NAD83_UTM_zone_3N 
	     && PCSCode <= PCS_NAD83_UTM_zone_23N )
    {
	Datum = GCS_NAD83;
	Proj = MapSys_UTM_North;
	nZone = PCSCode - PCS_NAD83_UTM_zone_3N + 3;
    }

    else if( PCSCode >= PCS_WGS72_UTM_zone_1N
	     && PCSCode <= PCS_WGS72_UTM_zone_60N )
    {
	Datum = GCS_WGS_72;
	Proj = MapSys_UTM_North;
	nZone = PCSCode - PCS_WGS72_UTM_zone_1N + 1;
    }
    else if( PCSCode >= PCS_WGS72_UTM_zone_1S
	     && PCSCode <= PCS_WGS72_UTM_zone_60S )
    {
	Datum = GCS_WGS_72;
	Proj = MapSys_UTM_South;
	nZone = PCSCode - PCS_WGS72_UTM_zone_1S + 1;
    }

    else if( PCSCode >= PCS_WGS72BE_UTM_zone_1N
	     && PCSCode <= PCS_WGS72BE_UTM_zone_60N )
    {
	Datum = GCS_WGS_72BE;
	Proj = MapSys_UTM_North;
	nZone = PCSCode - PCS_WGS72BE_UTM_zone_1N + 1;
    }
    else if( PCSCode >= PCS_WGS72BE_UTM_zone_1S
	     && PCSCode <= PCS_WGS72BE_UTM_zone_60S )
    {
	Datum = GCS_WGS_72BE;
	Proj = MapSys_UTM_South;
	nZone = PCSCode - PCS_WGS72BE_UTM_zone_1S + 1;
    }

    else if( PCSCode >= PCS_WGS84_UTM_zone_1N
	     && PCSCode <= PCS_WGS84_UTM_zone_60N )
    {
	Datum = GCS_WGS_84;
	Proj = MapSys_UTM_North;
	nZone = PCSCode - PCS_WGS84_UTM_zone_1N + 1;
    }
    else if( PCSCode >= PCS_WGS84_UTM_zone_1S
	     && PCSCode <= PCS_WGS84_UTM_zone_60S )
    {
	Datum = GCS_WGS_84;
	Proj = MapSys_UTM_South;
	nZone = PCSCode - PCS_WGS84_UTM_zone_1S + 1;
    }
    else if( PCSCode >= PCS_SAD69_UTM_zone_18N 
	     && PCSCode <= PCS_SAD69_UTM_zone_22N )
    {
	Datum = KvUserDefined;
	Proj = MapSys_UTM_North;
	nZone = PCSCode - PCS_SAD69_UTM_zone_18N + 18;
    }
    else if( PCSCode >= PCS_SAD69_UTM_zone_17S
	     && PCSCode <= PCS_SAD69_UTM_zone_25S )
    {
	Datum = KvUserDefined;
	Proj = MapSys_UTM_South;
	nZone = PCSCode - PCS_SAD69_UTM_zone_17S + 17;
    }

/* -------------------------------------------------------------------- */
/*      State Plane zones, first we translate any PCS_ codes to		*/
/*	a Proj_ code that we can get a handle on.			*/
/* -------------------------------------------------------------------- */
    for( i = 0; StatePlaneTable[i] != KvUserDefined; i += 2 )
    {
	if( StatePlaneTable[i] == PCSCode )
	    PCSCode = StatePlaneTable[i+1];
    }

    if( PCSCode <= 15900 && PCSCode >= 10000 )
    {
	if( (PCSCode % 100) >= 30 )
        {
            Proj = MapSys_State_Plane_83;
	    Datum = GCS_NAD83;
        }
	else
        {
            Proj = MapSys_State_Plane_27;
	    Datum = GCS_NAD27;
        }
	
	nZone = PCSCode - 10000;
	if( Datum == GCS_NAD83 )
	    nZone -= 30;
    }

    if( pDatum != NULL )
        *pDatum = Datum;

    if( pZone != NULL )
        *pZone = nZone;

    return( Proj );
}

/************************************************************************/
/*                          GTIFProjToMapSys()                          */
/************************************************************************/

/**
 * Translate a Proj_ code into a UTM or State Plane map system, and a zone
 * if possible.
 *
 * @param ProjCode The projection code (Proj_*) as would be stored in the
 * ProjectionGeoKey of a GeoTIFF file.
 * @param pZone Pointer to an integer into which the zone will be placed
 * if the function is successful.
 *
 * @return Returns either MapSys_UTM_North, MapSys_UTM_South,
 * MapSys_State_Plane_27, MapSys_State_Plane_83 or KvUserDefined.
 * KvUserDefined indicates that the
 * function failed to recognise the projection as UTM or State Plane.
 *
 * The zone value is only set if the return code is other than KvUserDefined.
 * For utm map system the returned zone will be between 1 and 60.  For
 * State Plane, the USGS state plane zone number is returned.  For instance,
 * Alabama East is zone 101.
 *
 * This function is useful to recognise UTM and State Plane coordinate
 * systems, and to extract zone numbers so the projections can be
 * represented as UTM rather than as the underlying projection method such
 * Transverse Mercator for instance.
 */

int GTIFProjToMapSys( int ProjCode, int * pZone )

{
    int		nZone = KvUserDefined;
    int		MapSys = KvUserDefined;

/* -------------------------------------------------------------------- */
/*      Handle UTM.                                                     */
/* -------------------------------------------------------------------- */
    if( ProjCode >= Proj_UTM_zone_1N && ProjCode <= Proj_UTM_zone_60N )
    {
	MapSys = MapSys_UTM_North;
	nZone = ProjCode - Proj_UTM_zone_1N + 1;
    }
    else if( ProjCode >= Proj_UTM_zone_1S && ProjCode <= Proj_UTM_zone_60S )
    {
	MapSys = MapSys_UTM_South;
	nZone = ProjCode - Proj_UTM_zone_1S + 1;
    }

/* -------------------------------------------------------------------- */
/*      Handle State Plane.  I think there are some anomolies in        */
/*      here, so this is a bit risky.                                   */
/* -------------------------------------------------------------------- */
    else if( ProjCode >= 10101 && ProjCode <= 15299 )
    {
        if( ProjCode % 100 >= 30 )
        {
            MapSys = MapSys_State_Plane_83;
            nZone = ProjCode - 10000 - 30;
        }
        else
        {
            MapSys = MapSys_State_Plane_27;
            nZone = ProjCode - 10000;
        }
    }
    
    if( pZone != NULL )
        *pZone = nZone;

    return( MapSys );
}