Hmm, turns out that when I asked nicely at my university they just gave me access to ArcGIS Pro, so I decided to try the transformation and re-calculation one more time using ArcGIS instead of GDAL, just to be sure...
And what do you know:
I finally got it to work!!! Okay, so here's what I did:
I started from the 32-bit floating-point DEM dataset (
DEM_LiDAR_1ft_2010_Improved_NYC.img).
As usual I started by projecting the raster from NAD83 to WGS84 UTM 18N. In ArcGIS I used the
Data Management Tools → Projections and Transformations → Raster → Project Raster tool to project the raster from NAD83 to WGS84 UTM 18N. This is a tool that is also available in ArcDesktop (using ArcMap og ArcCatalog), so you do not need full ArcGIS Pro for this. I googled a little and found
How To: Determine which NAD_1983_To_WGS_1984 transformation to use which suggests using the
WGS_1984_(ITRF00)_To_NAD_1983 transformation yielding a transformation accuracy of up to 0.10m. This leads me to the following parameters:
Code:
Parameters
Input Raster D:\Downloads\NYC Open Data\DEM_LiDAR_1ft_2010_Improved_NYC.img
Output Raster Dataset D:\Downloads\NYC Open Data\DEM_LiDAR_1ft_2010_Improved_NYC_wgs84.tif
Output Coordinate System PROJCS['WGS_1984_UTM_Zone_18N',GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Transverse_Mercator'],PARAMETER['False_Easting',500000.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',-75.0],PARAMETER['Scale_Factor',0.9996],PARAMETER['Latitude_Of_Origin',0.0],UNIT['Meter',1.0]]
Resampling Technique NEAREST
Output Cell Size 1 1
Geographic Transformation WGS_1984_(ITRF00)_To_NAD_1983
Registration Point
Input Coordinate System PROJCS['NAD_1983_StatePlane_New_York_Long_Island_FIPS_3104_Feet',GEOGCS['GCS_North_American_1983',DATUM['D_North_American_1983',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Lambert_Conformal_Conic'],PARAMETER['False_Easting',984250.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',-74.0],PARAMETER['Standard_Parallel_1',40.66666666666666],PARAMETER['Standard_Parallel_2',41.03333333333333],PARAMETER['Latitude_Of_Origin',40.16666666666666],UNIT['Foot_US',0.3048006096012192]]
Vertical NO_VERTICAL
Environments
Output Coordinate System PROJCS['WGS_1984_UTM_Zone_18N',GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Transverse_Mercator'],PARAMETER['False_Easting',500000.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',-75.0],PARAMETER['Scale_Factor',0.9996],PARAMETER['Latitude_Of_Origin',0.0],UNIT['Meter',1.0]]
Pyramid NONE
Messages
Start Time: onsdag 24. oktober 2018 23.37.47
Succeeded at onsdag 24. oktober 2018 23.49.13 (Elapsed Time: 11 minutes 26 seconds)
Since we still need to transform the cell values from feet to metres, I disabled building pyramids, we won't render the output raster in ArcMap anyways, and building Pyramids can take quite a lot of time. I actually did this transform twice, using different settings for the Output cell size. ArcGIS originally suggest the cell size to be 1ft x 1ft (i.e. 0.3048006096012192 x 0.3048006096012192), since the original is 1ft by 1ft. However, keeping a cell size of 1sqft did not work out. So I simply sacrifice some accuracy and set the cell size to 1sqm instead. If I understand correctly, this will merge ~9 cells into one. The result will still be pretty accurate, with an accuracy of ~1m. It also results in a much smaller file, remember the original image is 96GB in size, while the output (uncompressed) will yield ~8.73GB.
Okay, now I applied the recalculation of all elevation values from US Survey Feet to Metres using the
Spatial Analyst Tools → Math → Times tool using the following parameters. Note that the
Math tools are usable only with an ArcGIS Pro license, so ArcGIS Desktop won't do here. But since the calculation here is relatively benign, I guess the
gdal_calc.py utility I used earlier would do the trick as well.
Code:
Parameters
Input raster or constant value 1 D:\Downloads\NYC Open Data\DEM_LiDAR_1ft_2010_Improved_NYC_wgs84.tif
Input raster or constant value 2 0,304800609601219
Output raster D:\Downloads\NYC Open Data\DEM_LiDAR_1ft_2010_Improved_NYC_wgs84_metres.tif
Messages
Start Time: onsdag 24. oktober 2018 23.49.14
Succeeded at onsdag 24. oktober 2018 23.53.40 (Elapsed Time: 4 minutes 26 seconds)
And again, we have to remember to set the vertical coordinate system. As before I use
Unknown_height_system_(meters) as the vertical coordinate system.
Okay, moment of truth: Loading the file in TranzDEM. And yes, it finally works!
So, I have these theories to why this worked and the previous attempts did not:
- TransDEM does not like cell sizes that are less than 1sqm in GeoTIFF
- The 1sqft file is so large (96GB) that TransDEM fails early when loading it, whereas the 9GB large 1sqm file is just small enough that TransDEM will accept it
- Since ArcGIS has a UI, it might be that I set better settings or that ArcGIS is smarter and prevents me from doing errors as opposed to the command-line tools of GDAL.
Anyways, when I export the DEM to T:ANE everything looks fine, so now I can start filling this up with OSM Map tiles.