Merge 8140cfa into 21f618f

edwardoughton · Mar 6, 2020 · d60958f · d60958f
2 parents 21f618f + 8140cfa
commit d60958f
Show file tree

Hide file tree

Showing 8 changed files with 270 additions and 163 deletions.
diff --git a/data/S004E026_AVE_DSM.tif b/data/S004E026_AVE_DSM.tif
diff --git a/data/S004E026_AVE_DSM.tif.source.txt b/data/S004E026_AVE_DSM.tif.source.txt
@@ -0,0 +1,63 @@
+DOI: https://doi.org/10.5069/G94M92HB
+OT Collection ID: OT.112016.4326.2
+Dataset Name: ALOS World 3D - 30m Ellipsoidal
+Short Name: AW3D30_E
+Collection Platform: Satellite Data
+
+----------------------------------------------------------
+
+Dataset Overview: The ALOS Global Digital Surface Model (AW3D30) is a global dataset generated from images collected using the Panchromatic Remote-sensing  Instrument for Stereo Mapping (PRISM) aboard the Advanced Land Observing Satellite (ALOS) from 2006 to 2011. As described by the Japan Aerospace Exploration Agency: The Japan Aerospace Exploration Agency (JAXA) releases the global digital surface model (DSM) dataset with a horizontal resolution of approx. 30-meter mesh (1 arcsec) free of charge. The dataset has been compiled with images acquired by the Advanced Land Observing Satellite "DAICHI" (ALOS). The dataset is published based on the DSM dataset (5-meter mesh version) of the "World 3D Topographic Data", which is the most precise global-scale elevation data at this time, and its elevation precision is also at a world-leading level as a 30-meter mesh version. This dataset is expected to be useful for scientific research, education, as well as the private service sector that uses geospatial information.    Version: May 2016: Global terrestrial region (within approx. 82 deg. of N/S latitudes) of Version 1 released (approx. 22,100 tiles)  Note: JAXA provides two versions of AW3D30 created from the original 5-meter mesh using different downsampling methods: average (provided here) and median (not available from OpenTopography).
+
+Dataset Acknowledgement: J. Takaku, T. Tadono, K. Tsutsui : Generation of High Resolution Global DSM from ALOS PRISM, The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, pp.243-248, Vol. XL-4, ISPRS TC IV Symposium, Suzhou, China, 2014. [PDF file]
+
+  T. Tadono, H. Ishida, F. Oda, S. Naito, K. Minakawa, H. Iwamoto : Precise Global DEM Generation By ALOS PRISM, ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, pp.71-76, Vol.II-4, 2014.  [PDF file]
+
+Dataset Keywords: global, gdem, ellipsoidal, sar
+
+Survey Date: 01/01/2006 - 01/01/2011
+
+----------------------------------------------------------
+
+Data Provider and Roles:
+	1. Japan Aerospace Exploration Agency
+	Role: Collector
+	Url: https://global.jaxa.jp/
+
+
+----------------------------------------------------------
+
+Total lidar Points: N/A
+
+Area: N/A
+
+Point Density: N/A
+
+Raster Resolution: 30 meter
+
+Coordinates System:
+	Horizontal: WGS 1984 [EPSG: 4326]
+	Vertical: WGS84 (Ellipsoid)
+
+Survey Report: http://www.eorc.jaxa.jp/ALOS/en/aw3d30/
+
+----------------------------------------------------------
+
+Download and Access Products:
+	1. Raster
+		opentopoID: OTALOS.082017.4326.1
+		Data Processing: http://opentopo.sdsc.edu/raster?opentopoID=OTALOS.082017.4326.1
+
+	2. Raster
+		opentopoID: OTALOS.112016.4326.2
+		Data Processing: http://opentopo.sdsc.edu/raster?opentopoID=OTALOS.112016.4326.2
+
+
+----------------------------------------------------------
+
+Dataset Extent in KMZ format: https://cloud.sdsc.edu/v1/AUTH_opentopography/www/lidarkmls/raster/AW3D30_bounds.kmz
+
+Dataset Spatial Bounds:
+	North: 83.0°
+	South: -81.9999999999999°
+	East: 180.0°
+	West: -180.0°
diff --git a/data/S004E027_AVE_DSM.tif b/data/S004E027_AVE_DSM.tif
diff --git a/data/S_AVE_DSM.vrt b/data/S_AVE_DSM.vrt
@@ -0,0 +1,21 @@
+<VRTDataset rasterXSize="7200" rasterYSize="3600">
+  <SRS dataAxisToSRSAxisMapping="2,1">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"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]</SRS>
+  <GeoTransform>  2.6000000000000000e+01,  2.7777777777777778e-04,  0.0000000000000000e+00, -3.0000000000000000e+00,  0.0000000000000000e+00, -2.7777777777777778e-04</GeoTransform>
+  <VRTRasterBand dataType="Int16" band="1">
+    <ColorInterp>Gray</ColorInterp>
+    <SimpleSource>
+      <SourceFilename relativeToVRT="1">S004E026_AVE_DSM.tif</SourceFilename>
+      <SourceBand>1</SourceBand>
+      <SourceProperties RasterXSize="3600" RasterYSize="3600" DataType="Int16" BlockXSize="3600" BlockYSize="1" />
+      <SrcRect xOff="0" yOff="0" xSize="3600" ySize="3600" />
+      <DstRect xOff="0" yOff="0" xSize="3600" ySize="3600" />
+    </SimpleSource>
+    <SimpleSource>
+      <SourceFilename relativeToVRT="1">S004E027_AVE_DSM.tif</SourceFilename>
+      <SourceBand>1</SourceBand>
+      <SourceProperties RasterXSize="3600" RasterYSize="3600" DataType="Int16" BlockXSize="3600" BlockYSize="1" />
+      <SrcRect xOff="0" yOff="0" xSize="3600" ySize="3600" />
+      <DstRect xOff="3600" yOff="0" xSize="3600" ySize="3600" />
+    </SimpleSource>
+  </VRTRasterBand>
+</VRTDataset>
diff --git a/scripts/area.py b/scripts/area.py
@@ -202,28 +202,22 @@ def csv_writer(data, directory, filename):
     dem_path = BASE_PATH
     directory = DATA_PROCESSED
     directory_shapes = os.path.join(DATA_PROCESSED, 'shapes')
