-
Notifications
You must be signed in to change notification settings - Fork 0
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Merge branch 'master' of github.com:ansonl/DEM2STL
- Loading branch information
Showing
10 changed files
with
456 additions
and
26 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -1,15 +1,3 @@ | ||
| #buffer at 1k width | ||
| #rasterize 250x250m | ||
| #export in 102004 | ||
|
|
||
| # Upsample rasterized source hydro1k to 1500x1500m resolution | ||
| # source resolution 25x25m | ||
| # output resolution 750x750m | ||
|
|
||
| # 25x25 > 1000x1000 > 750x750 > 500x500 | ||
|
|
||
| #gdalwarp -overwrite -of GTiff -tr 1500.0 1500.0 ./sources/hydro1k_mask_25m_102004.tif ./hydrographic-masks/hydro1k_mask_1250m_102004.tif -r max -multi -co "TILED=YES" -co "COMPRESS=ZSTD" -co "PREDICTOR=2" -co "BIGTIFF=IF_SAFER" -co "NUM_THREADS=ALL_CPUS" -wo "NUM_THREADS=ALL_CPUS" -wo "USE_OPENCL=TRUE" | ||
|
|
||
| #gdalwarp -overwrite -of GTiff -tr 750.0 750.0 ./hydrographic-masks/hydro1k_mask_1250m_102004.tif ./hydrographic-masks/hydro1k_mask_750m_102004.tif -r max -multi -co "TILED=YES" -co "COMPRESS=ZSTD" -co "PREDICTOR=2" -co "BIGTIFF=IF_SAFER" -co "NUM_THREADS=ALL_CPUS" -wo "NUM_THREADS=ALL_CPUS" -wo "USE_OPENCL=TRUE" | ||
|
|
||
| #gdalwarp -overwrite -of GTiff -tr 250.0 250.0 ./hydrographic-masks/hydro1k_mask_750m_102004.tif ./hydrographic-masks/hydro1k_mask_250m_102004.tif -r bilinear -multi -co "TILED=YES" -co "COMPRESS=ZSTD" -co "PREDICTOR=2" -co "BIGTIFF=IF_SAFER" -co "NUM_THREADS=ALL_CPUS" -wo "NUM_THREADS=ALL_CPUS" -wo "USE_OPENCL=TRUE" | ||
| #buffer at 1000m width | ||
| #reproject to wgs84 | ||
| #rasterize at source DEM 0.0002777777777777777775 degrees resolution |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,62 @@ | ||
| [ | ||
| { | ||
| "name": "west coast pacific", | ||
| "extent": "-125.288780633,-116.880174030,32.284590810,49.269976148 [EPSG:4326]", | ||
| "start": "-123.69, 38.21", | ||
| "elevation": "0.5" | ||
| }, | ||
| { | ||
| "name": "northeast atlantic ME-NJ", | ||
| "extent": "-75.741066225,-65.524609202,36.961878233,45.454570902 [EPSG:4326]", | ||
| "start": "-69.9,42.62", | ||
| "elevation": "0.5" | ||
| }, | ||
| { | ||
| "name": "south atlantic DE-FL", | ||
| "extent": "-81.8,-75.293373615,24.167165057,40.774163098 [EPSG:4326]", | ||
| "start": "-78.19, 32.95", | ||
| "elevation": "0.5" | ||
| }, | ||
| { | ||
| "name": "south gulf coast FL-TX", | ||
| "extent": "-98.353977224,-81.368591886,24.083078991,31 [EPSG:4326]", | ||
| "start": "-83.81, 29.13", | ||
| "elevation": "0.5" | ||
| }, | ||
| { | ||
| "name": "LA Bon secour bay", | ||
| "extent": "-89.571805556,-89.105694444,30.288750000,30.531250000 [EPSG:4326]", | ||
| "start": "-89.30523, 30.3411", | ||
| "elevation": "0.5" | ||
| }, | ||
| { | ||
| "name": "FL Saint Johns River N", | ||
| "extent": "-81.758750000,-81.494027778,29.654861111,30.089305556 [EPSG:4326]", | ||
| "start": "-81.5568, 29.7321", | ||
| "elevation": "3.5" | ||
| }, | ||
| { | ||
| "name": "FL Saint Johns River S", | ||
| "extent": "-81.725147988,-81.587560973,29.388604094,29.659191890 [EPSG:4326]", | ||
| "start": "-81.621, 29.6323", | ||
| "elevation": "5.5" | ||
| }, | ||
| { | ||
| "name": "FL Indian River N", | ||
| "extent": "-80.886527778,-80.497638889,28.283472222,28.835972222 [EPSG:4326]", | ||
| "start": "-80.8053, 28.6989", | ||
| "elevation": "0.5" | ||
| }, | ||
| { | ||
| "name": "FL Banana River N", | ||
| "extent": "-80.886527778,-80.497638889,28.283472222,28.835972222 [EPSG:4326]", | ||
| "start": "-80.5988, 28.5107", | ||
| "elevation": "0.5" | ||
| }, | ||
| { | ||
| "name": "FL Banana River S", | ||
| "extent": "-80.886527778,-80.497638889,28.283472222,28.835972222 [EPSG:4326]", | ||
| "start": "-80.64126, 28.3807", | ||
| "elevation": "0.5" | ||
| } | ||
| ] |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,27 @@ | ||
| import re | ||
| def get_trailing_number(s): | ||
| m = re.search(r'\d+$', s) | ||
| return int(m.group()) if m else None | ||
|
|
||
| GTIFF_creation_options = { | ||
| 'TILES': 'YES', | ||
| 'COMPRESS': 'ZSTD', | ||
| 'PREDICTOR': '2', | ||
| 'BIGTIFF': 'IF_SAFER', | ||
| 'NUM_THREADS': 'ALL_CPUS' | ||
| } | ||
|
|
||
| GTIFF_write_options = { | ||
| 'NUM_THREADS': 'ALL_CPUS' | ||
| #USE_OPENCL: 'TRUE' | ||
| } | ||
|
|
||
| configuration_keywords = { | ||
| 'GDAL_VRT_ENABLE_PYTHON': 'YES' | ||
| } | ||
|
|
||
| GTIFF_creation_options_str = ' '.join([f'-co "{key}={value}"' for key, value in GTIFF_creation_options.items()]) | ||
| GTIFF_write_options_str = ' '.join([f'-wo "{key}={value}"' for key, value in GTIFF_creation_options.items()]) | ||
|
|
||
| # https://gdal.org/programs/raster_common_options.html#cmdoption-config | ||
| configuration_keywords_str = ' '.join([f'--config "{key} {value}"' for key, value in configuration_keywords.items()]) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,121 @@ | ||
| import os | ||
| import json | ||
|
|
||
| from constants import * | ||
|
|
||
| def generateCoastlineDEMCommands(sourceDEMFilename): | ||
|
|
||
| lakeClippedMasksPath = 'coastline/' | ||
| lakeMaskFilename = 'coastline_hydrographic_mask_merge.vrt' | ||
| lakeMaskInt16Filename = os.path.splitext(lakeMaskFilename)[0] + "_int16" + ".vrt" | ||
|
|
||
| #load lakes data from lake-locations.json | ||
| input_file = open ('../hydrographic-masks/lake-locations.json') | ||
| lakes = json.load(input_file) | ||
|
|
||
| for location in lakes: | ||
|
|
||
| featureName = location['feature'] | ||
| clippedDEMFilename = featureName.replace(' ', '-') + '.tif' | ||
|
|
||
| featureExtent = location['extent'] | ||
| #re | ||
| extentMatches = re.findall('([-\d.]+),([-\d.]+),([-\d.]+),([-\d.]+) \[([A-Z:0-9]*)\]', featureExtent) | ||
| if len(extentMatches) == 5: | ||
| extents = extentMatches[0:4] | ||
| extentsCRS = extentMatches[4] | ||
| else: | ||
| print(f'# Could not parse feature extent for {featureName}. Got {len(extentMatches)}/5 matches.') | ||
| continue | ||
|
|
||
| featureStart = location['start'] | ||
|
|
||
| featureElevation = location['elevation'] | ||
|
|
||
| clipRasterCmd = ' '.join(( | ||
| 'gdalwarp', | ||
| f'-te {extents[0]} {extents[1]} {extents[2]} {extents[3]}', | ||
| f'-te_srs {extentsCRS}', | ||
| f'{GTIFF_write_options_str}', | ||
| '-r cubicspline', | ||
| '-multi', | ||
| '-of GTiff', | ||
| f'{GTIFF_creation_options_str}', | ||
| '-overwrite', | ||
| f'{sourceDEMFilename}', | ||
| f'{clippedDEMFilename}' | ||
| )) | ||
|
|
||
| print(f'# Clip lake feature {featureName} - source DEM {sourceDEMFilename} → {clippedDEMFilename}') | ||
| print(clipRasterCmd) | ||
| print(f'# Run grass.lake with:\nStart: {featureStart}\nElevation: {featureElevation}\nLake raster location: {lakeClippedMasksPath}{clippedDEMFilename}') | ||
| print() | ||
|
|
||
| print(f'# Merge all lake rasters in {lakeClippedMasksPath}* → {lakeMaskFilename}') | ||
| print(f'gdalbuildvrt -resolution highest -overwrite {lakeMaskFilename} {lakeClippedMasksPath}/*.tif') | ||
| print() | ||
|
|
||
| print(f'# Set merged lake VRT data type from Float32 to Int16 to save space - {lakeMaskFilename} → {lakeMaskInt16Filename}') | ||
| print(f'gdal_translate -ot Int16 {lakeMaskFilename} {lakeMaskInt16Filename}') | ||
| print() | ||
|
|
||
| return lakeMaskInt16Filename | ||
|
|
||
| # all inputs must be WGS84 CRS | ||
| def generateHydrographicMaskFinalCommands(coastlinesFilename, streamsFilename, lakesFilename, targetSRS, modelResolution): | ||
|
|
||
| hydrographicMask_merge_filename = "north_america_hydrographic_mask_merge.vrt" | ||
|
|
||
| print(f'# Merge coastline, hydro1k, hydrolakes masks - {coastlinesFilename} {streamsFilename} {lakesFilename} → {hydrographicMask_merge_filename}') | ||
| print(f'gdalbuildvrt -resolution highest -overwrite {hydrographicMask_merge_filename} {coastlinesFilename} {streamsFilename} {lakesFilename}') | ||
| print() | ||
|
|
||
| hydrographicMask_reproject_filename = os.path.splitext(hydrographicMask_merge_filename)[0] + "_102004" + ".vrt" | ||
| hydrographicMask_reproject_resampling_method = "near" | ||
|
|
||
| hydrographicMask_reproject_cmd = ' '.join(( | ||
| 'gdalwarp', | ||
| f'-t_srs {targetSRS}', #ESRI:102004 | ||
| f'-r {hydrographicMask_reproject_resampling_method}', | ||
| '-multi', | ||
| '-overwrite', | ||
| f'{hydrographicMask_merge_filename}', | ||
| f'{hydrographicMask_reproject_filename}' | ||
| )) | ||
|
|
||
| print(f'# Reproject hydrographic mask merge to {targetSRS} CRS - {hydrographicMask_merge_filename} → {hydrographicMask_reproject_filename}') | ||
| print(hydrographicMask_reproject_cmd) | ||
| print() | ||
|
|
||
| hydrographicMask_reproject_250m_VRT_filename = os.path.splitext(hydrographicMask_reproject_filename)[0] + "_250m" + ".vrt" | ||
| hydrographicMask_reproject_250m_TIF_filename = os.path.splitext(hydrographicMask_reproject_250m_VRT_filename)[0] + ".tif" | ||
| hydrographicMask_reproject_250m_target_resolution = 250 | ||
| hydrographicMask_reproject_250m_resampling_method = "max" | ||
|
|
||
| hydrographicMask_reproject_250m_buildVRT_cmd = ' '.join(( | ||
| 'gdalwarp', | ||
| f'-tr {hydrographicMask_reproject_250m_target_resolution} {hydrographicMask_reproject_250m_target_resolution}', | ||
| f'-r {hydrographicMask_reproject_250m_resampling_method}', | ||
| '-multi', | ||
| '-overwrite', | ||
| f'{hydrographicMask_reproject_filename}', | ||
| f'{hydrographicMask_reproject_250m_VRT_filename}' | ||
| )) | ||
|
|
||
| hydrographicMask_reproject_250m_translate_cmd = ' '.join(( | ||
| 'gdal_translate', | ||
| f'-ot Int16', | ||
| f'{GTIFF_creation_options_str}', | ||
| f'{hydrographicMask_reproject_250m_VRT_filename}', | ||
| f'{hydrographicMask_reproject_250m_TIF_filename}' | ||
| )) | ||
|
|
||
| print(f'# Downscale hydrographic mask reproject to {hydrographicMask_reproject_250m_target_resolution}x{hydrographicMask_reproject_250m_target_resolution}m resolution. Create VRT first for previewing. - {hydrographicMask_reproject_filename} → {hydrographicMask_reproject_250m_VRT_filename}') | ||
| print(hydrographicMask_reproject_250m_buildVRT_cmd) | ||
| print() | ||
|
|
||
| print(f'# Create {hydrographicMask_reproject_250m_target_resolution}x{hydrographicMask_reproject_250m_target_resolution}m resolution TIF from VRT - {hydrographicMask_reproject_250m_VRT_filename} → {hydrographicMask_reproject_250m_TIF_filename}') | ||
| print(hydrographicMask_reproject_250m_translate_cmd) | ||
| print() | ||
|
|
||
| return hydrographicMask_reproject_250m_TIF_filename |
Oops, something went wrong.