-
Notifications
You must be signed in to change notification settings - Fork 0
/
ex4_import_data.py
executable file
·75 lines (51 loc) · 3.04 KB
/
ex4_import_data.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
#!/usr/bin/env python
import grass.script as gscript
# subprocess is a Python module that lets you execute processes in the command
# line from within your python script e.g. gdal_transform.
import subprocess
# os is a package for using operating system dependent functionality e.g. list files in a directory
import os
def main():
# Path to the folder containing all data to be imported
path_data = "C:/Users/nb152/Documents/OpenSourceGIS_exercise_4_data"
# 0. Import cities (survey capture areas)
path_cities = os.path.join(path_data, 'nz_cities', 'nz-survey-capture-areas-sca.shp')
gscript.run_command('v.in.ogr', input=path_cities, output='cities')
# 1. Importing railways (nz_railway.shp)
# ---------------------------------------
path_railways = os.path.join(path_data, 'nz_railway', 'nz_railway.shp')
gscript.run_command('v.in.ogr', input=path_railways, output='railways')
# 2. Importing rainfall data (avg_annual_rainfall.tif)
# ------------------------------------------------------
path_rainfall = os.path.join(path_data, 'avg_annual_rainfall.tif')
path_rainfall_corrected = os.path.join(path_data, 'avg_annual_rainfall_corrected.tif')
# Use gdal_translate to assign the correct crs to the avg_annual_rainfall.tif
subprocess.call(['gdal_translate', '-a_srs', 'EPSG:4326', path_rainfall, path_rainfall_corrected])
# Import the corrected avg_annual_rainfall_corrected.tif into GRASS
gscript.run_command('r.import', input=path_rainfall_corrected, output='rainfall')
# 3. Importing sealed roads (nz_road.shp)
# -----------------------------------------
path_roads = os.path.join(path_data, 'nz_road', 'nz_road.shp')
path_roads_reprojected = os.path.join(path_data, 'nz_road','nz_road_reprojected.shp')
# Repoject nz_roads.shp using OGR2OGR
subprocess.call(['ogr2ogr', '-f', 'ESRI Shapefile', '-t_srs', 'EPSG:32760', path_roads_reprojected, path_roads])
# Import the reprojected roads file into GRASS GIS
gscript.run_command('v.in.ogr', input=path_roads_reprojected, layer='nz_road_reprojected', output='sealed_roads', where="surface='sealed'")
# 4. Importing airports (nz_airports.csv)
# ------------------------------------------
path_airports = os.path.join(path_data, 'nz_airport.csv')
gscript.run_command('v.in.ascii', input=path_airports, output='airports', separator='comma', skip=1, x=4, y=5)
# 5. Set region to match rainfall dataset
# -----------------------------------------
gscript.run_command('g.region', rast='rainfall')
# Print some information about the mapset to check if everything worked
# Print projection of mapset
gscript.run_command('g.proj', flags='p')
# Print extent of region
gscript.run_command('g.region', flags='p')
# List all raster and vector layers
gscript.run_command('g.list', type='vector')
gscript.run_command('g.list', type='raster')
print('Importing data done.')
if __name__ == '__main__':
main()