-------- Original Message --------
Subject: Re: [Perldl] Proj4 and Transform::map
Date: Fri, 21 Sep 2007 14:30:06 -0400
From: Greg Ederer <[EMAIL PROTECTED]>
To: Judd Taylor <[EMAIL PROTECTED]>
References: <[EMAIL PROTECTED]>
<[EMAIL PROTECTED]>
Hi Judd:
I am just inverse projecting from Sinusoidal to Geographic. I couldn't
make it work for other projections, so I thought I'd first try and get
it to work for Geo. I did use the "inverse" function; the line in the
original email message got split by the email formatter, so maybe you
didn't see it.
So far, I've tried:
$tf = t_proj( proj_params => '+proj=sinu +inv +over +a=6371007.181'
)->inverse();
$result = $pdl->map($tf, $to_FITS_hdr, { method=>'s' });
which produces the same output as the unexpected.jpg that was generated by
$tf = t_proj( proj_params => '+proj=sinu +over +a=6371007.181' )->inverse();
$result = $pdl->map($tf, $to_FITS_hdr, { method=>'s' });
Just for kicks, I also tried the forward version using the Proj4 inverse
transform
$tf = t_proj( proj_params => '+proj=sinu +inv +over +a=6371007.181' );
$result = $pdl->map($tf, $to_FITS_hdr, { method=>'s' });
but this complains about:
E_fwd_trans_inplace(): Projection conversion failed at
(-7783653.637667, 5559752.598333): latitude or longitude exceeded limits
which is expected, and finally, the forward projection using the forward
transform
$tf = t_proj( proj_params => '+proj=sinu +over +a=6371007.181');
$result = $pdl->map($tf, $to_FITS_hdr, { method=>'s' });
which (rightly) makes the same complaint.
I also tried the Transform::Cartography module's
"t_sinusoidal()->inverse()" instead (after appropriately scaling the
input to radians), and got the same unexpected image. This makes it look
like either a FITS problem, or a problem with the "map" function, or I'm
still missing something.
--greg
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
_______________________________________________
Perldl mailing list
[email protected]
http://mailman.jach.hawaii.edu/mailman/listinfo/perldl
