TransDEM Forum

TransDEM News, Support, Hints and Resources
It is currently 28 Mar 2024 19:39

All times are UTC + 1 hour




Post new topic Reply to topic  [ 9 posts ] 
Author Message
PostPosted: 22 Oct 2018 22:14 
Offline

Joined: 23 Sep 2018 17:38
Posts: 15
Location: Tromsø, Norway
Hi,

I recently found out that New York City publishes a NYC 1foot Digital Elevation Model (https://data.cityofnewyork.us/City-Government/1-foot-Digital-Elevation-Model-DEM-/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.


Top
 Profile  
 
PostPosted: 22 Oct 2018 22:23 
Offline

Joined: 23 Sep 2018 17:38
Posts: 15
Location: Tromsø, Norway
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:
Code:
Error reading DEM!


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.)
Code:
> 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

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.


Top
 Profile  
 
PostPosted: 22 Oct 2018 22:32 
Offline

Joined: 23 Sep 2018 17:38
Posts: 15
Location: Tromsø, Norway
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-Government/1-foot-Digital-Elevation-Model-DEM-Integer-Raster/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:
Code:
ProjLinearUnitsGeoKey tag not Linear_Meter [Linear_Foot_US_Survey]


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:
Code:
> 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.

Top
 Profile  
 
PostPosted: 22 Oct 2018 22:48 
Offline

Joined: 23 Sep 2018 17:38
Posts: 15
Location: Tromsø, Norway
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/epsg/wgs-84-utm-zone-18n/ the proj4 string for UTM Zone 18N is:
Code:
+proj=utm +zone=18 +ellps=WGS84 +datum=WGS84 +units=m +no_defs

Note that it actually says: units=m (i.e. metric units)

Using this, the gdalwarp command-line looks like this:
Code:
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.tif

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:
Code:
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.tif

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:
Code:
VerticalUnitsGeoTagKey tag not found.

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.


Top
 Profile  
 
PostPosted: 23 Oct 2018 15:25 
Offline

Joined: 23 Sep 2018 17:38
Posts: 15
Location: Tromsø, Norway
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.
Code:
> 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.901829498744

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:
Code:
Errors reading file "DEM_LiDAR_1ft_2010_Improved_NYC_int_wgs84_metres.tif".

Okay, now I am out of my depths!!! What do I do now? And more importantly, how do I fix this? Any suggestions, please?


Top
 Profile  
 
PostPosted: 23 Oct 2018 17:16 
Offline

Joined: 05 Jan 2011 16:45
Posts: 1463
This is the first ever DEM I have encountered in 16 years all over the world that has an imperial grid :o . 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.


Top
 Profile  
 
PostPosted: 23 Oct 2018 19:46 
Offline

Joined: 23 Sep 2018 17:38
Posts: 15
Location: Tromsø, Norway
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.
Code:
> 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:
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.92583221638

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:
Errors reading file "DEM_LiDAR_1ft_2010_Improved_NYC_float_wgs84_metres.tif".


Top
 Profile  
 
PostPosted: 23 Oct 2018 22:56 
Offline

Joined: 05 Jan 2011 16:45
Posts: 1463
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.


Top
 Profile  
 
PostPosted: 25 Oct 2018 08:59 
Offline

Joined: 23 Sep 2018 17:38
Posts: 15
Location: Tromsø, Norway
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.


Top
 Profile  
 
Display posts from previous:  Sort by  
Post new topic Reply to topic  [ 9 posts ] 

All times are UTC + 1 hour


Who is online

Users browsing this forum: SemrushBot [Bot] and 5 guests


You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum

Search for:
Jump to:  
cron

Imprint & Privacy

Powered by phpBB® Forum Software © phpBB Group