-
-    old_crs = 'EPSG:4326'
-    new_crs = 'EPSG:3857'
     cell_range = 20000
 
-    #create new geojson for Crystal Palace radio transmitter
-    transmitter = {
-        'type': 'Feature',
-        'geometry': {
-            'type': 'Point',
-            'coordinates': (-0.07491679518573545, 51.42413477117786)
-            },
-        'properties': {
-            'id': 'Crystal Palace radio transmitter'
-            }
-        }
-
     #Terrain Irregularity Parameter delta h (in meters)
-    tip = terrain_area(dem_path, transmitter, cell_range, old_crs)
+    tip = terrain_area(
+        os.path.join(dem_path, 'ASTGTM2_N51W001_dem.tif'),
+        -0.074916,
+        51.424134,
+        cell_range)
+    print('TIP for AST DEM', tip)
 
     print('Running itmlogic')
     output = itmlogic_area(tip)
 
-    print('Writing results to .csv')
+    print('Writing results to ', os.path.join(directory, 'uarea_output.csv'))
     csv_writer(output, directory, 'uarea_output.csv')
+
+    tip = terrain_area(os.path.join(dem_path, 'S_AVE_DSM.vrt'), 26.9976, -3.5409, cell_range)
+    print('TIP for ALOS DEM', tip)
+    output = itmlogic_area(tip)
diff --git a/scripts/p2p.py b/scripts/p2p.py
@@ -17,14 +17,15 @@
 from collections import OrderedDict
 
 import fiona
+from fiona.crs import from_epsg
+from pyproj import Transformer
 from shapely.geometry import LineString, mapping
 from shapely.ops import transform
 
 from itmlogic.qerfi import qerfi
 from itmlogic.qlrpfl import qlrpfl
 from itmlogic.avar import avar
 from terrain_module import terrain_p2p
-from pyproj import Transformer
 
 # #set up file paths
 CONFIG = configparser.ConfigParser()
@@ -260,6 +261,28 @@ def write_shapefile(data, directory, filename, crs):
             sink.write(datum)
 
 
+def straight_line_from_points(a, b):
+    return {
+        'type': 'Feature',
+        'geometry': {
+            'type': 'LineString',
+            'coordinates': [
+                (
+                    a['geometry']['coordinates'][0],
+                    a['geometry']['coordinates'][1]
+                ),
+                (
+                    b['geometry']['coordinates'][0],
+                    b['geometry']['coordinates'][1]
+                ),
+            ]
+        },
+        'properties': {
+            'id': 'terrain path'
+        }
+    }
+
+
 if __name__ == '__main__':
 
     dem_folder = os.path.join(BASE_PATH)
