Opened 6 years ago

Closed 6 years ago

#4200 closed defect (fixed)

[PATCH] gdal_translate GeoTIFF to NetCDF results in coordinates shifting

Reported by: hsumanto Owned by: warmerdam
Priority: normal Milestone:
Component: default Version: unspecified
Severity: normal Keywords: netcdf, gdal_translate
Cc: Kyle Shannon

Description

Converting a WELD LANDSAT (in Transverse Mercator projection) Geotiff into NetCDF using gdal_translate (1.8.1) results in origin and corner coordinates shifting.

gdal_translate -of netCDF L71092086_08620110129_B50.TIF L71092086_08620110129_B50-netCDF.nc

The origin has shifted from Origin = (281384.999999999941792,-4040684.999999999534339) to Origin = (281385.000000000000000,-4040680.000000000000000)

The corner coordinates have shifted from Corner Coordinates: Upper Left ( 281385.000,-4040685.000) Lower Left ( 281385.000,-4254315.000) Upper Right ( 525015.000,-4040685.000) Lower Right ( 525015.000,-4254315.000) Center ( 403200.000,-4147500.000) to Corner Coordinates: Upper Left ( 281385.000,-4040680.000) Lower Left ( 281385.000,-4254310.000) Upper Right ( 525015.000,-4040680.000) Lower Right ( 525015.000,-4254310.000) Center ( 403200.000,-4147495.000)

You can easily notice around 5 metres shift in Y value.

Any idea what cause this shifting?

Even when I tried to assign/force the original Geotiff corner bounds to the output file using -a_ullr ulx uly lrx lry as below, gdal_translate -of netCDF -a_ullr 281385.000 -4040685.000 525015.000 -4254315.000 L71092086_08620110129_B50.TIF L71092086_08620110129_B50-forcedbound.nc the output netcdf file still have the origin and corner coordinates shifted (not using the assigned corner bounds.

================================================================== gdalinfo L71092086_08620110129_B50.TIF

Driver: GTiff/GeoTIFF Files: L71092086_08620110129_B50.TIF Size is 8121, 7121 Coordinate System is: PROJCS["WGS 84 / UTM zone 55N",

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"]],

PROJECTIONTransverse_Mercator?, PARAMETER["latitude_of_origin",0], PARAMETER["central_meridian",147], PARAMETER["scale_factor",0.9996], PARAMETER["false_easting",500000], PARAMETER["false_northing",0], UNIT["metre",1,

AUTHORITY["EPSG","9001"]],

AUTHORITY["EPSG","32655"]]

Origin = (281384.999999999941792,-4040684.999999999534339) Pixel Size = (30.000000000000000,-30.000000000000000) Metadata:

AREA_OR_POINT=Point

Image Structure Metadata:

INTERLEAVE=BAND

Corner Coordinates: Upper Left ( 281385.000,-4040685.000) (144d33'34.28"E, 36d29'11.56"S) Lower Left ( 281385.000,-4254315.000) (144d29'46.21"E, 38d24'37.27"S) Upper Right ( 525015.000,-4040685.000) (147d16'45.71"E, 36d30'40.26"S) Lower Right ( 525015.000,-4254315.000) (147d17'11.85"E, 38d26'12.33"S) Center ( 403200.000,-4147500.000) (145d54'19.15"E, 37d28' 9.55"S) Band 1 Block=8121x1 Type=Byte, ColorInterp?=Gray

================================================================== gdalinfo L71092086_08620110129_B50-netCDF.nc

Driver: netCDF/Network Common Data Format Files: L71092086_08620110129_B50-netCDF.nc

L71092086_08620110129_B50-netCDF.nc.aux.xml

Size is 8121, 7121 Coordinate System is: PROJCS["WGS 84 / UTM zone 55N",

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"]],

PROJECTIONTransverse_Mercator?, PARAMETER["latitude_of_origin",0], PARAMETER["central_meridian",147], PARAMETER["scale_factor",0.9996], PARAMETER["false_easting",500000], PARAMETER["false_northing",0], UNIT["metre",1,

AUTHORITY["EPSG","9001"]],

AUTHORITY["EPSG","32655"]]

Origin = (281385.000000000000000,-4040680.000000000000000) Pixel Size = (30.000000000000000,-30.000000000000000) Metadata:

NC_GLOBAL#Conventions=CF-1.0 Band1#grid_mapping=transverse_mercator Band1#long_name=GDAL Band Number 1 transverse_mercator#Northernmost_Northing=-4.04068e+06 transverse_mercator#Southernmost_Northing=-4.25432e+06 transverse_mercator#Easternmost_Easting=525015 transverse_mercator#Westernmost_Easting=281385 transverse_mercator#spatial_ref=PROJCS["WGS 84 / UTM zone 55N",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"]],PROJECTIONTransverse_Mercator?,PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",147],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","32655"]] transverse_mercator#GeoTransform?=281385 30 0 -4.04068e+06 0 -30 transverse_mercator#grid_mapping_name=transverse_mercator transverse_mercator#latitude_of_projection_origin=0.000000e+00 transverse_mercator#longitude_of_central_meridian=1.470000e+02 transverse_mercator#scale_factor_at_central_meridian=9.996000e-01 transverse_mercator#false_easting=5.000000e+05 transverse_mercator#false_northing=0.000000e+00 AREA_OR_POINT=Point

Corner Coordinates: Upper Left ( 281385.000,-4040680.000) (144d33'34.28"E, 36d29'11.40"S) Lower Left ( 281385.000,-4254310.000) (144d29'46.22"E, 38d24'37.11"S) Upper Right ( 525015.000,-4040680.000) (147d16'45.71"E, 36d30'40.10"S) Lower Right ( 525015.000,-4254310.000) (147d17'11.85"E, 38d26'12.17"S) Center ( 403200.000,-4147495.000) (145d54'19.15"E, 37d28' 9.38"S) Band 1 Block=8121x1 Type=Byte, ColorInterp?=Gray

NoData? Value=0 Metadata:

NETCDF_VARNAME=Band1

Attachments (6)

L71092086_08620110129_B50.gdalinfo (1.3 KB) - added by hsumanto 6 years ago.
Geotiff file gdalinfo
L71092086_08620110129_B50-netCDF.gdalinfo (2.7 KB) - added by hsumanto 6 years ago.
netCDF gdalinfo
L71092086_08620110129_B50-small.TIF (108.8 KB) - added by hsumanto 6 years ago.
Smaller portion of the Geotiff dataset
L71092086_08620110129_B50-netCDF-small.nc (109.6 KB) - added by hsumanto 6 years ago.
netCDF generated from the small Geotiff dataset
netcdf_fix4200.patch (749 bytes) - added by Even Rouault 6 years ago.
tile_x5_y2.tif (588.7 KB) - added by hsumanto 5 years ago.

Download all attachments as: .zip

Change History (14)

comment:1 Changed 6 years ago by hsumanto

I initially want to attach the GeoTIFF file but as it is 56MB I don't think I could attach in this ticket as an attachment but it is just a free LANDSAT data in GeoTIFF I got from USGS using

WRS-2 Path: 92 Row: 86 (Australia scene) from http://glovis.usgs.gov/

The gdalinfo for both Geotiff and netCDF files have been attached.

Changed 6 years ago by hsumanto

Geotiff file gdalinfo

Changed 6 years ago by hsumanto

netCDF gdalinfo

Changed 6 years ago by hsumanto

Smaller portion of the Geotiff dataset

Changed 6 years ago by hsumanto

netCDF generated from the small Geotiff dataset

comment:2 Changed 6 years ago by hsumanto

I have taken a smaller portion the dataset by doing

gdal_translate -projwin 281385.000 -4040685.000 291385.000 -4050685.000 -a_ullr 281385.000 -4040685.000 291385.000 -4050685.000 L71092086_08620110129_B50.TIF L71092086_08620110129_B50-small.TIF

and still can reproduce the same issue, so I guess we could just use this smaller dataset (L71092086_08620110129_B50-small.TIF) as testing input file.

comment:3 Changed 6 years ago by Even Rouault

Cc: Kyle Shannon added
Summary: gdal_translate GeoTIFF to NetCDF results in coordinates shifting[PATCH] gdal_translate GeoTIFF to NetCDF results in coordinates shifting

Well, the problem is an insufficient precision when writting the geotransform:

  transverse_mercator#GeoTransform=281385 30.03 0 -4.04068e+06 0 -30.03

due to a sprintf("%g") when formatting it. This can easily be fixed by the attached patch that forces the precision to 16 significant digits.

However, the nc_put_att_double() call used when writing the Northernmost_Northing and Southernmost_Northing parameters also seem to loose precision, but we don't have control on that code, so I'm not sure what could be done. Perhaps using nc_put_att_text() instead and formatting ourself ? But I've no idea if there's strong typing in NetCDF metadata, so I'll let that to the NetCDF specialists.

  transverse_mercator#Northernmost_Northing=-4.04068e+06
  transverse_mercator#Southernmost_Northing=-4.05068e+06

Changed 6 years ago by Even Rouault

Attachment: netcdf_fix4200.patch added

comment:4 Changed 6 years ago by hsumanto

Thank you kyle for your patch.

I have tried the patch, it definitely fixed the coordinates shifting caused by insufficient precision.

I tried to confirm your suspicion that the Northernmost_Northing and Southernmost_Northing parameters also seem to loose precision by looking directly through the newly created netcdf file using ncdump as shown below.

They don't actually loose any precision in the actual netcdf file but I wonder why the gdalinfo still print them with insufficient precision, my guess is it might still sprintf using "%g " format thus results (the full gdalinfo is at the bottom).

  transverse_mercator#Northernmost_Northing=-4.04068e+06
  transverse_mercator#Southernmost_Northing=-4.25432e+06

Hopefully the Northernmost_Northing and Southernmost_Northing parameters could be fixed too.

Will this patch (netcdf_fix4200.patch) be committed and released? If yes when or which release version will it be?

netcdf C:/tmp/L71092086_08620110129_B50-p1.nc {
 dimensions:
   x = 8121;
   y = 7121;
 variables:
   char transverse_mercator;
     :Northernmost_Northing = -4040684.9999999995; // double
     :Southernmost_Northing = -4254315.0; // double
     :Easternmost_Easting = 525015.0; // double
     :Westernmost_Easting = 281384.99999999994; // double
     :spatial_ref = "PROJCS[\"WGS 84 / UTM zone 55N\",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\",147],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",500000],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AUTHORITY[\"EPSG\",\"32655\"]]";
     :GeoTransform = "281384.9999999999 30 0 -4040685 0 -30 ";
     :grid_mapping_name = "transverse_mercator";
     :latitude_of_projection_origin = 0.0f; // float
     :longitude_of_central_meridian = 147.0f; // float
     :scale_factor_at_central_meridian = 0.9996f; // float
     :false_easting = 500000.0f; // float
     :false_northing = 0.0f; // float
   byte Band1(y=7121, x=8121);
     :grid_mapping = "transverse_mercator";
     :long_name = "GDAL Band Number 1";

 :Conventions = "CF-1.0";
}
Driver: netCDF/Network Common Data Format
Files: L71092086_08620110129_B50-p1.nc
       L71092086_08620110129_B50-p1.nc.aux.xml
Size is 8121, 7121
Coordinate System is:
PROJCS["WGS 84 / UTM zone 55N",
    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",147],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","32655"]]
Origin = (281384.999999999883585,-4040685.000000000000000)
Pixel Size = (30.000000000000000,-30.000000000000000)
Metadata:
  NC_GLOBAL#Conventions=CF-1.0
  Band1#grid_mapping=transverse_mercator
  Band1#long_name=GDAL Band Number 1
  transverse_mercator#Northernmost_Northing=-4.04068e+06
  transverse_mercator#Southernmost_Northing=-4.25432e+06
  transverse_mercator#Easternmost_Easting=525015
  transverse_mercator#Westernmost_Easting=281385
  transverse_mercator#spatial_ref=PROJCS["WGS 84 / UTM zone 55N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.017453292
  transverse_mercator#GeoTransform=281384.9999999999 30 0 -4040685 0 -30
  transverse_mercator#grid_mapping_name=transverse_mercator
  transverse_mercator#latitude_of_projection_origin=0.000000e+00
  transverse_mercator#longitude_of_central_meridian=1.470000e+02
  transverse_mercator#scale_factor_at_central_meridian=9.996000e-01
  transverse_mercator#false_easting=5.000000e+05
  transverse_mercator#false_northing=0.000000e+00
  AREA_OR_POINT=Point
Corner Coordinates:
Upper Left  (  281385.000,-4040685.000) (144d33'34.28"E, 36d29'11.56"S)
Lower Left  (  281385.000,-4254315.000) (144d29'46.21"E, 38d24'37.27"S)
Upper Right (  525015.000,-4040685.000) (147d16'45.71"E, 36d30'40.26"S)
Lower Right (  525015.000,-4254315.000) (147d17'11.85"E, 38d26'12.33"S)
Center      (  403200.000,-4147500.000) (145d54'19.15"E, 37d28' 9.55"S)
Band 1 Block=8121x1 Type=Byte, ColorInterp=Gray
  NoData Value=0
  Metadata:
    NETCDF_VARNAME=Band1

Regards, Hendy

comment:5 Changed 6 years ago by hsumanto

Thank you rouault.

comment:6 Changed 6 years ago by hsumanto

As for using nc_put_att_text() instead nc_put_att_double() call used when writing the Northernmost_Northing and Southernmost_Northing parameters, the unidata actually recommends the type safe versions of this function whenever possible, thus they are in favour of nc_put_att_double() for those parameters as the parameters are of type double.

http://www.unidata.ucar.edu/software/netcdf/docs/netcdf-c/nc_005fput_005fatt_005f-type.html

        status = nc_put_att_double( fpImage,
                                    NCDFVarID,
                                    "Northernmost_Northing",
                                    NC_DOUBLE,
                                    1,
                                    &dfNN );

comment:7 Changed 6 years ago by hsumanto

I tried mod-2-meta.txt patch (Ticket #2129) and it fixes the duplication of metadata and it also fixes the issue where Northernmost_Northing and Southernmost_Northing parameters also seem to loose precision.

Thank you for both rouault and etiennesky.

comment:8 Changed 6 years ago by Kyle Shannon

Resolution: fixed
Status: newclosed

I haven't been around, trying to catch up. Even, your patch looks good here, I applied and commited to the trunk in r22977 . I will pick through the other tickets (#2129 and #4204) as soon as I can. I just got a new job, so my time is a little sporadic. Thanks for your patience. Thanks again Even.

Changed 5 years ago by hsumanto

Attachment: tile_x5_y2.tif added
Note: See TracTickets for help on using tickets.