NYC Open Data 1 foot Digital Elevation Model
NYC Open Data 1 foot Digital Elevation Model
Hi,
I recently found out that New York City publishes a NYC 1foot Digital Elevation Model (https://data.cityofnewyork.us/City-Gove ... /dpc8-z3jc). So, I though to myself, what if we could take this pretty accurate DEM and make a Trainz route out of it. The same website also publishes KML-files containing the routes of the New York Subway, so getting OSM tiles along the subways is pretty easy. But first to the DEM.
During this process I encountered a series of problems, the following posts document what I have done and any input to my process is highly appreciated. I'll show the properties of the files using gdalinfo. I am currently stuck and do not know how to continue with this project.
Noob Disclaimer: This is my first time dealing with DEMs of any kind, so please bear with me for anything stupid I do. I am a Computer Scientist though, so I am pretty comfortable with command-line tools, different file formats and datatypes. Yes, starting with such a big file on such a large area might not be the best way to get started, but hey: Real World scenario. I did go through the Müngsten tutorial and that worked without any problems.
I recently found out that New York City publishes a NYC 1foot Digital Elevation Model (https://data.cityofnewyork.us/City-Gove ... /dpc8-z3jc). So, I though to myself, what if we could take this pretty accurate DEM and make a Trainz route out of it. The same website also publishes KML-files containing the routes of the New York Subway, so getting OSM tiles along the subways is pretty easy. But first to the DEM.
During this process I encountered a series of problems, the following posts document what I have done and any input to my process is highly appreciated. I'll show the properties of the files using gdalinfo. I am currently stuck and do not know how to continue with this project.
Noob Disclaimer: This is my first time dealing with DEMs of any kind, so please bear with me for anything stupid I do. I am a Computer Scientist though, so I am pretty comfortable with command-line tools, different file formats and datatypes. Yes, starting with such a big file on such a large area might not be the best way to get started, but hey: Real World scenario. I did go through the Müngsten tutorial and that worked without any problems.
Re: NYC Open Data 1 foot Digital Elevation Model
I downloaded the 1 foot Digital Elevation Model from the NYC Open Data website (link above). The file is 26 GB large, so I downloaded it while I was on the University Campus and had access to a high-speed internet connection.
After unzipping the files (inflating them to ~121 GB of data!), I tried to open the DEM in TranzDEM. TranzDEM simply throws a very unhelpful error message:
Okay, so let's look at the file with some other tools. (I have omitted irrelevant histogram binary data and the GDALRasterAttributeTable information from the output from gdalinfo below to keep the post short.)
Can anyone spot a reason why TranzDEM complains? I also opened the file in ArcMap (ArcGIS Desktop 10.6) and the DEM renders beautifully with amazing accuracy.
After unzipping the files (inflating them to ~121 GB of data!), I tried to open the DEM in TranzDEM. TranzDEM simply throws a very unhelpful error message:
Code: Select all
Error reading DEM!Code: Select all
> gdalinfo DEM_LiDAR_1ft_2010_Improved_NYC.img
Driver: HFA/Erdas Imagine Images (.img)
Files: DEM_LiDAR_1ft_2010_Improved_NYC.img
DEM_LiDAR_1ft_2010_Improved_NYC.ige
DEM_LiDAR_1ft_2010_Improved_NYC.rrd
DEM_LiDAR_1ft_2010_Improved_NYC.rde
Size is 158100, 156100
Coordinate System is:
PROJCS[
PROJECTION["Lambert_Conformal_Conic_2SP"],
UNIT["Meter",1]]
Origin = (910719.300000000046566,275160.674999999988358)
Pixel Size = (1.000000000000000,-1.000000000000000)
Corner Coordinates:
Upper Left ( 910719.300, 275160.675)
Lower Left ( 910719.300, 119060.675)
Upper Right ( 1068819.300, 275160.675)
Lower Right ( 1068819.300, 119060.675)
Center ( 989769.300, 197110.675)
Band 1 Block=512x512 Type=Float32, ColorInterp=Undefined
Description = Layer_1
Min=-73.828 Max=411.272
Minimum=-73.828, Maximum=411.272, Mean=40.993, StdDev=50.695
Overviews: 79050x78050, 39525x39025, 19762x19512, 9881x9756, 4940x4878, 2470x2439, 1235x1219, 617x609, 308x304, 154x152, 77x76, 38x38
Metadata:
LAYER_TYPE=athematic
STATISTICS_EXCLUDEDVALUES=0
STATISTICS_HISTOMAX=411.27249145508
STATISTICS_HISTOMIN=-73.827598571777
STATISTICS_HISTONUMBINS=256
STATISTICS_MAXIMUM=411.27249145508
STATISTICS_MEAN=40.993348964296
STATISTICS_MEDIAN=22.81343498826
STATISTICS_MINIMUM=-73.827598571777
STATISTICS_MODE=-3.7154761850834
STATISTICS_SKIPFACTORX=33
STATISTICS_SKIPFACTORY=33
STATISTICS_STDDEV=50.695343773293
Re: NYC Open Data 1 foot Digital Elevation Model
So since the ERDAS Imagine file did not work in TranzDEM I decided to try the alternate file that is published on the NYC Open Data website: 1 foot Digital Elevation Model (DEM) Integer Raster (https://data.cityofnewyork.us/City-Gove ... 7kuu-zah7/).
The file is also significantly smaller (~3.5 GB), obviously sacrificing precision by storing the elevation as 16-bit integers instead of 32-bit floats. Still, the file renders beautifully in ArcMap, so let's give it a try.
TranzDEM complains:
If I read this error message correctly, TranzDEM expects the input data to be metric instead of imperial. And sure enough according to both ArcGIS and GDAL the units used here are feet:
The file is also significantly smaller (~3.5 GB), obviously sacrificing precision by storing the elevation as 16-bit integers instead of 32-bit floats. Still, the file renders beautifully in ArcMap, so let's give it a try.
TranzDEM complains:
Code: Select all
ProjLinearUnitsGeoKey tag not Linear_Meter [Linear_Foot_US_Survey]Code: Select all
> gdalinfo DEM_LiDAR_1ft_2010_Improved_NYC_int.tif
Driver: GTiff/GeoTIFF
Files: DEM_LiDAR_1ft_2010_Improved_NYC_int.tif
DEM_LiDAR_1ft_2010_Improved_NYC_int.tif.aux.xml
Size is 158100, 156100
Coordinate System is:
PROJCS["NAD_1983_StatePlane_New_York_Long_Island_FIPS_3104_Feet",
GEOGCS["NAD83",
DATUM["North_American_Datum_1983",
SPHEROID["GRS 1980",6378137,298.2572221010042,
AUTHORITY["EPSG","7019"]],
AUTHORITY["EPSG","6269"]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433],
AUTHORITY["EPSG","4269"]],
PROJECTION["Lambert_Conformal_Conic_2SP"],
PARAMETER["standard_parallel_1",41.03333333333333],
PARAMETER["standard_parallel_2",40.66666666666666],
PARAMETER["latitude_of_origin",40.16666666666666],
PARAMETER["central_meridian",-74],
PARAMETER["false_easting",984250.0000000002],
PARAMETER["false_northing",0],
UNIT["US survey foot",0.3048006096012192,
AUTHORITY["EPSG","9003"]],
AUTHORITY["EPSG","2263"]]
Origin = (910719.300000000046566,275160.674999999988358)
Pixel Size = (1.000000000000000,-1.000000000000000)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
COMPRESSION=LZW
INTERLEAVE=BAND
Corner Coordinates:
Upper Left ( 910719.300, 275160.675) ( 74d15'57.85"W, 40d55'17.83"N)
Lower Left ( 910719.300, 119060.675) ( 74d15'51.71"W, 40d29'35.39"N)
Upper Right ( 1068819.300, 275160.675) ( 73d41'38.35"W, 40d55'17.47"N)
Lower Right ( 1068819.300, 119060.675) ( 73d41'45.41"W, 40d29'35.04"N)
Center ( 989769.300, 197110.675) ( 73d58'48.33"W, 40d42'27.71"N)
Band 1 Block=158100x7 Type=UInt16, ColorInterp=Gray
Description = Layer_1
Min=0.000 Max=411.000
Minimum=0.000, Maximum=411.000, Mean=18.435, StdDev=39.424
Overviews: 79050x78050, 39525x39025, 19763x19513, 9882x9757, 4941x4879, 2471x2440, 1236x1220, 618x610, 309x305, 155x153
Metadata:
DESCRIPTION=Layer_1
STATISTICS_COVARIANCES=1554.21862787829
STATISTICS_MAXIMUM=411
STATISTICS_MEAN=18.435255624912
STATISTICS_MINIMUM=0
STATISTICS_SKIPFACTORX=1
STATISTICS_SKIPFACTORY=1
STATISTICS_STDDEV=39.4235795924
Last edited by couven92 on 22 Oct 2018 22:51, edited 1 time in total.
Re: NYC Open Data 1 foot Digital Elevation Model
Since TranzDEM this time at least told us what is wrong, this times we can do something with the error. In the information above we can see that the GTiff uses a NAD83 projection with imperial units. So the file needs to be reprojected into a metrical coordinate system. I googled a little and found out that the UTM Zone 18N using WGS84 Datum is suitable for New York. TranzDEM is also more happy to accomodate UTM/WGS84 projections than NAD projections.
So, we need to use the gdalwarp utility to warp the raster to WGS84. According to http://www.spatialreference.org/ref/eps ... -zone-18n/ the proj4 string for UTM Zone 18N is:
Note that it actually says: units=m (i.e. metric units)
Using this, the gdalwarp command-line looks like this:
I added some extra parameters to use 4GB of memory for the cache and to use all CPUs in parallel for the warp. This speeds up the process considerably. I use --co to indicate a big TIFF file and I'd like to use LZW compression.
You could achieve the same using ArcGIS. I think the Project Raster tool is used for this. I did however find that ArcGIS is several magnitudes slower in the transformation.
Okay, so this takes care of the coordinates used in the file. However, this is a DEM, meaning we also have vertical units, i.e. the actual value at each pixel. These values are still in feet, as gdalwarp only projects along the x and y axes. All pixels need to be re-calculated from feet to metres. I finally found this answer after some googling: https://gis.stackexchange.com/a/259904
Form the previous post we can see that the format in the file defines a US Survey Foot as 0.3048006096012192 metres. Ergo, we have to multiply each pixel with 0.3048006096012192 using gdal_calc.py as suggested. Then we can use ArcGIS to change the metadata in the GeoTIFF file to say that we have a metric vertical units.
I ended up with the following command-line:
Note that I am using --type=Int16 to change the data type from unsigned to signed. This is because I on the way discovered that TranzDEM wants signed integers instead of unsigned integers. When looking at the original metadata the minimum value is 0ft and max is 411ft, so conversion to signed integers is lossless in this case.
BTW: As far as I understood ArcGIS would also be able to do the exact same operation using the Times operation in the Math Toolbox of the Spatial Analyst Tools. But you need an ArcGIS Pro license for that. So, I am not using that one.
Finally, TransDEM wants the DEM to actually specify the vertical units. Otherwise, it now fails with:
I haven't yet found out how to set that tag using GDAL, but in ArcGIS you very easily manipulate the metadata in the GeoTIFF file using ArcCatalog. I choose Unknown_height_system_(meters) for the Vertical Z Coordinate system.
So, we need to use the gdalwarp utility to warp the raster to WGS84. According to http://www.spatialreference.org/ref/eps ... -zone-18n/ the proj4 string for UTM Zone 18N is:
Code: Select all
+proj=utm +zone=18 +ellps=WGS84 +datum=WGS84 +units=m +no_defsUsing this, the gdalwarp command-line looks like this:
Code: Select all
gdalwarp -t_srs "+proj=utm +zone=18 +ellps=WGS84 +datum=WGS84 +units=m +no_defs" -wo "NUM_THREADS=ALL_CPUS" -wm 4096 -multi -of GTiff -co "COMPRESS=LZW" -co "BIGTIFF=YES" DEM_LiDAR_1ft_2010_Improved_NYC_int.tif DEM_LiDAR_1ft_2010_Improved_NYC_int_wgs84.tifYou could achieve the same using ArcGIS. I think the Project Raster tool is used for this. I did however find that ArcGIS is several magnitudes slower in the transformation.
Okay, so this takes care of the coordinates used in the file. However, this is a DEM, meaning we also have vertical units, i.e. the actual value at each pixel. These values are still in feet, as gdalwarp only projects along the x and y axes. All pixels need to be re-calculated from feet to metres. I finally found this answer after some googling: https://gis.stackexchange.com/a/259904
Form the previous post we can see that the format in the file defines a US Survey Foot as 0.3048006096012192 metres. Ergo, we have to multiply each pixel with 0.3048006096012192 using gdal_calc.py as suggested. Then we can use ArcGIS to change the metadata in the GeoTIFF file to say that we have a metric vertical units.
I ended up with the following command-line:
Code: Select all
gdal_calc.py -A DEM_LiDAR_1ft_2010_Improved_NYC_int_wgs84.tif --type=Int16 --co="COMPRESS=LZW" --co="BIGTIFF=YES" --calc="A*0.3048006096012192" --outfile=DEM_LiDAR_1ft_2010_Improved_NYC_int_wgs84_metres.tifBTW: As far as I understood ArcGIS would also be able to do the exact same operation using the Times operation in the Math Toolbox of the Spatial Analyst Tools. But you need an ArcGIS Pro license for that. So, I am not using that one.
Finally, TransDEM wants the DEM to actually specify the vertical units. Otherwise, it now fails with:
Code: Select all
VerticalUnitsGeoTagKey tag not found.Re: NYC Open Data 1 foot Digital Elevation Model
Okay, phew, after all that warping and re-calculating, let's see what we have now. This time, I'll also use the -approx_stats flag to gdalinfo to give an approximate statistics calculation (which only take a few seconds). Of course, you can run the command with -stats which will calculate accurate statistics, but which will read the entire file and take a while to process.
Okay, this looks promising. We see a UTM Zone 18N WGS84 projection and it also shows metre as the unit type in Band 1 (at the bottom). Before we had feet between 0ft and 411ft, now we have 0m to 125m. And sure enough, 411ft is 125.273m. YAIH!
Okay, let's load it into TranzDEM… And it fails:
Okay, now I am out of my depths!!! What do I do now? And more importantly, how do I fix this? Any suggestions, please?
Code: Select all
> gdalinfo -approx_stats DEM_LiDAR_1ft_2010_Improved_NYC_int_wgs84_metres.tif
Driver: GTiff/GeoTIFF
Files: DEM_LiDAR_1ft_2010_Improved_NYC_int_wgs84_metres.tif
DEM_LiDAR_1ft_2010_Improved_NYC_int_wgs84_metres.tif.aux.xml
Size is 159866, 157890
Coordinate System is:
PROJCS["WGS_1984_UTM_Zone_18N",
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",-75],
PARAMETER["scale_factor",0.9996],
PARAMETER["false_easting",500000],
PARAMETER["false_northing",0],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]],
AUTHORITY["EPSG","32618"]]
Origin = (561797.957299843896180,4530866.120325331576169)
Pixel Size = (0.304707187527485,-0.304707187527485)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
COMPRESSION=LZW
INTERLEAVE=BAND
Corner Coordinates:
Upper Left ( 561797.957, 4530866.120) ( 74d15'57.65"W, 40d55'35.69"N)
Lower Left ( 561797.957, 4482755.902) ( 74d16'14.73"W, 40d29'35.54"N)
Upper Right ( 610510.277, 4530866.120) ( 73d41'15.23"W, 40d55'17.21"N)
Lower Right ( 610510.277, 4482755.902) ( 73d41'45.75"W, 40d29'17.34"N)
Center ( 586154.117, 4506811.011) ( 73d58'48.34"W, 40d42'27.75"N)
Band 1 Block=159866x1 Type=Int16, ColorInterp=Gray
Min=0.000 Max=125.000
Minimum=0.000, Maximum=125.000, Mean=5.483, StdDev=11.902
NoData Value=-32767
Unit Type: metre
Metadata:
STATISTICS_APPROXIMATE=YES
STATISTICS_MAXIMUM=125
STATISTICS_MEAN=5.4827286885428
STATISTICS_MINIMUM=0
STATISTICS_STDDEV=11.901829498744Okay, let's load it into TranzDEM… And it fails:
Code: Select all
Errors reading file "DEM_LiDAR_1ft_2010_Improved_NYC_int_wgs84_metres.tif".Re: NYC Open Data 1 foot Digital Elevation Model
This is the first ever DEM I have encountered in 16 years all over the world that has an imperial grid
. And you are right, TransDEM does not like this. Conversion to UTM/WGS84 should be okay. For the vertical units, it should be metres and the format should be 16bit floating point. (Integer does not make much sense with hi-res LIDAR.)
Unfortunately, I am not able to look into the details at the moment as I am away on a holiday trip this week.
Unfortunately, I am not able to look into the details at the moment as I am away on a holiday trip this week.
Re: NYC Open Data 1 foot Digital Elevation Model
Okay I ran the same commands to warp and recalculate the floating-point data set, turns out GDAL deals as well with Imagine Images as with GeoTIFF.
I agree that it makes more sense that high-resolution LiDAR should use floats rather than ints. But when I load my new floating-point file (metric, WGS84/UTM 18N) in TranzDEM, it issues the same error as with the Integer-Raster:
Code: Select all
> gdalwarp -t_srs "+proj=utm +zone=18 +ellps=WGS84 +datum=WGS84 +units=m +no_defs" -wo "NUM_THREADS=ALL_CPUS" -wm 4096 -multi -of GTiff -ot Float32 -co "COMPRESS=LZW" -co "BIGTIFF=YES" DEM_LiDAR_1ft_2010_Improved_NYC.img DEM_LiDAR_1ft_2010_Improved_NYC_float_wgs84.tif
> gdal_calc.py -A DEM_LiDAR_1ft_2010_Improved_NYC_float_wgs84.tif --type=Float32 --co="COMPRESS=LZW" --co="BIGTIFF=YES" --calc="A*0.3048006096012192" --outfile=DEM_LiDAR_1ft_2010_Improved_NYC_float_wgs84_metres.tif
> gdalinfo -approx_stats DEM_LiDAR_1ft_2010_Improved_NYC_float_wgs84_metres.tif
Code: Select all
Driver: GTiff/GeoTIFF
Files: DEM_LiDAR_1ft_2010_Improved_NYC_float_wgs84_metres.tif
Size is 159866, 157890
Coordinate System is:
PROJCS["WGS 84 / UTM zone 18N",
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,
AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4326"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",-75],
PARAMETER["scale_factor",0.9996],
PARAMETER["false_easting",500000],
PARAMETER["false_northing",0],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]],
AXIS["Easting",EAST],
AXIS["Northing",NORTH],
AUTHORITY["EPSG","32618"]]
Origin = (561797.957299843896180,4530866.120325328782201)
Pixel Size = (0.304707187527488,-0.304707187527488)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
COMPRESSION=LZW
INTERLEAVE=BAND
Corner Coordinates:
Upper Left ( 561797.957, 4530866.120) ( 74d15'57.65"W, 40d55'35.69"N)
Lower Left ( 561797.957, 4482755.902) ( 74d16'14.73"W, 40d29'35.54"N)
Upper Right ( 610510.277, 4530866.120) ( 73d41'15.23"W, 40d55'17.21"N)
Lower Right ( 610510.277, 4482755.902) ( 73d41'45.75"W, 40d29'17.34"N)
Center ( 586154.117, 4506811.011) ( 73d58'48.34"W, 40d42'27.75"N)
Band 1 Block=159866x1 Type=Float32, ColorInterp=Gray
Minimum=-20.035, Maximum=124.841, Mean=5.432, StdDev=11.926
NoData Value=3.40282346600000016e+38
Metadata:
STATISTICS_APPROXIMATE=YES
STATISTICS_MAXIMUM=124.84114837646
STATISTICS_MEAN=5.4315597119415
STATISTICS_MINIMUM=-20.035221099854
STATISTICS_STDDEV=11.92583221638Code: Select all
Errors reading file "DEM_LiDAR_1ft_2010_Improved_NYC_float_wgs84_metres.tif".Re: NYC Open Data 1 foot Digital Elevation Model
I will have to make this conversion myself and analyse which of the many variables in TIFF/GeoTIFF triggers the error, but I need to be back home for this.
Re: NYC Open Data 1 foot Digital Elevation Model
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:
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.
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:
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: Select all
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)
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: Select all
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)
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.