Data preparation of land subsidence susceptibility mapping analysis in Xuzhou

Introduction

Data Preparing

DEM

An advanced space borne therman emission and reflection radiometer (ASTER) digital elevation model (DEM) was obtained (1 arc second or approximately 30 meter resolution) and processed in the GIS environment to produce topography-related factor layers (altitude, slope angle, slope aspect, plan curvature, profile curvature and TWI).

Download

ASTGTM v003 (ASTER Global Digital Elevation Model 1 arc second) can be accessed using different tools. In this study, we use NASA Earthdata Search tool to download DEM in Xuzhou.

Merging Multiple DEM in ArcGIS pro

Mosaic to New Raster → Select downloaded two DEMs N34E116_dem.tiff and N34E117_dem.tiff → Raster Dataset Name with Extension: elevation_xuzhou.dat → Pixel Type: 32-bot float → Number of Bands: 1 → Run

Export Data

Export DEM to features with clip option.

Surface Parameter

Search surface parameter in ArcGIS pro, set parameters to1:

Slope, aspect, curvature

• Quadratic
• Use adaptive neighborhood
• Degrees

OpenStreetMap

Download road map

Use overpass tool to download roadmap of a specific city, such as Xuhzou2.

1. Get city ID by overpass-api

2. Insert code in Overpass API Query Form

<osm-script>
  <query type="relation">
    <has-kv k="boundary" v="administrative"/>
	<has-kv k="name:zh" v="徐州市"/>
  </query>
  <print/></osm-script>

3. Download file interprept, open with Notepad, get city ID in

4. Insert code again in Overpass API Query Form

<osm-script timeout="1800" element-limit="100000000">
  <union>
    <area-query ref="3603218572"/>
    <recurse type="node-relation" into="rels"/>
    <recurse type="node-way"/>
    <recurse type="way-relation"/>
  </union>
  <union>
    <item/>
    <recurse type="way-node"/>
  </union>
  <print mode="body"/>
</osm-script>

5. Download interprepter file, this is the roadmap of Xuzhou city.

6. Open with QGIS and convert OSM type file to shp file.

7. Add OpenStreetMap layer in QGIS

8. Select interpreter_lines to load in QGIS

9. Right click on layer and select Export → Save feature as → ESRI shapefile → Save to path D:\Data\Roads\OSM\Roads_xuzhou_OSM.shp

Distance to Roads and Distance to stream

For preparing these two types of data, we use distance accumaluation tool in ArcGIS pro.

Distance to Roads

1. Search tool distance accumulation

2. Insert parameter as following fig

3. In environment tab, set Extent as specific layer xuzhou_center_dissolve

4. Run same process for study area Fengpei

Distance to stream

Same process as Distance to roads

TWI

1. Fill

2. Slope

3. Flow Direction

4. Flow Accumulation

5. Calculate SCA:

• Raster Calculator

"Land Subsidence Susceptibility Mapping\DEM\eleva_fp_fill"*900 / Con("Land Subsidence Susceptibility Mapping\Topography_related\flow_dir_xz" == 1,30,Con("Land Subsidence Susceptibility Mapping\Topography_related\flow_dir_xz"  == 4,30,Con("Land Subsidence Susceptibility Mapping\Topography_related\flow_dir_xz"  == 16,30,Con("Land Subsidence Susceptibility Mapping\Topography_related\flow_dir_xz"  == 64,30,Con("Land Subsidence Susceptibility Mapping\Topography_related\flow_dir_xz"  == 2,30 * SquareRoot(2),Con("Land Subsidence Susceptibility Mapping\Topography_related\flow_dir_xz"  == 8,30 * SquareRoot(2),Con("Land Subsidence Susceptibility Mapping\Topography_related\flow_dir_xz"  == 32,30 * SquareRoot(2),Con("Land Subsidence Susceptibility Mapping\Topography_related\flow_dir_xz"  == 128,30 * SquareRoot(2)))))))))

"Land Subsidence Susceptibility Mapping\Topography_related\flow_accum_fp"*900 / Con("Land Subsidence Susceptibility Mapping\Topography_related\flow_dir_fp.tif" == 1,30,Con("Land Subsidence Susceptibility Mapping\Topography_related\flow_dir_fp.tif" == 4,30,Con("Land Subsidence Susceptibility Mapping\Topography_related\flow_dir_fp.tif" == 16,30,Con("Land Subsidence Susceptibility Mapping\Topography_related\flow_dir_fp.tif" == 64,30,Con("Land Subsidence Susceptibility Mapping\Topography_related\flow_dir_fp.tif" == 2,30 * SquareRoot(2),Con("Land Subsidence Susceptibility Mapping\Topography_related\flow_dir_fp.tif" == 8,30 * SquareRoot(2),Con("Land Subsidence Susceptibility Mapping\Topography_related\flow_dir_fp.tif" == 32,30 * SquareRoot(2),Con("Land Subsidence Susceptibility Mapping\Topography_related\flow_dir_fp.tif" == 128,30 * SquareRoot(2)))))))))

6. Calculate TWI:

• Raster Calculator

Ln("Land Subsidence Susceptibility Mapping\Topography_related\SCA_xz" / Tan(Con("Land Subsidence Susceptibility Mapping\Topography_related\slope_xz_fill" <= 0,0.00001,Con("Land Subsidence Susceptibility Mapping\Topography_related\slope_xz_fill" > 0,"Land Subsidence Susceptibility Mapping\Topography_related\slope_xz_fill" * 3.1415926 / 180))))