From ee4f0ef3d0f3422e03c6296641e85a2ec6344037 Mon Sep 17 00:00:00 2001 From: Shawn Laffan Date: Wed, 7 Aug 2019 17:10:46 +1000 Subject: [PATCH] Support export of cluster outputs to shapefile Need to check code works for additional lists, and add some proper tests (that needs thought). Updates #161 --- lib/Biodiverse/Cluster.pm | 232 +++++++++++++++++++++++++++++++++++++- 1 file changed, 231 insertions(+), 1 deletion(-) diff --git a/lib/Biodiverse/Cluster.pm b/lib/Biodiverse/Cluster.pm index 5a16eb679..232bcab65 100644 --- a/lib/Biodiverse/Cluster.pm +++ b/lib/Biodiverse/Cluster.pm @@ -3,6 +3,8 @@ package Biodiverse::Cluster; use 5.010; +our $VERSION = '2.99_004'; + use Carp; use strict; use warnings; @@ -13,9 +15,11 @@ use Scalar::Util qw/blessed/; use Time::HiRes qw /gettimeofday tv_interval time/; use List::Util qw /first reduce min max/; use List::MoreUtils qw /any natatime/; +use Scalar::Util qw /looks_like_number/; use Time::HiRes qw /time/; +use Sort::Key::Natural qw /natsort/; -our $VERSION = '2.99_004'; +use Geo::GDAL::FFI; use Biodiverse::Matrix; use Biodiverse::Matrix::LowMem; @@ -24,6 +28,7 @@ use Biodiverse::SpatialConditions; use Biodiverse::Progress; use Biodiverse::Indices; use Biodiverse::Exception; +use Biodiverse::Utilities qw/sort_list_with_tree_names_aa/; use Ref::Util qw { :all }; @@ -306,6 +311,231 @@ sub export_matrices { return; } + +my $shape_export_comment_text = <<'END_OF_SHAPE_COMMENT' +Note: If you export a list then each shape (point or polygon) +will be repeated for each list item. + +Choose the __no_list__ option to not do this, +in which case to attach any lists you will need to run a second +export to the delimited text format and then join them. +This is needed because shapefile field names can only be +11 characters long and cannot contain non-alphanumeric characters. + +Note also that shapefiles do not have an undefined value +so any undefined values will be converted to zeroes. + +Export of array lists to shapefiles is not supported. +END_OF_SHAPE_COMMENT + ; + +sub get_metadata_export_shapefile { + my $self = shift; + # get the available lists + my @lists = $self->get_lists_for_export (no_array_lists => 1); + unshift @lists, '__no_list__'; + + # nodata won't have much effect until we make the outputs symmetric + #my @nodata_meta = $self->get_nodata_export_metadata; + + my @parameters = ( + { # GUI supports just one of these + name => 'file', + type => 'file' + }, + { + name => 'list', + label_text => 'List to export', + type => 'choice', + choices => \@lists, + default => 0, + }, + { + name => 'shapetype', + label_text => 'Shape type', + type => 'choice', + choices => [qw /POLYGON POINT/], + default => 0, + }, + #@nodata_meta, + { + type => 'comment', + label_text => $shape_export_comment_text, + }, + ); + foreach (@parameters) { + bless $_, $parameter_metadata_class; + } + + my %args = ( + format => 'Shapefile', + parameters => \@parameters, + ); + + return wantarray ? %args : \%args; +} + +sub export_shapefile { + my $self = shift; + my %args = (nodata_value => -2**128, @_); + + $args{file} =~ s/\.shp$//i; + my $file = $args{file}; + + my $list_name = $args{list}; + if (defined $list_name && $list_name eq '__no_list__') { + $list_name = undef; + } + + my $nodata = $args{nodata_value}; + if (!looks_like_number $nodata) { + $nodata = -1 * 2**128; + } + + # we are writing as 2D or 3D points or polygons, + # only Point, PointZ or Polygon are used + my $shape_type = uc ($args{shapetype} // 'POLYGON'); + croak "Invalid shapetype for shapefile export\n" + if $shape_type ne 'POINT' and $shape_type ne 'POLYGON'; + + say "Exporting to shapefile $file"; + + my $def_query = $args{def_query}; + + my $bd = $self->get_basedata_ref; + my $gp = $bd->get_groups_ref; + my @cell_sizes = $bd->get_cell_sizes; # get a copy + my @axes_to_use = (0, 1); + my $use_z; + if (scalar @cell_sizes > 2) { + @axes_to_use = (0, 1, 2); # we use Z in this case + $use_z = 1; + } + + my $half_csizes = []; + foreach my $size (@cell_sizes[@axes_to_use]) { + push @$half_csizes, $size > 0 ? $size / 2 : 0.5; + } + + my @list_col_specs_gdal; + if (defined $list_name) { # repeated polys per list item + push @list_col_specs_gdal, + {Name => 'KEY', Type => 'String'}, + {Name => 'VALUE', Type => 'Real'}, + } + + my $layer = Geo::GDAL::FFI::GetDriver('ESRI Shapefile') + ->Create($file . '.shp') + ->CreateLayer({ + Name => 'export', + GeometryType => ucfirst (lc $shape_type), + Fields => [ + { + Name => 'BRANCH', + Type => 'String' + }, + @list_col_specs_gdal, + ], + }); + + # need the list of branches + my $branches = $self->get_node_hash; + + NODE: + foreach my $branch_name (sort_list_with_tree_names_aa ([keys %$branches])) { + my $branch = $branches->{$branch_name}; + my $terminals = $branch->get_terminal_elements; + + my $wkt; + my $z_wkt = $use_z ? ' Z' : ''; + if ($shape_type eq 'POLYGON') { + $wkt = "MULTIPOLYGON$z_wkt ("; + } + else { + $wkt = "MULTIPOINT$z_wkt ("; + } + + # add a shape for each of the terminals + foreach my $element (natsort keys %$terminals) { + + my $coord_axes = $gp->get_element_name_coord (element => $element); + #my $name_axes = $gp->get_element_name_as_array (element => $element); + + if ($shape_type eq 'POLYGON') { + my ($x, $y) = @{$coord_axes}[@axes_to_use]; + my ($width, $height) = @{$half_csizes}[@axes_to_use]; + my $min_x = $x - $width; + my $max_x = $x + $width; + my $min_y = $y - $height; + my $max_y = $y + $height; + my $z = $use_z + ? $coord_axes->[$axes_to_use[2]] + : ''; + + $wkt .= "((" + . "$min_x $min_y $z, " + . "$min_x $max_y $z, " + . "$max_x $max_y $z, " + . "$max_x $min_y $z, " + . "$min_x $min_y $z" + . ')), '; + } + elsif ($shape_type eq 'POINT') { + my @pt = @{$coord_axes}[@axes_to_use]; + $wkt .= "(" + . join (' ', @pt) + . "), "; + } + } + + # close off the geometry + $wkt =~ s/, $//; + $wkt .= ")"; + + my $branch_name = $branch->get_name; + + my @data_for_gdal_layer; + if ($list_name) { + my %list_data = $branch->get_list_values ( + #element => $branch_name, + list => $list_name, + ); + + # write a separate shape for each label + foreach my $key (natsort keys %list_data) { + my %data = ( + BRANCH => $branch_name, + KEY => $key, + VALUE => ($list_data{$key} // $nodata), + ); + push @data_for_gdal_layer, \%data; + } + } + else { + my %data = ( + BRANCH => $branch_name, + ); + push @data_for_gdal_layer, \%data; + } + + foreach my $data_hr (@data_for_gdal_layer) { + my $f = Geo::GDAL::FFI::Feature->new($layer->GetDefn); + foreach my $key (keys %$data_hr) { + $f->SetField(uc ($key) => $data_hr->{$key}); + } + my $g = Geo::GDAL::FFI::Geometry->new(WKT => $wkt); + $f->SetGeomField($g); + $layer->CreateFeature($f); + } + } + + # close off + $layer = undef; + + return; +} + + # some of this can be refactored wth Spatial::get_spatial_conditions_ref sub process_spatial_conditions_and_def_query { my $self = shift;