Opened 13 years ago
Closed 13 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)
Change History (14)
comment:1 by , 13 years ago
by , 13 years ago
Attachment: | L71092086_08620110129_B50-small.TIF added |
---|
Smaller portion of the Geotiff dataset
by , 13 years ago
Attachment: | L71092086_08620110129_B50-netCDF-small.nc added |
---|
netCDF generated from the small Geotiff dataset
comment:2 by , 13 years ago
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 by , 13 years ago
Cc: | 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
by , 13 years ago
Attachment: | netcdf_fix4200.patch added |
---|
comment:4 by , 13 years ago
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:6 by , 13 years ago
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 by , 13 years ago
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 by , 13 years ago
Resolution: | → fixed |
---|---|
Status: | new → closed |
by , 12 years ago
Attachment: | tile_x5_y2.tif added |
---|
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
The gdalinfo for both Geotiff and netCDF files have been attached.