@@ -284,69 +307,46 @@ def write_shapefile(data, directory, filename, crs):
         116, 107, 104, 101,  98,  95, 103,  91,  97, 102, 107, 107, 107, 103,  98,
         94,  91, 105, 122, 122, 122, 122, 122, 137, 137, 137, 137, 137, 137, 137,
         137, 140, 144, 147, 150, 152, 159
-        ]
+    ]
 
     #create new geojson for Crystal Palace radio transmitter
     transmitter = {
         'type': 'Feature',
         'geometry': {
             'type': 'Point',
             'coordinates': (-0.07491679518573545, 51.42413477117786)
-            },
+        },
         'properties': {
             'id': 'Crystal Palace radio transmitter'
-            }
         }
+    }
 
     #create new geojson for Mursley
     receiver = {
         'type': 'Feature',
         'geometry': {
             'type': 'Point',
             'coordinates': (-0.8119433954872186, 51.94972494521946)
-            },
+        },
         'properties': {
             'id': 'Mursley'
-            }
         }
+    }
 
     #create new geojson for terrain path
-    line = {
-        'type': 'Feature',
-        'geometry': {
-            'type': 'LineString',
-            'coordinates': [
-                    (
-                        transmitter['geometry']['coordinates'][0],
-                        transmitter['geometry']['coordinates'][1]
-                    ),
-                    (
-                        receiver['geometry']['coordinates'][0],
-                        receiver['geometry']['coordinates'][1]
-                    ),
-                ]
-            },
-        'properties': {
-            'id': 'terrain path'
-            }
-        }
-
-    current_crs = 'EPSG:4326'
+    line = straight_line_from_points(transmitter, receiver)
 
     #run terrain module
     measured_terrain_profile, distance_km, points = terrain_p2p(
-        dem_folder, line, current_crs
-        )
-    print('Distance is {}'.format(distance_km))
+        os.path.join(dem_folder, 'ASTGTM2_N51W001_dem.tif'), line)
+    print('Distance is {}km'.format(distance_km))
 
     #check (out of interest) how many measurements are in each profile
     print('len(measured_terrain_profile) {}'.format(len(measured_terrain_profile)))
     print('len(original_surface_profile_m) {}'.format(len(original_surface_profile_m)))
 
     #run model and get output
-    output = itmlogic_p2p(
-        original_surface_profile_m, original_distance
-        )
+    output = itmlogic_p2p(original_surface_profile_m, original_distance)
 
     #grab coordinates for transmitter and receiver for writing to .csv
     transmitter_x = transmitter['geometry']['coordinates'][0]
@@ -368,3 +368,52 @@ def write_shapefile(data, directory, filename, crs):
     write_shapefile(points, directory_shapes, 'points.shp', new_crs)
 
     print('Completed run')
+
+
+    print('Start two-tile test')
+    transmitter = {
+        'type': 'Feature',
+        'geometry': {
+            'type': 'Point',
+            'coordinates': (26.676,-3.513)
+        },
+        'properties': {
+            'id': 'A'
+        }
+    }
+
+    receiver = {
+        'type': 'Feature',
+        'geometry': {
+            'type': 'Point',
+            'coordinates': (27.594,-3.514)
+        },
+        'properties': {
+            'id': 'B'
+        }
+    }
+
+    #create new geojson for terrain path
+    line = straight_line_from_points(transmitter, receiver)
+
+    #run terrain module
+    measured_terrain_profile, distance_km, points = terrain_p2p(
+        os.path.join(dem_folder, 'S_AVE_DSM.vrt'), line)
+
+    print("Profile [", measured_terrain_profile[0], ",  ... ,", measured_terrain_profile[-1], "]")
+    print("Distance is {}km".format(distance_km))
+    print("Sampled", len(points), "points")
+
+    schema = {
+        'geometry': 'Point',
+        'properties': {'elevation': 'float'}
+    }
+    crs = from_epsg(4326)
+    with fiona.open('data/processed/shapes/two_tile_points.shp',
+                    'w',
+                    driver='ESRI Shapefile',
+                    crs=crs,
+                    schema=schema) as fh:
+        for point in points:
+            fh.write(point)
+    print("Wrote two-tile profile to data/processed/shapes/two_tile_points.shp")