From 9429f8a8bd7a0e9ea406d2577e127ec9b6c932a4 Mon Sep 17 00:00:00 2001 From: Ed J Date: Wed, 25 Sep 2024 23:20:55 +0000 Subject: [PATCH] deprecate t_raster2fits, add raster2fits --- Basic/Core/Core.xs | 45 ++++++------- Changes | 1 + Graphics/TriD/DemoTriD1.pm | 2 +- Libtmp/Transform/Cartography/Cartography.pm | 73 +++++++++++++-------- Libtmp/Transform/Proj4/t/proj_transform.t | 2 +- Libtmp/Transform/t/transform.t | 4 +- 6 files changed, 70 insertions(+), 57 deletions(-) diff --git a/Basic/Core/Core.xs b/Basic/Core/Core.xs index bb9cefce3..6bcfb861d 100644 --- a/Basic/Core/Core.xs +++ b/Basic/Core/Core.xs @@ -1232,34 +1232,29 @@ sethdr(p,h) SV * hdr(p) pdl *p - CODE: - pdl_barf_if_error(pdl_make_physdims(p)); - - /* Make sure that in the undef case we return not */ - /* undef but an empty hash ref. */ - - if((p->hdrsv==NULL) || (p->hdrsv == &PL_sv_undef)) { - p->hdrsv = (void*) newRV_noinc( (SV*)newHV() ); - } - - RETVAL = newRV( (SV*) SvRV((SV*)p->hdrsv) ); - - OUTPUT: - RETVAL +CODE: + pdl_barf_if_error(pdl_make_physdims(p)); + /* Make sure that in the undef case we return not */ + /* undef but an empty hash ref. */ + if((p->hdrsv==NULL) || (p->hdrsv == &PL_sv_undef)) { + p->hdrsv = (void*) newRV_noinc( (SV*)newHV() ); + } + RETVAL = newRV( (SV*) SvRV((SV*)p->hdrsv) ); +OUTPUT: + RETVAL SV * gethdr(p) - pdl *p - CODE: - pdl_barf_if_error(pdl_make_physdims(p)); - - if((p->hdrsv==NULL) || (p->hdrsv == &PL_sv_undef)) { - RETVAL = &PL_sv_undef; - } else { - RETVAL = newRV( (SV*) SvRV((SV*)p->hdrsv) ); - } - OUTPUT: - RETVAL + pdl *p +CODE: + pdl_barf_if_error(pdl_make_physdims(p)); + if((p->hdrsv==NULL) || (p->hdrsv == &PL_sv_undef)) { + RETVAL = &PL_sv_undef; + } else { + RETVAL = newRV( (SV*) SvRV((SV*)p->hdrsv) ); + } +OUTPUT: + RETVAL SV * unpdl(x) diff --git a/Changes b/Changes index f9c67eb5b..fe0f69c24 100644 --- a/Changes +++ b/Changes @@ -5,6 +5,7 @@ - add Transform::Cartography::clean_lines options: orange, filter_nan - Transform::Proj4 captures projection information at build not runtime - move Stats::diff to Ufunc::diffover +- deprecate t_raster2fits, add raster2fits 2.092 2024-09-07 - add Type::howbig diff --git a/Graphics/TriD/DemoTriD1.pm b/Graphics/TriD/DemoTriD1.pm index fefa1f831..1cd3b56df 100644 --- a/Graphics/TriD/DemoTriD1.pm +++ b/Graphics/TriD/DemoTriD1.pm @@ -233,7 +233,7 @@ my @demo = ( use PDL::Transform::Cartography; eval { # this is in case no NetPBM, i.e. can't load Earth images $shape = earth_shape(); - $floats = t_raster2float()->apply($shape->mv(2,0)); + $floats = t_raster2float()->apply($shape->dummy(0,3)); $radius = $floats->slice('(2)'); # r g b all same $radius *= float((6377.09863 - 6370.69873) / 6371); $radius += float(6370.69873 / 6371); diff --git a/Libtmp/Transform/Cartography/Cartography.pm b/Libtmp/Transform/Cartography/Cartography.pm index 9d3d11849..6e1858477 100644 --- a/Libtmp/Transform/Cartography/Cartography.pm +++ b/Libtmp/Transform/Cartography/Cartography.pm @@ -258,6 +258,7 @@ our $VERSION = "0.6"; $VERSION = eval $VERSION; our @EXPORT_OK = qw( graticule earth_image earth_coast earth_shape clean_lines t_unit_sphere + raster2fits t_orthographic t_rot_sphere t_caree t_carree t_mercator t_utm t_sin_lat t_sinusoidal t_conic t_albers t_lambert t_stereographic t_gnomonic t_az_eqd t_az_eqa t_vertical t_perspective t_hammer t_aitoff @@ -266,6 +267,9 @@ our @EXPORT_OK = qw( our @EXPORT = @EXPORT_OK; our %EXPORT_TAGS = (Func=>\@EXPORT_OK); +our @PLATE_CARREE = ([qw(Longitude Latitude RGB)], + [qw(degrees degrees index)], [[-180,180], [-90,90], [0,2]]); + ############################## # Steal _opt from PDL::Transform. *PDL::Transform::Cartography::_opt = \&PDL::Transform::_opt; @@ -472,7 +476,7 @@ sub earth_image { barf("earth_image: $f not found in \@INC\n") if !$found; barf("earth_image: couldn't load $f; you may need to install netpbm.\n") unless defined($im); - t_raster2fits()->apply($im); + raster2fits($im, @PLATE_CARREE); } =head2 earth_shape @@ -535,8 +539,42 @@ sub earth_shape { unless defined($found); barf("earth_shape: couldn't load $f; you may need to install netpbm.\n") unless defined($im); - $im = $im->dummy(0,3); # fake RGB - t_raster2fits()->apply($im); + raster2fits($im, @PLATE_CARREE); +} + +=head2 raster2fits + +=for usage + + $pdl_fits = raster2fits($pdl, \@axislabels, \@axisunits, \@axisranges); + $pdl_fits = raster2fits($pdl, @PDL::Transform::Cartography::PLATE_CARREE); + +=for ref + +Convert a raster ([3,]x,y) to FITS (x,y[,3]), with a suitable header +for the given parameters. + +=cut + +sub raster2fits { + die "Usage: raster2fits(\$d, \\\@axislabels, \\\@axisunits, \\\@axisranges)\n" if @_ != 4; + my ($d, $axislabels, $axisunits, $axisranges) = @_; + my $is_single_plane = (my $ndims = $d->ndims) <= 2; + my $out = $is_single_plane ? $d : $d->mv(0,2); + my $h = $out->fhdr; + $h->{SIMPLE} = 'T'; + $h->{NAXIS} = $ndims; + local $_; + my @dims = $out->dims; + $h->{"NAXIS".($_+1)} = $dims[$_] for 0..$ndims-1; + $h->{"CRVAL".($_+1)} = 0 for 0..$ndims-1; + $h->{"CRPIX".($_+1)} = $_<2 ? ($dims[$_]+1)/2 : 1 for 0..$ndims-1; + $h->{"CTYPE".($_+1)} = $axislabels->[$_] for 0..$ndims-1; + $h->{"CUNIT".($_+1)} = $axisunits->[$_] for 0..$ndims-1; + $h->{"CDELT".($_+1)} = ($axisranges->[$_][1]-$axisranges->[$_][0])/$dims[$_] + for 0..$ndims-1; + $h->{HISTORY}='PDL conversion from raster image', + $out; } =head2 clean_lines @@ -842,47 +880,26 @@ sub t_raster2float { =head2 t_raster2fits -=for usage - - $t = t_raster2fits(); - =for ref -(Cartography) Convert a raster (3,x,y) to FITS plate carree (x,y,3) - -Adds suitable C. Assumes degrees. Used by L. +Deprecated as it's not actually a transformation, which operates +on coordinates. Use L. =cut sub t_raster2fits { my ($me) = _new(@_, 'Raster to FITS plate carree conversion'); - $me->{odim} = 3; $me->{params}->{itype} = ['RGB','X','Y']; $me->{params}->{iunit} = ['RGB','pixels','pixels']; $me->{params}->{otype} = ['X','Y','RGB']; $me->{params}->{ounit} = ['pixels','pixels','RGB']; $me->{func} = sub { my($d,$o) = @_; - my $out = $d->mv(0,2); - my $h = $out->fhdr; - $h->{SIMPLE} = 'T'; - $h->{NAXIS} = $d->ndims; - local $_; - $h->{"NAXIS".($_+1)} = $out->dim($_) for 0..$out->ndims-1; - $h->{"CRVAL".($_+1)} = 0 for 0..$out->ndims-1; - $h->{"CRPIX".($_+1)} = $_<2 ? ($out->dim($_)+1)/2 : 1 for 0..$out->ndims-1; - my ($lon, $lat) = $out->dims; - $h->{CTYPE1}='Longitude'; $h->{CUNIT1}='degrees'; $h->{CDELT1}=360/$lon; - $h->{CTYPE2}='Latitude'; $h->{CUNIT2}='degrees'; $h->{CDELT2}=180/$lat; - $h->{CTYPE3}='RGB'; $h->{CUNIT3}='index'; $h->{CDELT3}=1.0; - $h->{COMMENT}='Plate Carree Projection'; - $h->{HISTORY}='PDL conversion from raster image', - $out->hdrcpy(1); - $out; + raster2fits($d, @PLATE_CARREE); }; $me->{inv} = sub { my($d,$o) = @_; - $d->mv(2,0); + $d->ndims > 2 ? $d->mv(2,0) : $d; }; $me; } diff --git a/Libtmp/Transform/Proj4/t/proj_transform.t b/Libtmp/Transform/Proj4/t/proj_transform.t index 73e8e485b..e0d198c90 100644 --- a/Libtmp/Transform/Proj4/t/proj_transform.t +++ b/Libtmp/Transform/Proj4/t/proj_transform.t @@ -53,7 +53,7 @@ SKIP: { EOF my $shape = earth_shape(); - $got = t_raster2float()->apply($shape->mv(2,0)); + $got = t_raster2float()->apply($shape->dummy(0,3)); my $lonlatradius = $got->slice('0:2'); # r g b all same $lonlatradius->slice('(2)') *= float((6377.09863 - 6370.69873) / 6371); $lonlatradius->slice('(2)') += float(6370.69873 / 6371); diff --git a/Libtmp/Transform/t/transform.t b/Libtmp/Transform/t/transform.t index 32ea1a4de..6988d0a25 100644 --- a/Libtmp/Transform/t/transform.t +++ b/Libtmp/Transform/t/transform.t @@ -208,8 +208,8 @@ EOF { use PDL::Transform::Cartography; -my $pa = t_raster2fits()->apply(sequence(byte, 3, 10, 10)); -eval { $pa->match([100,100,3]) }; +my $pa = raster2fits(sequence(byte, 10, 10), @PDL::Transform::Cartography::PLATE_CARREE); +eval { $pa->match([100,100]) }; is $@, '', 't_fits invertible'; is earth_coast()->nbad, 0, 'earth_coast no BAD';