Skip to content

Commit

Permalink
Add provider for SRTM 1-arc second files (#3)
Browse files Browse the repository at this point in the history 
* add srtm 1-arc second provider

* fix srtm 1-arc second provider for multiple rows per strip

* fix code style

* remove hhvm from travis build matrix

* add srtm 1-arc second test files
  • Loading branch information
@laufhannes
laufhannes committed Jun 9, 2017
1 parent 3574b98 commit acc7a48
Show file tree
Hide file tree
Showing 5 changed files with 270 additions and 11 deletions.
1 change: 0 additions & 1 deletion .travis.yml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,6 @@ php:
- 5.5
- 5.6
- 7.0
- hhvm

sudo: false

Expand Down
8 changes: 8 additions & 0 deletions bin/download-testfiles.sh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,4 +11,12 @@ do
unzip -o tests/testfiles/${file}.zip ${file}.tif -d tests/testfiles
rm -rf tests/testfiles/${file}.zip
fi
done

for file in n40_w075 n49_e007 n51_w001 s23_e017 s34_e151
do
if [ ! -f tests/testfiles/${file}_1arc_v3.tif ]
then
curl https://cdn.runalyze.com/static/srtm1/${file}_1arc_v3.tif > tests/testfiles/${file}_1arc_v3.tif
fi
done
43 changes: 33 additions & 10 deletions src/Provider/GeoTIFF/GeoTIFFReader.php
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,9 @@
use Runalyze\DEM\Exception\RuntimeException;
use Runalyze\DEM\Provider\AbstractResourceReader;

/**
* @see http://www.awaresystems.be/imaging/tiff/tifftags.html
*/
class GeoTIFFReader extends AbstractResourceReader
{
/**
Expand All @@ -22,13 +25,6 @@ class GeoTIFFReader extends AbstractResourceReader
*/
const LEN_OFFSET = 4;

/**
* The number of bytes containing each item of elevation data
* (= BitsPerSample tag value / 8).
* @var int
*/
const BYTES_PER_SAMPLE = 2;

/**
* Magic number located at bytes 2-3 which identifies a TIFF file.
* @var int
Expand All @@ -44,6 +40,15 @@ class GeoTIFFReader extends AbstractResourceReader
/** @var int */
const TIFF_CONST_IMAGE_LENGTH = 257;

/** @var int */
const TIFF_CONST_BITS_PER_SAMPLE = 258;

/** @var int */
const TIFF_CONST_SAMPLES_PER_PIXEL = 277;

/** @var int */
const TIFF_CONST_ROWS_PER_STRIP = 278;

/** @var int */
const TIFF_CONST_STRIPBYTECOUNTS = 279;

Expand All @@ -59,6 +64,15 @@ class GeoTIFFReader extends AbstractResourceReader
/** @var int */
const UNKNOWN = -32768;

/** @var int The number of bytes containing each item of elevation data */
protected $BytesPerSample = 2;

/** @var int The number of components per pixel */
protected $SamplesPerPixel = 1;

/** @var int The number of rows per strip */
protected $RowsPerStrip = 1;

/** @var int */
protected $NumDataRows;

Expand Down Expand Up @@ -130,6 +144,15 @@ protected function readIFDEntries()
case static::TIFF_CONST_IMAGE_LENGTH:
$this->NumDataRows = $constData['offset'];
break;
case static::TIFF_CONST_BITS_PER_SAMPLE:
$this->BytesPerSample = (int) ($constData['offset'] / 8);
break;
case static::TIFF_CONST_SAMPLES_PER_PIXEL:
$this->SamplesPerPixel = $constData['offset'];
break;
case static::TIFF_CONST_ROWS_PER_STRIP:
$this->RowsPerStrip = $constData['offset'];
break;
case static::TIFF_CONST_STRIPOFSETS:
$this->StripOffsets = $constData['offset'];
break;
Expand Down Expand Up @@ -165,13 +188,13 @@ public function getExactRowAndColFor($relativeLatitude, $relativeLongitude)
*/
public function getElevationFor($row, $col)
{
fseek($this->FileResource, $this->StripOffsets + ($row * static::LEN_OFFSET));
fseek($this->FileResource, $this->StripOffsets + ceil($row / $this->RowsPerStrip) * static::LEN_OFFSET);

$firstColumnData = unpack('Voffset', fread($this->FileResource, static::LEN_OFFSET));

fseek($this->FileResource, $firstColumnData['offset'] + $col * static::BYTES_PER_SAMPLE);
fseek($this->FileResource, $firstColumnData['offset'] + ($this->NumDataCols * ($row % $this->RowsPerStrip) + $col) * $this->BytesPerSample);

$elevation = unpack('velevation', fread($this->FileResource, static::BYTES_PER_SAMPLE))['elevation'];
$elevation = unpack('velevation', fread($this->FileResource, $this->BytesPerSample))['elevation'];

return ($elevation <= self::UNKNOWN || $elevation >= -self::UNKNOWN) ? false : $elevation;
}
Expand Down
68 changes: 68 additions & 0 deletions src/Provider/GeoTIFF/SRTM1ArcSecondProvider.php
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,68 @@
<?php

/*
* This file is part of the Runalyze DEM Reader.
*
* (c) RUNALYZE <mail@runalyze.com>
*
* This source file is subject to the MIT license that is bundled
* with this source code in the file LICENSE.
*/

namespace Runalyze\DEM\Provider\GeoTIFF;

use Runalyze\DEM\Exception\InvalidArgumentException;

class SRTM1ArcSecondProvider extends SRTM4Provider
{
/** @var float */
const DEGREES_PER_TILE = 1.0;

/** @var string */
protected $FilenameFormat = '%s%02d_%s%03d.tif';

/**
* @param float $latitude
* @param float $longitude
* @return string
* @throws InvalidArgumentException
*/
protected function getFilenameFor($latitude, $longitude)
{
$this->loadTileReferencesFor($latitude, $longitude);

return sprintf(
$this->FilenameFormat,
$this->CurrentTileVerticalReference < 0 ? 's' : 'n',
abs($this->CurrentTileVerticalReference),
$this->CurrentTileHorizontalReference < 0 ? 'w' : 'e',
abs($this->CurrentTileHorizontalReference)
);
}

/**
* @param float $latitude
* @param float $longitude
*/
protected function loadTileReferencesFor($latitude, $longitude)
{
$this->CurrentTileVerticalReference = floor($latitude / static::DEGREES_PER_TILE);
$this->CurrentTileHorizontalReference = floor($longitude / static::DEGREES_PER_TILE);

$this->CurrentTileLatitude = $this->CurrentTileVerticalReference;
$this->CurrentTileLongitude = $this->CurrentTileHorizontalReference;
}

/**
* @param float $latitude
* @param float $longitude
* @return float[] array(row, col)
*/
protected function getExactRowAndColFor($latitude, $longitude)
{
return $this->ResourceReader->getExactRowAndColFor(
abs($this->CurrentTileLatitude + 1 - $latitude) / static::DEGREES_PER_TILE,
abs($longitude - $this->CurrentTileLongitude) / static::DEGREES_PER_TILE
);
}
}
161 changes: 161 additions & 0 deletions src/Tests/Provider/GeoTIFF/SRTM1ArcSecondProviderTest.php
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,161 @@
<?php

/*
* This file is part of the Runalyze DEM Reader.
*
* (c) RUNALYZE <mail@runalyze.com>
*
* This source file is subject to the MIT license that is bundled
* with this source code in the file LICENSE.
*/

namespace Runalyze\DEM\Tests\Provider\GeoTIFF;

use Runalyze\DEM\Interpolation\BilinearInterpolation;
use Runalyze\DEM\Provider\GeoTIFF\SRTM1ArcSecondProvider;

class SRTM1ArcSecondProviderTest extends \PHPUnit_Framework_TestCase
{
/** @var string */
const PATH_TO_FILES = '/../../../../tests/testfiles';

/**
* @var SRTM1ArcSecondProvider
*/
protected $Provider;

public function setUp()
{
$this->Provider = new SRTM1ArcSecondProvider(__DIR__.self::PATH_TO_FILES);
$this->Provider->setFilenameFormat('%s%02d_%s%03d_1arc_v3.tif');
$this->Provider->setInterpolation(new BilinearInterpolation());
}

/**
* @param string $filename
*/
protected function checkFile($filename)
{
if (!$this->fileIsThere($filename)) {
$this->markTestSkipped('Required testfile "'.$filename.'" is not available.');
}
}

/**
* @param string $filename
* @return bool
*/
protected function fileIsThere($filename)
{
return file_exists(__DIR__.self::PATH_TO_FILES.'/'.$filename);
}

public function testTileBoundary()
{
$this->checkFile('n49_e007_1arc_v3.tif');

if (
!$this->fileIsThere('n49_e006_1arc_v3.tif') &&
!$this->fileIsThere('n49_e008_1arc_v3.tif') &&
!$this->fileIsThere('n48_e007_1arc_v3.tif') &&
!$this->fileIsThere('n50_e007_1arc_v3.tif')
) {
$this->assertTrue($this->Provider->hasDataFor([[49.00000, 7.00000]]));
$this->assertTrue($this->Provider->hasDataFor([[49.99999, 7.00000]]));
$this->assertTrue($this->Provider->hasDataFor([[49.99999, 7.99999]]));
$this->assertTrue($this->Provider->hasDataFor([[49.00000, 7.99999]]));

$this->assertFalse($this->Provider->hasDataFor([[49.00000, 6.99999]]));
$this->assertFalse($this->Provider->hasDataFor([[49.00000, 8.00000]]));
$this->assertFalse($this->Provider->hasDataFor([[48.99999, 7.99999]]));
} else {
$this->markTestSkipped('Can\'t check boundaries of tile as too many files are there.');
}
}

/**
* @group germany
*/
public function testSingleElevationInGermany()
{
$this->checkFile('n49_e007_1arc_v3.tif');

// SRTM4: 237
$this->assertEquals(239, $this->Provider->getElevation(49.444722, 7.768889));
}

public function testMultipleElevationsInSydney()
{
$this->checkFile('s34_e151_1arc_v3.tif');

// SRTM4: 4, 4, 3
$this->assertEquals(
[3, 3, 2],
$this->Provider->getElevations(
[-33.8706555, -33.8705667, -33.8704860],
[151.1486918, 151.1486337, 151.1485585]
)
);
}

/**
* @group london
*/
public function testLocationsInLondon()
{
$this->checkFile('n51_w001_1arc_v3.tif');

// SRTM4: 20, 19, 19
$this->assertEquals(
[18, 19, 18],
$this->Provider->getElevations(
[51.5073509, 51.5074509, 51.5075509],
[-0.1277583, -0.1278583, -0.1279583]
)
);
}

public function testLocationsInWindhoek()
{
$this->checkFile('s23_e017_1arc_v3.tif');

// SRTM4: 1668, 1671, 1672
$this->assertEquals(
[1665, 1672, 1671],
$this->Provider->getElevations(
[-22.5700, -22.5705, -22.5710],
[17.0836, 17.0841, 17.0846]
)
);
}

public function testNewYork()
{
$this->checkFile('n40_w075_1arc_v3.tif');

// SRTM4: 26, 29, 39
$this->assertEquals(
[28, 27, 46],
$this->Provider->getElevations(
[40.7127, 40.7132, 40.7137],
[-74.0059, -74.0064, -74.0069]
)
);
}

public function testNewYorkWithoutInterpolation()
{
$this->checkFile('n40_w075_1arc_v3.tif');

$this->Provider->removeInterpolation();

// SRTM4: 26, 32, 32
$this->assertEquals(
[28, 28, 48],
$this->Provider->getElevations(
[40.7127, 40.7132, 40.7137],
[-74.0059, -74.0064, -74.0069]
)
);
}
}

0 comments on commit acc7a48

Please sign in to comment.