[Liblas-devel] Re: GetProj4() string and NAD83
Martin Isenburg
martin.isenburg at gmail.com
Wed Jan 11 12:45:30 EST 2012
hi adam,
should be an easy fix in that you can do in your GDAL code look how
the projection setting is done in GeoProjectionConverter.cpp from
LASlib. There it only sets up a user defined projection if all
parameters for that projection were provided. This happens in the
set_projection_from_geo_keys() function (see below).
cheers,
martin
bool GeoProjectionConverter::set_projection_from_geo_keys(int
num_geo_keys, GeoProjectionGeoKeys* geo_keys, char* geo_ascii_params,
double* geo_double_params, char* description)
{
bool user_defined_ellipsoid = false;
int user_defined_projection = 0;
int offsetProjStdParallel1GeoKey = -1;
int offsetProjStdParallel2GeoKey = -1;
int offsetProjNatOriginLatGeoKey = -1;
int offsetProjFalseEastingGeoKey = -1;
int offsetProjFalseNorthingGeoKey = -1;
int offsetProjCenterLongGeoKey = -1;
int offsetProjScaleAtNatOriginGeoKey = -1;
bool has_projection = false;
int ellipsoid = -1;
this->num_geo_keys = num_geo_keys;
this->geo_keys = geo_keys;
this->geo_ascii_params = geo_ascii_params;
this->geo_double_params = geo_double_params;
for (int i = 0; i < num_geo_keys; i++)
{
switch (geo_keys[i].key_id)
{
case 1024: // GTModelTypeGeoKey
if (geo_keys[i].value_offset == 2) // ModelTypeGeographic
{
has_projection = set_longlat_projection(description);
}
break;
case 2048: // GeographicTypeGeoKey
switch (geo_keys[i].value_offset)
{
case 32767: // user-defined GCS
user_defined_ellipsoid = true;
break;
case 4001: // GCSE_Airy1830
ellipsoid = 1;
break;
case 4002: // GCSE_AiryModified1849
ellipsoid = 16;
break;
case 4003: // GCSE_AustralianNationalSpheroid
ellipsoid = 2;
break;
case 4004: // GCSE_Bessel1841
case 4005: // GCSE_Bessel1841Modified
ellipsoid = 3;
break;
case 4006: // GCSE_BesselNamibia
ellipsoid = 4;
break;
case 4008: // GCSE_Clarke1866
case 4009: // GCSE_Clarke1866Michigan
ellipsoid = GEO_ELLIPSOID_NAD27;
break;
case 4010: // GCSE_Clarke1880_Benoit
case 4011: // GCSE_Clarke1880_IGN
case 4012: // GCSE_Clarke1880_RGS
case 4013: // GCSE_Clarke1880_Arc
case 4014: // GCSE_Clarke1880_SGA1922
case 4034: // GCSE_Clarke1880
ellipsoid = 6;
break;
case 4015: // GCSE_Everest1830_1937Adjustment
case 4016: // GCSE_Everest1830_1967Definition
case 4017: // GCSE_Everest1830_1975Definition
ellipsoid = 7;
break;
case 4018: // GCSE_Everest1830Modified
ellipsoid = 17;
break;
case 4019: // GCSE_GRS1980
ellipsoid = GEO_ELLIPSOID_NAD83;
break;
case 4020: // GCSE_Helmert1906
ellipsoid = 12;
break;
case 4022: // GCSE_International1924
case 4023: // GCSE_International1967
ellipsoid = 14;
break;
case 4024: // GCSE_Krassowsky1940
ellipsoid = 15;
break;
case 4030: // GCSE_WGS84
ellipsoid = GEO_ELLIPSOID_WGS84;
break;
case 4267: // GCS_NAD27
ellipsoid = GEO_ELLIPSOID_NAD27;
break;
case 4269: // GCS_NAD83
ellipsoid = GEO_ELLIPSOID_NAD83;
break;
case 4322: // GCS_WGS_72
ellipsoid = GEO_ELLIPSOID_WGS72;
break;
case 4326: // GCS_WGS_84
ellipsoid = GEO_ELLIPSOID_WGS84;
break;
default:
fprintf(stderr, "GeographicTypeGeoKey: look-up for %d not
implemented\n", geo_keys[i].value_offset);
}
break;
case 2050: // GeogGeodeticDatumGeoKey
switch (geo_keys[i].value_offset)
{
case 32767: // user-defined GCS
user_defined_ellipsoid = true;
break;
case 6202: // Datum_Australian_Geodetic_Datum_1966
case 6203: // Datum_Australian_Geodetic_Datum_1984
ellipsoid = 2;
break;
case 6267: // Datum_North_American_Datum_1927
ellipsoid = GEO_ELLIPSOID_NAD27;
break;
case 6269: // Datum_North_American_Datum_1983
ellipsoid = GEO_ELLIPSOID_NAD83;
break;
case 6322: // Datum_WGS72
ellipsoid = GEO_ELLIPSOID_WGS72;
break;
case 6326: // Datum_WGS84
ellipsoid = GEO_ELLIPSOID_WGS84;
break;
case 6001: // DatumE_Airy1830
ellipsoid = 1;
break;
case 6002: // DatumE_AiryModified1849
ellipsoid = 16;
break;
case 6003: // DatumE_AustralianNationalSpheroid
ellipsoid = 2;
break;
case 6004: // DatumE_Bessel1841
case 6005: // DatumE_BesselModified
ellipsoid = 3;
break;
case 6006: // DatumE_BesselNamibia
ellipsoid = 4;
break;
case 6008: // DatumE_Clarke1866
case 6009: // DatumE_Clarke1866Michigan
ellipsoid = GEO_ELLIPSOID_NAD27;
break;
case 6010: // DatumE_Clarke1880_Benoit
case 6011: // DatumE_Clarke1880_IGN
case 6012: // DatumE_Clarke1880_RGS
case 6013: // DatumE_Clarke1880_Arc
case 6014: // DatumE_Clarke1880_SGA1922
case 6034: // DatumE_Clarke1880
ellipsoid = 6;
break;
case 6015: // DatumE_Everest1830_1937Adjustment
case 6016: // DatumE_Everest1830_1967Definition
case 6017: // DatumE_Everest1830_1975Definition
ellipsoid = 7;
break;
case 6018: // DatumE_Everest1830Modified
ellipsoid = 17;
break;
case 6019: // DatumE_GRS1980
ellipsoid = GEO_ELLIPSOID_NAD83;
break;
case 6020: // DatumE_Helmert1906
ellipsoid = 12;
break;
case 6022: // DatumE_International1924
case 6023: // DatumE_International1967
ellipsoid = 14;
break;
case 6024: // DatumE_Krassowsky1940
ellipsoid = 15;
break;
case 6030: // DatumE_WGS84
ellipsoid = GEO_ELLIPSOID_WGS84;
break;
default:
fprintf(stderr, "GeogGeodeticDatumGeoKey: look-up for %d not
implemented\n", geo_keys[i].value_offset);
}
break;
case 2052: // GeogLinearUnitsGeoKey
switch (geo_keys[i].value_offset)
{
case 9001: // Linear_Meter
set_coordinates_in_meter();
break;
case 9002: // Linear_Foot
set_coordinates_in_feet();
break;
case 9003: // Linear_Foot_US_Survey
set_coordinates_in_survey_feet();
break;
default:
fprintf(stderr, "GeogLinearUnitsGeoKey: look-up for %d not
implemented\n", geo_keys[i].value_offset);
}
break;
case 2056: // GeogEllipsoidGeoKey
switch (geo_keys[i].value_offset)
{
case 7001: // Ellipse_Airy_1830
ellipsoid = 1;
break;
case 7002: // Ellipse_Airy_Modified_1849
ellipsoid = 16;
break;
case 7003: // Ellipse_Australian_National_Spheroid
ellipsoid = 2;
break;
case 7004: // Ellipse_Bessel_1841
case 7005: // Ellipse_Bessel_Modified
ellipsoid = 3;
break;
case 7006: // Ellipse_Bessel_Namibia
ellipsoid = 4;
break;
case 7008: // Ellipse_Clarke_1866
case 7009: // Ellipse_Clarke_1866_Michigan
ellipsoid = GEO_ELLIPSOID_NAD27;
break;
case 7010: // Ellipse_Clarke1880_Benoit
case 7011: // Ellipse_Clarke1880_IGN
case 7012: // Ellipse_Clarke1880_RGS
case 7013: // Ellipse_Clarke1880_Arc
case 7014: // Ellipse_Clarke1880_SGA1922
case 7034: // Ellipse_Clarke1880
ellipsoid = 6;
break;
case 7015: // Ellipse_Everest1830_1937Adjustment
case 7016: // Ellipse_Everest1830_1967Definition
case 7017: // Ellipse_Everest1830_1975Definition
ellipsoid = 7;
break;
case 7018: // Ellipse_Everest1830Modified
ellipsoid = 17;
break;
case 7019: // Ellipse_GRS_1980
ellipsoid = GEO_ELLIPSOID_NAD83;
break;
case 7020: // Ellipse_Helmert1906
ellipsoid = 12;
break;
case 7022: // Ellipse_International1924
case 7023: // Ellipse_International1967
ellipsoid = 14;
break;
case 7024: // Ellipse_Krassowsky1940
ellipsoid = 15;
break;
case 7030: // Ellipse_WGS_84
ellipsoid = GEO_ELLIPSOID_WGS84;
break;
default:
fprintf(stderr, "GeogEllipsoidGeoKey: look-up for %d not
implemented\n", geo_keys[i].value_offset);
}
break;
case 3072: // ProjectedCSTypeGeoKey
if (geo_keys[i].value_offset != 32767)
has_projection =
set_ProjectedCSTypeGeoKey(geo_keys[i].value_offset, description);
break;
case 3075: // ProjCoordTransGeoKey
user_defined_projection = 0;
switch (geo_keys[i].value_offset)
{
case 1: // CT_TransverseMercator
user_defined_projection = 1;
break;
case 8: // CT_LambertConfConic_2SP
user_defined_projection = 8;
break;
case 2: // CT_TransvMercator_Modified_Alaska
fprintf(stderr, "ProjCoordTransGeoKey:
CT_TransvMercator_Modified_Alaska not implemented\n");
break;
case 3: // CT_ObliqueMercator
fprintf(stderr, "ProjCoordTransGeoKey: CT_ObliqueMercator not
implemented\n");
break;
case 4: // CT_ObliqueMercator_Laborde
fprintf(stderr, "ProjCoordTransGeoKey:
CT_ObliqueMercator_Laborde not implemented\n");
break;
case 5: // CT_ObliqueMercator_Rosenmund
fprintf(stderr, "ProjCoordTransGeoKey:
CT_ObliqueMercator_Rosenmund not implemented\n");
break;
case 6: // CT_ObliqueMercator_Spherical
fprintf(stderr, "ProjCoordTransGeoKey:
CT_ObliqueMercator_Spherical not implemented\n");
break;
case 7: // CT_Mercator
fprintf(stderr, "ProjCoordTransGeoKey: CT_Mercator not implemented\n");
break;
case 9: // CT_LambertConfConic_Helmert
fprintf(stderr, "ProjCoordTransGeoKey:
CT_LambertConfConic_Helmert not implemented\n");
break;
case 10: // CT_LambertAzimEqualArea
fprintf(stderr, "ProjCoordTransGeoKey: CT_LambertAzimEqualArea
not implemented\n");
break;
case 11: // CT_AlbersEqualArea
fprintf(stderr, "ProjCoordTransGeoKey: CT_AlbersEqualArea not
implemented\n");
break;
case 12: // CT_AzimuthalEquidistant
fprintf(stderr, "ProjCoordTransGeoKey: CT_AzimuthalEquidistant
not implemented\n");
break;
case 13: // CT_EquidistantConic
fprintf(stderr, "ProjCoordTransGeoKey: CT_EquidistantConic not
implemented\n");
break;
case 14: // CT_Stereographic
fprintf(stderr, "ProjCoordTransGeoKey: CT_Stereographic not
implemented\n");
break;
case 15: // CT_PolarStereographic
fprintf(stderr, "ProjCoordTransGeoKey: CT_PolarStereographic
not implemented\n");
break;
case 16: // CT_ObliqueStereographic
fprintf(stderr, "ProjCoordTransGeoKey: CT_ObliqueStereographic
not implemented\n");
break;
case 17: // CT_Equirectangular
fprintf(stderr, "ProjCoordTransGeoKey: CT_Equirectangular not
implemented\n");
break;
case 18: // CT_CassiniSoldner
fprintf(stderr, "ProjCoordTransGeoKey: CT_CassiniSoldner not
implemented\n");
break;
case 19: // CT_Gnomonic
fprintf(stderr, "ProjCoordTransGeoKey: CT_Gnomonic not implemented\n");
break;
case 20: // CT_MillerCylindrical
fprintf(stderr, "ProjCoordTransGeoKey: CT_MillerCylindrical
not implemented\n");
break;
case 21: // CT_Orthographic
fprintf(stderr, "ProjCoordTransGeoKey: CT_Orthographic not
implemented\n");
break;
case 22: // CT_Polyconic
fprintf(stderr, "ProjCoordTransGeoKey: CT_Polyconic not implemented\n");
break;
case 23: // CT_Robinson
fprintf(stderr, "ProjCoordTransGeoKey: CT_Robinson not implemented\n");
break;
case 24: // CT_Sinusoidal
fprintf(stderr, "ProjCoordTransGeoKey: CT_Sinusoidal not
implemented\n");
break;
case 25: // CT_VanDerGrinten
fprintf(stderr, "ProjCoordTransGeoKey: CT_VanDerGrinten not
implemented\n");
break;
case 26: // CT_NewZealandMapGrid
fprintf(stderr, "ProjCoordTransGeoKey: CT_NewZealandMapGrid
not implemented\n");
break;
case 27: // CT_TransvMercator_SouthOriented
fprintf(stderr, "ProjCoordTransGeoKey:
CT_TransvMercator_SouthOriented not implemented\n");
break;
default:
fprintf(stderr, "ProjCoordTransGeoKey: look-up for %d not
implemented\n", geo_keys[i].value_offset);
}
break;
case 3076: // ProjLinearUnitsGeoKey
set_ProjLinearUnitsGeoKey(geo_keys[i].value_offset);
break;
case 3078: // ProjStdParallel1GeoKey
offsetProjStdParallel1GeoKey = geo_keys[i].value_offset;
break;
case 3079: // ProjStdParallel2GeoKey
offsetProjStdParallel2GeoKey = geo_keys[i].value_offset;
break;
case 3081: // ProjNatOriginLatGeoKey
offsetProjNatOriginLatGeoKey = geo_keys[i].value_offset;
break;
case 3082: // ProjFalseEastingGeoKey
offsetProjFalseEastingGeoKey = geo_keys[i].value_offset;
break;
case 3083: // ProjFalseNorthingGeoKey
offsetProjFalseNorthingGeoKey = geo_keys[i].value_offset;
break;
case 3088: // ProjCenterLongGeoKey
offsetProjCenterLongGeoKey = geo_keys[i].value_offset;
break;
case 3092: // ProjScaleAtNatOriginGeoKey
offsetProjScaleAtNatOriginGeoKey = geo_keys[i].value_offset;
break;
case 4096: // VerticalCSTypeGeoKey
set_VerticalCSTypeGeoKey(geo_keys[i].value_offset);
break;
case 4099: // VerticalUnitsGeoKey
set_VerticalUnitsGeoKey(geo_keys[i].value_offset);
break;
}
}
if (ellipsoid != -1)
{
set_reference_ellipsoid(ellipsoid);
}
if (!has_projection)
{
if (user_defined_projection == 1)
{
if ((offsetProjFalseEastingGeoKey >= 0) &&
(offsetProjFalseNorthingGeoKey >= 0) &&
(offsetProjNatOriginLatGeoKey >= 0) &&
(offsetProjCenterLongGeoKey >= 0) &&
(offsetProjScaleAtNatOriginGeoKey >= 0))
{
double falseEastingMeter =
geo_double_params[offsetProjFalseEastingGeoKey] * coordinates2meter;
double falseNorthingMeter =
geo_double_params[offsetProjFalseNorthingGeoKey] * coordinates2meter;
double latOriginDeg = geo_double_params[offsetProjNatOriginLatGeoKey];
double longMeridianDeg = geo_double_params[offsetProjCenterLongGeoKey];
double scaleFactor =
geo_double_params[offsetProjScaleAtNatOriginGeoKey];
set_transverse_mercator_projection(falseEastingMeter,
falseNorthingMeter, latOriginDeg, longMeridianDeg, scaleFactor);
if (description)
{
sprintf(description, "generic traverse mercator");
}
has_projection = true;
}
}
else if (user_defined_projection == 8)
{
if ((offsetProjFalseEastingGeoKey >= 0) &&
(offsetProjFalseNorthingGeoKey >= 0) &&
(offsetProjNatOriginLatGeoKey >= 0) &&
(offsetProjCenterLongGeoKey >= 0) &&
(offsetProjStdParallel1GeoKey >= 0) &&
(offsetProjStdParallel2GeoKey >= 0))
{
double falseEastingMeter =
geo_double_params[offsetProjFalseEastingGeoKey] * coordinates2meter;
double falseNorthingMeter =
geo_double_params[offsetProjFalseNorthingGeoKey] * coordinates2meter;
double latOriginDeg = geo_double_params[offsetProjNatOriginLatGeoKey];
double longOriginDeg = geo_double_params[offsetProjCenterLongGeoKey];
double firstStdParallelDeg =
geo_double_params[offsetProjStdParallel1GeoKey];
double secondStdParallelDeg =
geo_double_params[offsetProjStdParallel2GeoKey];
set_lambert_conformal_conic_projection(falseEastingMeter,
falseNorthingMeter, latOriginDeg, longOriginDeg, firstStdParallelDeg,
secondStdParallelDeg);
if (description)
{
sprintf(description, "generic lambert conformal conic");
}
has_projection = true;
}
}
}
return has_projection;
}
More information about the Liblas-devel
mailing list