I'm working on a project in gdal where I'm generating standardized, Height-Above-Elipsoid(HAE)-elevation DEMs in UTM. Since pretty much all of my source data is in Mean Sea Level (MSL), that means I need to convert the vertical datum from MSL to HAE. The easy way to do this is by going to this site (https://www.ngs.noaa.gov/PC_PROD/GEOID09/dnpc09u.shtml) and downloading the GEOID09 tile that covers the area that you want to convert. I've been using Global Mapper to do this in the past, and have had no issues with either the download or the conversion, but I want to move away from using GM for other reasons.
To start the process, I try to slim down the data that I'm working on (each GEOID09 tile is huge, and there's no need to have all of the extra data that I don't want in there). To do this, I try to use
gdalwarp to crop the file down to within certain lat/long extents. I run this command:
subprocess.call("gdalwarp GEOID09/g2009u07.bin tmp/cutGeoid.tif -r bilinear -overwrite -te %d %d %d %d -s_srs %s -t_srs %s" %(corners, corners, corners, corners)
In this, "corners" is an array of the corner coordinates that I want to cut to, in lat/long. However, the GEOID09 .bin files use a type of longitude where they don't have "East" and "West", or positive and negative values on either side of the Prime Meridian; they instead go from 0 to 180 from the Prime Meridian east, out to the International Date Line, then continue from 181 to 360 from there eastward across the Americas and back to the Prime. Turning those back into "normal" coordinates is trivial, but since the .bin uses this "positives only" format, it doesn't understand if you try to crop using "normal" coordinates. Finding this, I simply cropped using the "positive only" coordinate system, which it was fine with. Then, just to check that it had done it properly, I loaded both the original and the cropped files into Global Mapper, and discovered that Global Mapper now thought that the cropped image wasn't converting the coordinates anymore; it was putting the cropped image where it thought the longitude line '234 East' should be, waaaayyy off to the side, while the original .bin was still right where it should be at '-93 East' or '93 West'.
Now very confused, I looked at the metadata for both files using
gdalinfo to see what it thought the differences between the systems were. The .bin file's metadata read as this:
Driver: NGSGEOID/NOAA NGS Geoid Height Grids Files: GEOID09/g2009u07.bin Size is 1141, 1081 Coordinate System is: GEOGCS["WGS 84", DATUM["WGS_1984", SPHEROID["WGS 84",6378137,298.257223563, AUTHORITY["EPSG","7030"]], TOWGS84[0,0,0,0,0,0,0], AUTHORITY["EPSG","6326"]], PRIMEM["Greenwich",0, AUTHORITY["EPSG","8901"]], UNIT["degree",0.0174532925199433, AUTHORITY["EPSG","9108"]], AUTHORITY["EPSG","4326"]] Origin = (263.991666666666674,42.008333333336942) Pixel Size = (0.016666666666670,-0.016666666666670) Corner Coordinates: Upper Left ( 263.992, 42.008) (263d59'30.00"E, 42d 0'30.00"N) Lower Left ( 263.992, 23.992) (263d59'30.00"E, 23d59'30.00"N) Upper Right ( 283.008, 42.008) (283d 0'30.00"E, 42d 0'30.00"N) Lower Right ( 283.008, 23.992) (283d 0'30.00"E, 23d59'30.00"N) Center ( 273.500, 33.000) (273d30' 0.00"E, 33d 0' 0.00"N) Band 1 Block=1141x1 Type=Float32, ColorInterp=Undefined Unit Type: m
While the metadata for the .tif I had generated with
gdalwarp read as such
Driver: GTiff/GeoTIFF Files: tmp/cutGeoid.tif Size is 120, 120 Coordinate System is: 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"]] Origin = (266.000000000000000,38.000000000000000) Pixel Size = (0.016666666666667,-0.016666666666667) Metadata: AREA_OR_POINT=Area Image Structure Metadata: INTERLEAVE=BAND Corner Coordinates: Upper Left ( 266.000, 38.000) (266d 0' 0.00"E, 38d 0' 0.00"N) Lower Left ( 266.000, 36.000) (266d 0' 0.00"E, 36d 0' 0.00"N) Upper Right ( 268.000, 38.000) (268d 0' 0.00"E, 38d 0' 0.00"N) Lower Right ( 268.000, 36.000) (268d 0' 0.00"E, 36d 0' 0.00"N) Center ( 267.000, 37.000) (267d 0' 0.00"E, 37d 0' 0.00"N) Band 1 Block=120x17 Type=Float32, ColorInterp=Gray Unit Type: m
You'll notice that the coordinate system has been altered - the .tif is missing some of the EPSG definitions that the .bin had. So far as I can tell, something that it lost there was what allowed Global Mapper to parse the longitude data correctly, but I can't figure out what it could have been. Has anyone ever encountered this problem before? Is there something unique about the GEOID09 data that gdal can't figure out or throws away?
EDIT: I believe I've found why the altered metadata is causing problems; EPSG 9108 defines degrees specifically in integers, which would explain why the files with it in their metadata can understand the longitude conversion while the files without it can't. However, if that is the case, why is it being removed at all? It seems like it would be a pretty important piece of information to have, yet
gdalwarp trims it off.