Hey, all,
Here's hoping, Greg, that your inverse-flag issue raised by Judd is
the right
thing. I'm a bit puzzled that the apply() was working OK but not the
map().
They do work in opposite directions: map() uses the inverse of the
transform
you give it (to map the integer-valued positions in output space back
to real
valued positions in input space). So if there were a problem with
inverses in
the proj4 layer, it could yield the symptoms you describe.
You might also try mapping with the inverse of the transform you
describe:
just add a bang before the transform argument to map(), like so:
$a->map(!$transform);
If the problem is you not specifying that the transformation is an
inverse, then
you should get more reasonable output. If the problem is that
inverses aren't
being handled correctly, then you should get the same (or at least
similar)
output either way.
Cheers,
Craig
On Sep 20, 2007, at 5:05 PM, Judd Taylor wrote:
Greg,
From looking at your proj parameters string, I don't see anything
denoting it's an inverse projection. Please double check the proj docs
(google for proj4), to make sure you're aren't missing a parameter
like
+inv or -I or something like that.
Also, since you're starting from a projected dataset, you may need
to do
an inverse to get to "geographic coords" before projecting it into the
resulting image. Maybe you were doing this step before, but I don't
see
a proj command for it, nor an image below.
As a last note, maybe Craig DeForest will chime in on this thread.
He is
responsible for the general PDL::Transform code, whereas I just
added an
interface between that and a more general interface to the Proj4
library. I'll still help as best as I can, however.
-Judd
On Thu, 2007-09-20 at 15:11 -0400, Greg Ederer wrote:
Hi All:
I am starting to work with the PDL::Transform and Proj4 modules,
and am
running into some things I don't understand.
I have a file containing data in the Sinusoidal projection. The
corner
points are at
upper left: [-7783653.637667, 5559752.598333] meters
lower right: [-6671703.118, 4447802.078667] meters
pixel size: 926.62543305 meters
image size: 1200x1200
A jpeg image of the original, sinusoidal projection data is attached,
and named "small_original.jpg" (you can see parts of the Great Lakes,
especially the lower part of Lake Michigan and western part of Lake
Superior in black clearly in the image):
I "apply" the inverse sinusoidal transform to the corner points,
and get
upper left: [-108.900667860276, 49.9999999955069] degrees
lower right: [-78.324437348786, 39.9999999964109] degrees
pixel size: 0.00833366871395212 degrees
image size: 3669x1200
These corner points and sizes seem reasonable to me, and match what I
get when I use an alternative C program to do the reprojection
(see below).
However, when I "map" the data points using the inverse sinusoidal
transform, it looks like the image gets warped the wrong way -- east
instead of west. The jpeg of the image when I use the Proj4
transform is
attached and named "small_unexpected.jpg"
The FITS header that I assign to the input data is:
CDELT1: 926.625433055
CDELT2: 926.625433055
COMMENT: SNSOID
CRPIX1: 601
CRPIX2: 601
CRVAL1: -7227678.3778335
CRVAL2: 5003777.3385
CTYPE1: Longitude
CTYPE2: Latitude
CUNIT1: meters
CUNIT2: meters
HISTORY: FITS header created by MapGeometry.pm
NAXIS: 2
NAXIS1: 1200
NAXIS2: 1200
SIMPLE: T
The template (to) FITS header is:
CDELT1: 0.00833366871395212
CDELT2: 0.00833333333257998
COMMENT: GEO
CRPIX1: 1835.5
CRPIX2: 601
CRVAL1: -93.6125526045312
CRVAL2: 44.9999999959589
CTYPE1: Longitude
CTYPE2: Latitude
CUNIT1: degrees
CUNIT2: degrees
HISTORY: FITS header created by MapGeometry.pm
NAXIS: 2
NAXIS1: 3669
NAXIS2: 1200
SIMPLE: T
the function call looks like:
$tf = t_proj( proj_params => '+proj=sinu +over
+a=6371007.181'
)->inverse();
$result = $pdl->map($tf, $to_FITS_hdr, { method=>'s' });
For reference, I attached a jpeg (named named
"small_expected.jpg".) of
what I expected to see (as generated from an old C program I have
that
uses the GCTP library).
There are no other transforms applied to the data, and no other data
manipulation, except to read the data in and write them out. Am I
generating FITS headers incorrectly, or am I using the map function
incorrectly, or is there a problem with the map function? I am
using the
most recent CPAN version of PDL.
Thanks for your help,
--greg
NOTE: The attached jpegs were scaled smaller from the original jpegs
using a separate jpeg manipulation program. I wanted to be sure I was
only using one transform when running my test. The smaller sizes are
easier to see in their entirety and take up less space.
Also, MapGeometry.pm (mentioned in the FITS comments) is the
module I am
writing.
_______________________________________________
Perldl mailing list
[email protected]
http://mailman.jach.hawaii.edu/mailman/listinfo/perldl
--
____________________________
Judd Taylor
Software Engineer
Orbital Systems, Ltd.
3807 Carbon Rd.
Irving, TX 75038-3415
[EMAIL PROTECTED]
(972) 915-3669 x127
_______________________________________________
Perldl mailing list
[email protected]
http://mailman.jach.hawaii.edu/mailman/listinfo/perldl
_______________________________________________
Perldl mailing list
[email protected]
http://mailman.jach.hawaii.edu/mailman/listinfo/perldl
