"Fossies" - the Fresh Open Source Software Archive  

Source code changes of the file "Libtmp/Transform/Cartography/Cartography.pm" between
PDL-2.077.tar.gz and PDL-2.078.tar.gz

About: PDL (Perl Data Language) aims to turn perl into an efficient numerical language for scientific computing (similar to IDL and MatLab).

Cartography.pm  (PDL-2.077):Cartography.pm  (PDL-2.078)
skipping to change at line 34 skipping to change at line 34
Sun, the Celestial sphere, various moons and planets, distant stars, Sun, the Celestial sphere, various moons and planets, distant stars,
etc. They also are useful for interpreting scientific images, which etc. They also are useful for interpreting scientific images, which
are themselves generally projections of a sphere onto a flat focal are themselves generally projections of a sphere onto a flat focal
plane (e.g. the L</t_gnomonic> projection). plane (e.g. the L</t_gnomonic> projection).
Unless otherwise noted, all the transformations in this file convert Unless otherwise noted, all the transformations in this file convert
from (theta,phi) coordinates on the unit sphere (e.g. (lon,lat) on a from (theta,phi) coordinates on the unit sphere (e.g. (lon,lat) on a
planet or (RA,dec) on the celestial sphere) into some sort of planet or (RA,dec) on the celestial sphere) into some sort of
projected coordinates, and have inverse transformations that convert projected coordinates, and have inverse transformations that convert
back to (theta,phi). This is equivalent to working from the back to (theta,phi). This is equivalent to working from the
equidistant cylindrical (or L<"plate caree"|/t_caree>) projection, if equidistant cylindrical (or L<"plate carree"|/t_carree>) projection, if
you are a cartography wonk. you are a cartography wonk.
The projected coordinates are generally in units of body radii The projected coordinates are generally in units of body radii
(radians), so that multiplying the output by the scale of the map (radians), so that multiplying the output by the scale of the map
yields physical units that are correct wherever the scale is correct yields physical units that are correct wherever the scale is correct
for that projection. For example, areas should be correct everywhere for that projection. For example, areas should be correct everywhere
in the authalic projections; and linear scales are correct along in the authalic projections; and linear scales are correct along
meridians in the equidistant projections and along the standard meridians in the equidistant projections and along the standard
parallels in all the projections. parallels in all the projections.
skipping to change at line 258 skipping to change at line 258
use Exporter (); use Exporter ();
use PDL::LiteF; use PDL::LiteF;
use PDL::Math; use PDL::Math;
use PDL::Transform; use PDL::Transform;
use PDL::MatrixOps; use PDL::MatrixOps;
use Carp; use Carp;
our @ISA = ( 'Exporter','PDL::Transform' ); our @ISA = ( 'Exporter','PDL::Transform' );
our $VERSION = "0.6"; our $VERSION = "0.6";
$VERSION = eval $VERSION; $VERSION = eval $VERSION;
our @EXPORT_OK = qw(graticule earth_image earth_coast clean_lines t_unit_sphere our @EXPORT_OK = qw(
t_orthographic t_rot_sphere t_caree t_mercator t_utm t_sin_lat t_sinusoidal t_co graticule earth_image earth_coast earth_shape clean_lines t_unit_sphere
nic t_albers t_lambert t_stereographic t_gnomonic t_az_eqd t_az_eqa t_vertical t t_orthographic t_rot_sphere t_caree t_carree t_mercator t_utm t_sin_lat
_perspective t_hammer t_aitoff); 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
t_raster2fits t_raster2float
);
our @EXPORT = @EXPORT_OK; our @EXPORT = @EXPORT_OK;
our %EXPORT_TAGS = (Func=>\@EXPORT_OK); our %EXPORT_TAGS = (Func=>\@EXPORT_OK);
############################## ##############################
# Steal _opt from PDL::Transform. # Steal _opt from PDL::Transform.
*PDL::Transform::Cartography::_opt = \&PDL::Transform::_opt; *PDL::Transform::Cartography::_opt = \&PDL::Transform::_opt;
use overload '""' => \&_strval; use overload '""' => \&_strval;
our $PI = $PDL::Transform::PI; our $PI = $PDL::Transform::PI;
our $DEG2RAD = $PDL::Transform::DEG2RAD; our $DEG2RAD = $PDL::Transform::DEG2RAD;
skipping to change at line 390 skipping to change at line 396
=for usage =for usage
$x = earth_coast() $x = earth_coast()
=for ref =for ref
(Cartography) PDL constructor - coastline map of Earth (Cartography) PDL constructor - coastline map of Earth
Returns a vector coastline map based on the 1987 CIA World Coastline Returns a vector coastline map based on the 1987 CIA World Coastline
database (see author information). The vector coastline data are in database (see author information). The vector coastline data are in
plate caree format so they can be converted to other projections via plate carree format so they can be converted to other projections via
the L<apply|PDL::Transform/apply> method and cartographic transforms, the L<apply|PDL::Transform/apply> method and cartographic transforms,
and are suitable for plotting with the and are suitable for plotting with the
L<lines|PDL::Graphics::PGPLOT::Window/lines> L<lines|PDL::Graphics::PGPLOT::Window/lines>
method in the PGPLOT method in the PGPLOT
output library: the first dimension is (X,Y,pen) with breaks having output library: the first dimension is (X,Y,pen) with breaks having
a pen value of 0 and hairlines having negative pen values. The second a pen value of 0 and hairlines having negative pen values. The second
dimension broadcasts over all the points in the data set. dimension broadcasts over all the points in the data set.
The vector map includes lines that pass through the antipodean The vector map includes lines that pass through the antipodean
meridian, so if you want to plot it without reprojecting, you should meridian, so if you want to plot it without reprojecting, you should
run it through L</clean_lines> first: run it through L</clean_lines> first:
$w = pgwin(); $w = pgwin();
$w->lines(earth_coast->clean_lines); # plot plate caree map of world $w->lines(earth_coast->clean_lines); # plot plate carree map of world
$w->lines(earth_coast->apply(t_gnomonic))# plot gnomonic map of world $w->lines(earth_coast->apply(t_gnomonic))# plot gnomonic map of world
C<earth_coast> is just a quick-and-dirty way of loading the file C<earth_coast> is just a quick-and-dirty way of loading the file
"earth_coast.vec.fits" that is part of the normal installation tree. "earth_coast.vec.fits" that is part of the normal installation tree.
=cut =cut
sub earth_coast { sub earth_coast {
my $fn = "PDL/Transform/Cartography/earth_coast.vec.fits"; my $fn = "PDL/Transform/Cartography/earth_coast.vec.fits";
local $_; local $_;
skipping to change at line 436 skipping to change at line 442
$rgb = earth_image() $rgb = earth_image()
=for ref =for ref
(Cartography) PDL constructor - RGB pixel map of Earth (Cartography) PDL constructor - RGB pixel map of Earth
Returns an RGB image of Earth based on data from the MODIS instrument Returns an RGB image of Earth based on data from the MODIS instrument
on the NASA EOS/Terra satellite. (You can get a full-resolution on the NASA EOS/Terra satellite. (You can get a full-resolution
image from L<http://earthobservatory.nasa.gov/Newsroom/BlueMarble/>). image from L<http://earthobservatory.nasa.gov/Newsroom/BlueMarble/>).
The image is a plate caree map, so you can convert it to other The image is a plate carree map, so you can convert it to other
projections via the L<map|PDL::Transform/map> method and cartographic projections via the L<map|PDL::Transform/map> method and cartographic
transforms. transforms.
This is just a quick-and-dirty way of loading the earth-image files that This is just a quick-and-dirty way of loading the earth-image files that
are distributed along with PDL. are distributed along with PDL.
=cut =cut
sub earth_image { sub earth_image {
my($nd) = shift; my($nd) = shift;
my $f; my $f;
require PDL::IO::Pic; require PDL::IO::Pic;
my $dir = "PDL/Transform/Cartography/earth_"; my $dir = "PDL/Transform/Cartography/earth_";
$f = ($nd =~ m/^n/i) ? "${dir}night.jpg" : "${dir}day.jpg"; $f = (($nd//'') =~ m/^n/i) ? "${dir}night.jpg" : "${dir}day.jpg";
local $_; local $_;
my $im; my $im;
my $found = 0; my $found = 0;
foreach(@INC) { foreach(@INC) {
my $file = "$_/$f"; my $file = "$_/$f";
if(-e $file) { if(-e $file) {
$found = 1; $found = 1;
$im = PDL::IO::Pic::rpic($file)->mv(0,-1); $im = PDL::IO::Pic::rpic($file);
} }
last if defined($im); last if defined($im);
} }
barf("earth_image: $f not found in \@INC\n") barf("earth_image: $f not found in \@INC\n")
unless defined($found); unless defined($found);
barf("earth_image: couldn't load $f; you may need to install netpbm.\n") barf("earth_image: couldn't load $f; you may need to install netpbm.\n")
unless defined($im); unless defined($im);
t_raster2fits()->apply($im);
}
=head2 earth_shape
=for usage
$fits_shape = earth_shape()
=for ref
my $h = $im->fhdr; (Cartography) PDL constructor - height map of Earth
$h->{SIMPLE} = 'T'; Returns a height map of Earth based on data from the General
$h->{NAXIS} = 3; Bathymetric Chart of the Oceans (GEBCO) produced by the British
$h->{NAXIS1}=2048; $h->{CRPIX1}=1024.5; $h->{CRVAL1}=0; Oceanographic Data Centre. (You can get a full-resolution image from
$h->{NAXIS2}=1024; $h->{CRPIX2}=512.5; $h->{CRVAL2}=0; L<http://visibleearth.nasa.gov/view.php?id=73934>). The image is a
$h->{NAXIS3}=3, $h->{CRPIX3}=1; $h->{CRVAL3}=0; plate carree map, so you can convert it to other projections via the
$h->{CTYPE1}='Longitude'; $h->{CUNIT1}='degrees'; $h->{CDELT1}=180/1024.0; L<map|PDL::Transform/map> method and cartographic transforms.
$h->{CTYPE2}='Latitude'; $h->{CUNIT2}='degrees'; $h->{CDELT2}=180/1024.0; The data is from 8-bit grayscale (so only 256 levels), but is returned
$h->{CTYPE3}='RGB'; $h->{CUNIT3}='index'; $h->{CDELT3}=1.0; in a similar format to L</earth_image>. The range represents a span of
$h->{COMMENT}='Plate Caree Projection'; 6400m, so Everest and the Marianas Trench are not accurately represented.
$h->{HISTORY}='PDL Distribution Image, derived from NASA/MODIS data',
To turn this into a C<float>, (C<lonlatradius,x,y>) with C<x>
and C<y> in radians, and the radius as a C<float> as a proportion of the
Earth's mean radius, use L</t_raster2float>.
The Earth is treated here as a perfect sphere with sea
level at radius 6,371km.
Value Hex value Float From centre in km Float as radius
Base 00 0.0 6370.69873km 0.99995
Sea level 0C 0.04705 6371km 1.0
Highest FF 1.0 6377.09863km 1.00096
Code:
$shape = earth_shape();
$floats = t_raster2float()->apply($shape->mv(2,0));
$lonlatradius = $floats->slice('0:2'); # r g b all same
$lonlatradius->slice('(2)') *= float((6377.09863 - 6370.69873) / 6371);
$lonlatradius->slice('(2)') += float(6370.69873 / 6371);
$im->hdrcpy(1); =cut
$im;
sub earth_shape {
my($nd) = shift;
require PDL::IO::Pic;
my $f = "PDL/Transform/Cartography/earth_height-2048x1024.jpg";
local $_;
my $im;
my $found = 0;
foreach (@INC) {
my $file = "$_/$f";
if(-e $file) {
$found = 1;
$im = PDL::IO::Pic::rpic($file);
}
last if defined($im);
}
barf("earth_shape: $f not found in \@INC\n")
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);
} }
=head2 clean_lines =head2 clean_lines
=for usage =for usage
$x = clean_lines(t_mercator->apply(scalar(earth_coast()))); $x = clean_lines(t_mercator->apply(scalar(earth_coast())));
$x = clean_lines($line_pen, [threshold]); $x = clean_lines($line_pen, [threshold]);
$x = $lines->clean_lines; $x = $lines->clean_lines;
skipping to change at line 711 skipping to change at line 764
$out->{iunit} = $me->{iunit}; $out->{iunit} = $me->{iunit};
$out->{otype} = $me->{otype}; $out->{otype} = $me->{otype};
$out->{ounit} = $me->{ounit}; $out->{ounit} = $me->{ounit};
$out->{odim} = 2; $out->{odim} = 2;
$out->{idim} = 2; $out->{idim} = 2;
return $out; return $out;
} }
return $me; return $me;
} }
=head2 t_raster2float
=for usage
$t = t_raster2float();
=for ref
(Cartography) Convert a raster (3,x,y) to C<float> (lonlatrgb,x,y)
Assumes C<bytes> input, and radians and C<float> output, with the first
2 coordinates suitable for use as plate carree.
=cut
sub t_raster2float {
my ($me) = _new(@_, 'Raster FITS plate carree to OpenGL-ready float conversion
');
$me->{odim} = 3;
$me->{params}->{itype} = ['RGB','X','Y'];
$me->{params}->{iunit} = ['index','pixels','pixels'];
$me->{odim} = 2;
$me->{params}->{otype} = ['LonLatRGB','X','Y'];
$me->{params}->{ounit} = ['Float','index','index'];
$me->{func} = sub {
my($d,$o) = @_;
my (undef, $x, $y, @otherdims) = $d->dims;
my $type = float;
my $out_xy = zeroes(byte, $x, $y, @otherdims);
my $out = zeroes($type, 5, $x, $y, @otherdims);
$out->slice($_->[0]) .= $_->[1]
for ['(0)', $out_xy->xlinvals(-$PI, $PI)],
['(1)', $out_xy->ylinvals(-$PI/2, $PI/2)],
['2:4', $d->convert($type) / 255];
$out;
};
$me->{inv} = sub {
my($d,$o) = @_;
my $type = byte;
my (undef, $x, $y, @otherdims) = $d->dims;
my $out = zeroes($type, 3, $x, $y, @otherdims);
$out .= $d->slice('2:4') * 255;
$out;
};
$me;
}
=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<hdr>. Assumes degrees. Used by L</earth_image>.
=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;
};
$me->{inv} = sub {
my($d,$o) = @_;
$d->mv(2,0);
};
$me;
}
###################################################################### ######################################################################
=head2 t_unit_sphere =head2 t_unit_sphere
=for usage =for usage
$t = t_unit_sphere(<options>); $t = t_unit_sphere(<options>);
=for ref =for ref
skipping to change at line 852 skipping to change at line 998
(Cartography) Generate oblique projections (Cartography) Generate oblique projections
You feed in the origin in (theta,phi) and a roll angle, and you get back You feed in the origin in (theta,phi) and a roll angle, and you get back
out (theta', phi') coordinates. This is useful for making oblique or out (theta', phi') coordinates. This is useful for making oblique or
transverse projections: just compose t_rot_sphere with your favorite transverse projections: just compose t_rot_sphere with your favorite
projection and you get an oblique one. projection and you get an oblique one.
Most of the projections automagically compose themselves with t_rot_sphere Most of the projections automagically compose themselves with t_rot_sphere
if you feed in an origin or roll angle. if you feed in an origin or roll angle.
t_rot_sphere converts the base plate caree projection (straight lon, straight t_rot_sphere converts the base plate carree projection (straight lon, straight
lat) to a Cassini projection. lat) to a Cassini projection.
OPTIONS OPTIONS
=over 3 =over 3
=item STANDARD POSITIONAL OPTIONS =item STANDARD POSITIONAL OPTIONS
=back =back
skipping to change at line 1036 skipping to change at line 1182
$out->slice("(2)") *= -1 if($o->{m} == 2); $out->slice("(2)") *= -1 if($o->{m} == 2);
$o->{t_int}->invert($out); $o->{t_int}->invert($out);
}; };
$me; $me;
} }
###################################################################### ######################################################################
=head2 t_caree =head2 t_carree
=for usage =for usage
$t = t_caree(<options>); $t = t_carree(<options>);
=for ref =for ref
(Cartography) Plate Caree projection (cylindrical; equidistant) (Cartography) Plate Carree projection (cylindrical; equidistant)
This is the simple Plate Caree projection -- also called a "lat/lon plot". This is the simple Plate Carree projection -- also called a "lat/lon plot".
The horizontal axis is theta; the vertical axis is phi. This is a no-op The horizontal axis is theta; the vertical axis is phi. This is a no-op
if the angular unit is radians; it is a simple scale otherwise. if the angular unit is radians; it is a simple scale otherwise.
OPTIONS OPTIONS
=over 3 =over 3
=item STANDARD POSITIONAL OPTIONS =item STANDARD POSITIONAL OPTIONS
=item s, std, standard, Standard (default 0) =item s, std, standard, Standard (default 0)
The standard parallel where the transformation is conformal. Conformality The standard parallel where the transformation is conformal. Conformality
is achieved by shrinking of the horizontal scale to match the is achieved by shrinking of the horizontal scale to match the
vertical scale (which is correct everywhere). vertical scale (which is correct everywhere).
=back =back
=cut =cut
@PDL::Transform::Cartography::Caree::ISA = ('PDL::Transform::Cartography'); @PDL::Transform::Cartography::Carree::ISA = ('PDL::Transform::Cartography');
sub t_caree { *t_caree = *t_carree; # compat alias for poor spellers
my($me) = _new(@_,'Plate Caree Projection'); sub t_carree {
my($me) = _new(@_,'Plate Carree Projection');
my $p = $me->{params}; my $p = $me->{params};
$me->{otype} = ['projected longitude','latitude']; $me->{otype} = ['projected longitude','latitude'];
$me->{ounit} = ['proj. body radii','body radii']; $me->{ounit} = ['proj. body radii','body radii'];
$p->{stretch} = cos($p->{std}); $p->{stretch} = cos($p->{std});
$me->{func} = sub { $me->{func} = sub {
my($d,$o) = @_; my($d,$o) = @_;
my($out) = $d->is_inplace ? $d : $d->copy; my($out) = $d->is_inplace ? $d : $d->copy;
skipping to change at line 1536 skipping to change at line 1683
=item STANDARD POSITIONAL OPTIONS =item STANDARD POSITIONAL OPTIONS
=item s, std, Standard (default 29.5, 45.5) =item s, std, Standard (default 29.5, 45.5)
The locations of the standard parallel(s) (where the cone intersects The locations of the standard parallel(s) (where the cone intersects
the surface of the sphere). If you specify only one then the other is the surface of the sphere). If you specify only one then the other is
taken to be the nearest pole. If you specify both of them to be one taken to be the nearest pole. If you specify both of them to be one
pole then you get an equidistant azimuthal map. If you specify both pole then you get an equidistant azimuthal map. If you specify both
of them to be opposite and equidistant from the equator you get a of them to be opposite and equidistant from the equator you get a
Plate Caree projection. Plate Carree projection.
=back =back
=cut =cut
sub t_conic { sub t_conic {
my($me) = _c_new(@_,"Simple Conic Projection",[29.5,45.5]); my($me) = _c_new(@_,"Simple Conic Projection",[29.5,45.5]);
my($p) = $me->{params}; my($p) = $me->{params};
if($p->{cylindrical}) { if($p->{cylindrical}) {
print STDERR "Simple conic: degenerate case; using Plate Caree\n" print STDERR "Simple conic: degenerate case; using Plate Carree\n"
if($PDL::verbose); if($PDL::verbose);
return t_caree($me->{options}); return t_carree($me->{options});
} }
$p->{n} = ((cos($p->{std}->slice("(0)")) - cos($p->{std}->slice("(1)"))) $p->{n} = ((cos($p->{std}->slice("(0)")) - cos($p->{std}->slice("(1)")))
/ /
($p->{std}->slice("(1)") - $p->{std}->slice("(0)"))); ($p->{std}->slice("(1)") - $p->{std}->slice("(0)")));
$p->{G} = cos($p->{std}->slice("(0)"))/$p->{n} + $p->{std}->slice("(0)"); $p->{G} = cos($p->{std}->slice("(0)"))/$p->{n} + $p->{std}->slice("(0)");
$me->{otype} = ['Conic X','Conic Y']; $me->{otype} = ['Conic X','Conic Y'];
$me->{ounit} = ['Proj. radians','Proj. radians']; $me->{ounit} = ['Proj. radians','Proj. radians'];
 End of changes. 23 change blocks. 
36 lines changed or deleted 181 lines changed or added

Home  |  About  |  Features  |  All  |  Newest  |  Dox  |  Diffs  |  RSS Feeds  |  Screenshots  |  Comments  |  Imprint  |  Privacy  |  HTTP(S)