Как установить гео-тэг GeogPrimeMeridianGeoKey

[TigerPVR]

Здравствуйте!
Создаю с помощью GDAL на C#.net файл GeoTiff.
Никак не получается сделать так, чтобы в выходном файле прописался гео-тэг:
GeogPrimeMeridianGeoKey (Short,1): PM_Greenwich

Гео-данные пишу в файл следующим образом:
oSrs.SetUTM(30, 1);
Double convertToRadians = Double.Parse(Osr.SRS_UA_DEGREE_CONV,
CultureInfo.InvariantCulture.NumberFormat);
oSrs.SetGeogCS(string.Empty, //"UTM 30 (WGS84) in northern hemisphere.",
Osr.SRS_DN_WGS84, Osr.SRS_DN_WGS84,
Osr.SRS_WGS84_SEMIMAJOR, Osr.SRS_WGS84_INVFLATTENING,
"Greenwich", 0.0, "degree", convertToRadians);

В чем может быть загвозка?

[Дмитрий Барышников]

У вас GDAL (gdalinfo) читает такой файл правильно?

[TigerPVR]

У вас GDAL (gdalinfo) читает такой файл правильно?

Кажется, да:

Код: Выделить всё

Driver: GTiff/GeoTIFF
Files: T0001S001_L1B_ver5.tif
Size is 12702, 7480
Coordinate System is:
PROJCS["WGS 84 / UTM zone 36N",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",33],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","32636"]]
Origin = (6659180.000000000000000,6404044.000000000000000)
Pixel Size = (5.884479999542236,-5.884479999542236)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_COPYRIGHT=(C) COPYRIGHT 2012 GEOMIR
  TIFFTAG_DATETIME=18.06.2012 10:50:41
  TIFFTAG_DOCUMENTNAME=PAN1A.120515_123257
  TIFFTAG_IMAGEDESCRIPTION=Geo-referenced image converted from HDF5
  TIFFTAG_MAXSAMPLEVALUE=4095
  TIFFTAG_MINSAMPLEVALUE=0
  TIFFTAG_SOFTWARE=GeoConverter
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( 6659180.000, 6404044.000) (146d 2' 2.59"W,102d 2'50.95"N)
Lower Left  ( 6659180.000, 6360028.090) (124d41'27.45"W, 95d59'40.78"N)
Upper Right ( 6733924.665, 6404044.000) (169d 6'37.63"W,109d10'51.21"N)
Lower Right ( 6733924.665, 6360028.090) (145d43' 7.88"W,102d29'27.38"N)
Center      ( 6696552.332, 6382036.045) (145d52'33.41"W,102d16' 3.17"N)
Band 1 Block=12702x90 Type=UInt16, ColorInterp=Gray

А вот необходимые мне ГЕО- теги не устанавливаются.
И почему-то некоторые из них прописаны не верно.
Например, должно быть: GeogCitationGeoKey (Ascii,6): "WGS84"
А получается: GeogCitationGeoKey (Ascii,7): "WGS 84" (с пробелом).

[Дмитрий Барышников]

А вот необходимые мне ГЕО- теги не устанавливаются.
И почему-то некоторые из них прописаны не верно.
Например, должно быть: GeogCitationGeoKey (Ascii,6): "WGS84"
А получается: GeogCitationGeoKey (Ascii,7): "WGS 84" (с пробелом).

А где это написано, что должно быть WGS84 без пробела?
По поводу не устанавливаемых геотегов - я уже писал, что GDAL представляет обобщенную модель доступа к растровым данным в различных форматах. Видимо в рамках такого подхода не все теги возможно вписать. Тоже хотелось бы почитать что-то типа стандарта, где прописаны обязательные теги для TIFF.

А вариант с SetMetadata не помог?

[TigerPVR]

Суть проблемы вот в чем.
Нужно создать GeoTIFF файл с привязкой. Файл с помощью GDAL создается нормально, но отчего-то он не правильно открывается в геодезической системе координат. Неправильно рассчитываются координаты в градусах.

У меня есть эталонный файл, с которым все ОК:

Код: Выделить всё

Geotiff_Information:
   Version: 1
   Key_Revision: 1.0
   Tagged_Information:
      ModelTiepointTag (2,3):
         0.5              0.5              0                
         489799.5         4802716.5        0                
      ModelPixelScaleTag (1,3):
         0.6              0.6              0                
      End_Of_Tags.
   Keyed_Information:
      GTModelTypeGeoKey (Short,1): ModelTypeProjected
      GTRasterTypeGeoKey (Short,1): RasterPixelIsArea
      GTCitationGeoKey (Ascii,24): "Projected Coordinates |"
      GeographicTypeGeoKey (Short,1): GCS_WGS_84
      GeogCitationGeoKey (Ascii,6): "WGS84"
      GeogPrimeMeridianGeoKey (Short,1): PM_Greenwich
      GeogLinearUnitsGeoKey (Short,1): Linear_Meter
      GeogAngularUnitsGeoKey (Short,1): Angular_Degree
      GeogEllipsoidGeoKey (Short,1): Ellipse_WGS_84
      GeogSemiMajorAxisGeoKey (Double,1): 6378137          
      GeogSemiMinorAxisGeoKey (Double,1): 6356752.31       
      ProjectedCSTypeGeoKey (Short,1): PCS_WGS84_UTM_zone_30N
      PCSCitationGeoKey (Ascii,12): "UTM / WGS84"
      ProjectionGeoKey (Short,1): Proj_UTM_zone_30N
      ProjLinearUnitsGeoKey (Short,1): Linear_Meter
      End_Of_Keys.
   End_Of_Geotiff.

PCS = 32630 (name unknown)
Projection = 16030 ()
Projection Method: CT_TransverseMercator
   ProjNatOriginLatGeoKey: 0.000000 (  0d 0' 0.00"N)
   ProjNatOriginLongGeoKey: -3.000000 (  3d 0' 0.00"W)
   ProjScaleAtNatOriginGeoKey: 0.999600
   ProjFalseEastingGeoKey: 500000.000000 m
   ProjFalseNorthingGeoKey: 0.000000 m
GCS: 4326/WGS 84
Datum: 6326/World Geodetic System 1984
Ellipsoid: 7030/WGS 84 (6378137.00,6356752.31)
Prime Meridian: 8901/Greenwich (0.000000/  0d 0' 0.00"E)
Projection Linear Units: 9001/metre (1.000000m)

Corner Coordinates:
Upper Left    (  489799.200, 4802716.800)
Lower Left    (  489799.200, 4796419.200)
Upper Right   (  493291.200, 4802716.800)
Lower Right   (  493291.200, 4796419.200)
Center        (  491545.200, 4799568.000)

И мой пробеммный файл:

Код: Выделить всё

Geotiff_Information:
   Version: 1
   Key_Revision: 1.0
   Tagged_Information:
      ModelTiepointTag (2,3):
         0                0                0                
         6659180          6404044          0                
      ModelPixelScaleTag (1,3):
         5.88448          5.88448          0                
      End_Of_Tags.
   Keyed_Information:
      GTModelTypeGeoKey (Short,1): ModelTypeProjected
      GTRasterTypeGeoKey (Short,1): RasterPixelIsArea
      GTCitationGeoKey (Ascii,24): "WGS_1984 / UTM zone 36N"
      GeogCitationGeoKey (Ascii,7): "WGS 84"
      GeogAngularUnitsGeoKey (Short,1): Angular_Degree
      ProjectedCSTypeGeoKey (Short,1): PCS_WGS84_UTM_zone_36N
      ProjLinearUnitsGeoKey (Short,1): Linear_Meter
      End_Of_Keys.
   End_Of_Geotiff.

PCS = 32636 (name unknown)
Projection = 16036 ()
Projection Method: CT_TransverseMercator
   ProjNatOriginLatGeoKey: 0.000000 (  0d 0' 0.00"N)
   ProjNatOriginLongGeoKey: 33.000000 ( 33d 0' 0.00"E)
   ProjScaleAtNatOriginGeoKey: 0.999600
   ProjFalseEastingGeoKey: 500000.000000 m
   ProjFalseNorthingGeoKey: 0.000000 m
GCS: 4326/WGS 84
Datum: 6326/World Geodetic System 1984
Ellipsoid: 7030/WGS 84 (6378137.00,6356752.31)
Prime Meridian: 8901/Greenwich (0.000000/  0d 0' 0.00"E)
Projection Linear Units: 9001/metre (1.000000m)

Corner Coordinates:
Upper Left    ( 6659180.000, 6404044.000)
Lower Left    ( 6659180.000, 6360028.090)
Upper Right   ( 6733924.665, 6404044.000)
Lower Right   ( 6733924.665, 6360028.090)
Center        ( 6696552.332, 6382036.045)

Почему так происходит? Не хватает какого-то тэга?

[Дмитрий Барышников]

Повторяю свои вопросы: почему вы считаете файл эталонным- он соответствует какому то стандарту? Или переформулируйте вопрос: как мне при помощи gdal получить файл аналогичный верхнему.

[ Сообщение с мобильного устройства ]

А еще у вас у верхнего и нижнего файла разные зоны 30 и 36 поэтому и координаты разные.

Код: Выделить всё

oSrs.SetUTM(36, 1);

[Ariki]

По-моему, вам нужно или убрать шестерку перед первой координатой, или пользоваться проекцией с false easting=6500000 вместо 500000.
Шестерка, кстати, меня смущает. У вас точно UTM, а не проекция Гаусса---Крюгера?

[TigerPVR]

По-моему, вам нужно или убрать шестерку перед первой координатой

Скажите, пожалуйста, перед какой координатой? Вместо 6659180 поставить 659180?

[Ariki]

Да, и перед всеми остальными абсциссами. Хотя я сам что-то запутался в том, что у вас происходит. Почему в коде у вас 30 зона UTM, а в тэгах GeoTIFF возникает 36-я? В то же время номер зоны 6, приписанный к абсциссе, наводит на мысль, что координаты углов растра, использованные вами для привязки, были вовсе не в 36 зоне UTM, а в 6 зоне Гаусса---Крюгера. Географически эти зоны совпадают (осевой меридиан один), но параметры проекции (ProjScale и FalseEasting) разные, а значит, привязка все равно будет неверной.

[TigerPVR]

Не получилось. И еще, при импорте правильного файла в EWI в MapInfo там стоит:
Proj: UTM, Zone 30N
UL Geo: 3'733.32''W, 43'22''38.10''N
UL Map: 489799.200, 4802716.800

А у GDAL-овского Proj какой-то странный, и градусы неверные:
Proj: GEOTIFF (Trandverse Mercator)
UL Geo: 16'28'21.48''W, 337'45'37.57''N
UL Map: 659180.000, 36404044.000

Параметры последней версии GDAL-овского файла:

Код: Выделить всё

Geotiff_Information:
   Version: 1
   Key_Revision: 1.0
   Tagged_Information:
      ModelTiepointTag (2,3):
         0                0                0                
         659180           36404044         0                
      ModelPixelScaleTag (1,3):
         5.88448          5.88448          0                
      End_Of_Tags.
   Keyed_Information:
      GTModelTypeGeoKey (Short,1): ModelTypeProjected
      GTRasterTypeGeoKey (Short,1): RasterPixelIsArea
      GTCitationGeoKey (Ascii,24): "WGS_1984 / UTM zone 36N"
      GeographicTypeGeoKey (Short,1): GCS_WGS_84
      GeogCitationGeoKey (Ascii,7): "WGS 84"
      GeogAngularUnitsGeoKey (Short,1): Angular_Degree
      GeogSemiMajorAxisGeoKey (Double,1): 6378137          
      GeogInvFlatteningGeoKey (Double,1): 298.257224       
      ProjectedCSTypeGeoKey (Short,1): User-Defined
      ProjectionGeoKey (Short,1): User-Defined
      ProjCoordTransGeoKey (Short,1): CT_TransverseMercator
      ProjLinearUnitsGeoKey (Short,1): Linear_Meter
      ProjNatOriginLongGeoKey (Double,1): 33               
      ProjNatOriginLatGeoKey (Double,1): 0                
      ProjFalseEastingGeoKey (Double,1): 6500000          
      ProjFalseNorthingGeoKey (Double,1): 0                
      ProjScaleAtNatOriginGeoKey (Double,1): 0.9996           
      End_Of_Keys.
   End_Of_Geotiff.

Projection Method: CT_TransverseMercator
   ProjNatOriginLatGeoKey: 0.000000 (  0d 0' 0.00"N)
   ProjNatOriginLongGeoKey: 33.000000 ( 33d 0' 0.00"E)
   ProjScaleAtNatOriginGeoKey: 0.999600
   ProjFalseEastingGeoKey: 6500000.000000 m
   ProjFalseNorthingGeoKey: 0.000000 m
GCS: 4326/WGS 84
Datum: 6326/World Geodetic System 1984
Ellipsoid: 7030/WGS 84 (6378137.00,6356752.31)
Prime Meridian: 8901/Greenwich (0.000000/  0d 0' 0.00"E)
Projection Linear Units: 9001/metre (1.000000m)

Corner Coordinates:
Upper Left    (  659180.000,36404044.000)
Lower Left    (  659180.000,36360028.090)
Upper Right   (  733924.665,36404044.000)
Lower Right   (  733924.665,36360028.090)
Center        (  696552.332,36382036.045)

[TigerPVR]

Почему в коде у вас 30 зона UTM, а в тэгах GeoTIFF возникает 36-я?

Эталонный файл в 30-й зоне, а мой - в 36. Это разные снимки. Просто эталонный открывается нормально, а мой - нет.

[TigerPVR]

Кажется, я понял, в чем проблема.
В правильном файле этот гео-тэг установлен: ProjectedCSTypeGeoKey (Short,1): PCS_WGS84_UTM_zone_30N
А в моем: ProjectedCSTypeGeoKey (Short,1): User-Defined

Почему так? Как его установить парвильно?

[Ariki]

Вы меня запутали еще больше. Или я вас...
Под абсциссой я здесь понимаю координату, направленную на восток. В каком порядке их GDAL показывает - не могу разобраться, может быть, и наоборот.
Если устанавливать 36 зону UTM, то False Easting для нее будет +500000, а перед координатой, направленной на восток, не должно быть номера зоны. Если же False Easting +36500000, то и перед координатой надо ставить 36. В проекции Гаусса---Крюгера (6 зона) будет соответственно +6500000 и шестерка перед координатой, а масштабный коэффициент равен 1 вместо 0.9996.

[Ariki]

И вообще, я подозреваю, что GDAL просто не понимает, как можно использовать SetGeogCS() с UTM. Этот метод вроде бы для географических систем координат, а углы у вас заданы в проекции UTM или что там у вас на самом деле.

[Ariki]

Мне кажется, вам нужен SetProjection(), которому можно передать параметры проекции в формате WKT или proj4. Когда-то использовал его (правда, из Python, а не C#) - проблем вроде не возникало.