Instantly share code, notes, and snippets.
Last active
September 1, 2017 23:10
-
Star
(0)
0
You must be signed in to star a gist -
Fork
(0)
0
You must be signed in to fork a gist
-
Save marioroy/61f72949568ac6bf2e98431ba172a55e to your computer and use it in GitHub Desktop.
Geo::GDAL + MCE::Shared + PDL parallel demonstration
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#!/usr/bin/env perl | |
# http://osgeo-org.1560.x6.nabble.com/gdal-dev-Experiments-with-multiprocessing-td5311397.html | |
# Demonstration - August 28, 2017 | |
# Use MCE::Shared 1.827 or later | |
# Usage: geo-shared-minmax.pl image.tif | |
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ | |
package main; | |
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ | |
use strict; | |
use warnings; | |
use feature 'state'; | |
use PDL; # must load PDL before MCE if passing piddles via IPC | |
use PDL::Types; | |
use MCE::Hobo; | |
use MCE::Mutex; | |
use MCE::Shared 1.827; | |
use Time::HiRes 'time'; | |
sub BufToPDL { | |
my ( $buf, $type, $xdim, $ydim ) = @_; | |
state $DATATYPE2PDL = { | |
$Geo::GDAL::Const::GDT_Byte => $PDL::Types::PDL_B, | |
$Geo::GDAL::Const::GDT_Int16 => $PDL::Types::PDL_S, | |
$Geo::GDAL::Const::GDT_UInt16 => $PDL::Types::PDL_US, | |
$Geo::GDAL::Const::GDT_Int32 => $PDL::Types::PDL_L, | |
$Geo::GDAL::Const::GDT_Float32 => $PDL::Types::PDL_F, | |
$Geo::GDAL::Const::GDT_Float64 => $PDL::Types::PDL_D | |
}; | |
my $pdl = PDL->new; | |
$pdl->set_datatype( $DATATYPE2PDL->{$type} // -1 ); | |
$pdl->setdims([ $xdim, $ydim ]); | |
my $data = $pdl->get_dataref(); | |
$$data = $buf; | |
$pdl->upd_data(); | |
return $pdl; | |
} | |
my ( $filename, $update ) = ( shift, 'ReadOnly' ); | |
print("Usage: $0 image.tif\n"), exit(1) unless $filename; | |
print "\nComputing: $filename\n"; | |
my $start = time; | |
my $iter = MCE::Shared->share( | |
{ module => 'Geo::GDAL::Iter' }, $filename, $update | |
); | |
my $mutex = MCE::Mutex->new(); | |
my $minv = MCE::Shared->scalar(255); | |
my $maxv = MCE::Shared->scalar(0); | |
sub task { | |
my ( $hobo_min, $hobo_max ) = ( 255, 0 ); | |
while ( my ($buf, $type, $xoff, $yoff, $xsize, $ysize) = $iter->Next ) { | |
my $pdl = BufToPDL($buf, $type, $xsize, $ysize); | |
my ( $min, $max ) = minmax $pdl; | |
$hobo_min = $min if ( $min < $hobo_min ); | |
$hobo_max = $max if ( $max > $hobo_max ); | |
} | |
$mutex->enter( sub { | |
$minv->set($hobo_min) if ( $hobo_min < $minv->get() ); | |
$maxv->set($hobo_max) if ( $hobo_max > $maxv->get() ); | |
}); | |
} | |
MCE::Hobo->create(\&task) for 1 .. MCE::Util::get_ncpu; | |
MCE::Hobo->wait_all; | |
$iter->Reset(); # resets iterator, optional | |
printf "Duration : %0.03f seconds\n", time - $start; | |
printf "Min value : %d\n", $minv->get(); | |
printf "Max value : %d\n", $maxv->get(); | |
exit 0; | |
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ | |
package Geo::GDAL::Iter; | |
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ | |
use strict; | |
use warnings; | |
use feature 'state'; | |
use Geo::GDAL; | |
use ProgressBar::Stack; | |
sub new { | |
my ( $class, $filename, $update ) = @_; | |
my %self; $self{band} = Geo::GDAL::Open($filename, $update)->Band(); | |
bless \%self, $class; | |
} | |
# $iter->Next() | |
sub Next { | |
my ( $self ) = @_; | |
state $min = sub { $_[0] < $_[1] ? $_[0] : $_[1] }; | |
state ($band, $type, $w, $h, $W, $H, $xoff, $yoff); | |
state $xfact = 4; state $yfact = 4; | |
if ( !defined $self->{done} ) { | |
$self->{done} = 0, init_progress(); | |
$band = $self->{band}; | |
$type = $band->{DataType}; | |
($w, $h) = $band->GetBlockSize; | |
($W, $H) = $band->Size; | |
($xoff, $yoff) = (0, 0); | |
} | |
return if $self->{done}; | |
if ( $xoff >= $W ) { | |
update_progress( sprintf "%0.1f", $yoff / $H * 100 ); | |
$xoff = 0, $yoff += ( $h * $yfact ); | |
if ( $yoff >= $H ) { | |
update_progress('100.0'), print "\n"; | |
$self->{done} = 1; | |
return; | |
} | |
} | |
my @args = ( $xoff, $yoff, $min->($W - $xoff, $w * $xfact), | |
$min->($H - $yoff, $h * $yfact) ); | |
my $buf = $band->ReadRaster( @args ); | |
$xoff += ( $w * $xfact ); | |
return ( $buf, $type, @args ); | |
} | |
# $iter->Reset() | |
sub Reset { | |
shift->{done} = undef; | |
return; | |
} | |
# $iter->Update( $pdl, $xoff, $yoff ) | |
sub Update { | |
shift->{band}->Piddle(@_); | |
return; | |
} | |
1; |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment