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;
}
_______________________________________________
Liblas-devel mailing list
[email protected]
http://lists.osgeo.org/mailman/listinfo/liblas-devel

Reply via email to