Use NASA SRTM data to estimate elevations

Project: Chankillo

Views: ¹⁴³

%md

# Chankillo Solar Observatory Worksheet \#2

Background: 

We will need the elevations of the observation points and the hill over which
sunrise will be observed.

## Get Elevation of Observation Points and Towers

Use NASA Shuttle Radar Topographic Mission (SRTM) online data to
obtain approximate elevation for the two observation points and the hill on which
the 13 towers rest.

Use longitude and latitude obtained by manually selecting points in Google Earth.

## bounding rectangle of the ruins:

- lat of northmost tower: -9.560356965861502
- lat of southmost tower: -9.562037172307594
- lon of west observation point: -78.2296179913435
- lon of east observation point: -78.22573325657659

## SRTMv4 data for 5-degree rectangle containing the ruins

### Download the SRTM file

Files are here: http://www.cgiar-csi.org/data/srtm-90m-digital-elevation-database-v4-1

Enter extremes of the bounding rectangle using the download GUI here http://srtm.csi.cgiar.org/SELECTION/inputCoord.asp

<img src="cgiar-select.png" width="80%">

Search returns one file. Download that.

<img src="cgiar-select2.png" width="80%">


### File format

Based on info at [ESRI ASCII Raster format](http://resources.esri.com/help/9.3/arcgisdesktop/com/gp_toolref/spatial_analyst_tools/esri_ascii_raster_format.htm),

Row 1 of the data is at the top of the raster, row 2 is just under row 1, and so on.

### File header
```
$ head -6 srtm_21_14.asc 
ncols         6001
nrows         6001
xllcorner     -80.000416424533
yllcorner     -10.000416206603
cellsize      0.00083333333333333
NODATA_value  -9999

```

### Data values represent elevation in meters

http://wiki.openstreetmap.org/wiki/SRTM

> Each file is a series of 16-bit integers giving the height of each cell in meters arranged from west to east and then north to south

Chankillo Solar Observatory Worksheet #2

Get Elevation of Observation Points and Towers

Use NASA Shuttle Radar Topographic Mission (SRTM) online data to obtain approximate elevation for the two observation points and the hill on which the 13 towers rest.

Use longitude and latitude obtained by manually selecting points in Google Earth.

bounding rectangle of the ruins:

lat of northmost tower: -9.560356965861502
lat of southmost tower: -9.562037172307594
lon of west observation point: -78.2296179913435
lon of east observation point: -78.22573325657659

SRTMv4 data for 5-degree rectangle containing the ruins

Download the SRTM file

Files are here: http://www.cgiar-csi.org/data/srtm-90m-digital-elevation-database-v4-1

Enter extremes of the bounding rectangle using the download GUI here http://srtm.csi.cgiar.org/SELECTION/inputCoord.asp

Search returns one file. Download that.

File format

Based on info at ESRI ASCII Raster format,

Row 1 of the data is at the top of the raster, row 2 is just under row 1, and so on.

File header

$ head -6 srtm_21_14.asc 
ncols         6001
nrows         6001
xllcorner     -80.000416424533
yllcorner     -10.000416206603
cellsize      0.00083333333333333
NODATA_value  -9999

Data values represent elevation in meters

http://wiki.openstreetmap.org/wiki/SRTM

Each file is a series of 16-bit integers giving the height of each cell in meters arranged from west to east and then north to south

# download files cover 5 degree range on both axes
# and have elevation for 6001 x 6001 grid
xllcorner = -80.000416424533
yllcorner = -10.000416206603
cellsize = 5.0/6000

# use x for ruins longitude, y for ruins latitude
xmax = -78.22573325657659
xmin = -78.2296179913435

ymax = -9.560356965861502
ymin = -9.562037172307594

# given longitude and latitude inside the SRTM rectangle,
# return 1-based row and column for that point

def get_row_and_col(pt):
    longitude, latitude = pt
    ncols = 6001
    nrows = 6001

    xllcorner = -80.000416424533
    yllcorner = -10.000416206603
    cellsize = 5.0/6000

    row_offset = (latitude - yllcorner) / cellsize
    row = nrows - int(row_offset) + 1
    
    col_offset = (longitude - xllcorner) / cellsize
    col = int(col_offset) + 1
    return row, col

# northmost tower
ptower1 = (-78.22748775951226,-9.560356965861502)
get_row_and_col(ptower1)
# allowing for 6-line header, we want line number 5480 from the SRTM file
# line 5480, col 2128 -> 204

(5474, 2128)

# southmost tower
ptower13 = (-78.22760136797875,-9.562037172307594)
get_row_and_col(ptower13)
# line 5482, col 2128 -> 209

(5476, 2128)

# west observation point
pwestop = (-78.2296179913435,-9.561116373529829)
get_row_and_col(pwestop)
# line 5481, col 2125 -> 175

(5475, 2125)

# east observation point
peastop = (-78.22573325657659,-9.561143170696472)
get_row_and_col(peastop)
# line 5481, col 2130 -> 175

(5475, 2130)

%sh
# print the elevations in meters for the bounding rectangle by extracting info
# from the SRTM file
# rows of matrix top to bottom represent north to south
# columns of matrix left to right represent west to east
unzip -p srtm_21_14.zip srtm_21_14.asc |
  sed -n 5480,5482p |
  cut -d' ' -f2125-2130 |
  tee ruins_elevations.txt

176 179 204 184 167
179 187 210 192 175
174 189 209 187 174

mean([204, 210, 209]).n()

207.666666666667

Approximate elevations in m above sea level:

W and E observation points: 175
Crest of hill (average, running North-South) with 13 towers: 208