[postgis-users] Ellipses and WGS84

[David W Talmage - CONTRACTOR]

Bruce Rindahl wrote:
> The example creates an ellipse directly in a projected coordinate 
> space.  What you need is the difference between latitude and longitude 
> of 1 meter at the location you want.  Your code should be:
> translate(rotate(scale(buffer(GeomFromText('Point(0.0 0.0)', 4326), 
> r),a,b),alpha),cx,cy)
> 
> Where r is the radius in decimal degrees of 1 meter at the equator, b =1 
> (the 1 meter distance does not change in the N/S direction), alpha is 0 
> since you will not rotate the ellipse, and cx and cy are you final 
> lat/long locations.  The a parameter is tricky.  You need to know how 
> much larger 1 degree of longitude is than 1 degree of latitude is at 
> cx,cy.  You might try:
> ...

Thank you for helping me!

I should have worded the problem differently.  I need to make an ellipse
with a and b greater than 1.  Following the code, I start with a 1m
circle, then scale it up to a-meters by b-meters.

So the r, a, and b are all in degrees?  Then I have to figure out how 
many meters per degree at latitude cy and divide r, a, and b by that 
number, right?

If that's right, then something strange is going on.  I'm creating the 
ellipse that way and my SELECT asks for all things whose bounding box 
intersects that of the ellipse.  PostGIS is returning far fewer things 
that I expect.  Here's my SELECT:


SELECT rf.name FROM raster_entry re, raster_file rf WHERE 
re.raster_data_set_id = rf.raster_data_set_id AND rf.type = 'main' AND 
(re.ground_geom && translate(rotate(scale(buffer(GeomFromText('Point(0.0 
0.0)', 4326), ?), ?, ?), ?), ?, ?));


I'm doing this in Java with a PreparedStatement.  The '?' correspond to 
r, a, b, alpha, cx, and cy, respectively.


> David W Talmage - CONTRACTOR wrote:
>> Would someone please explain the units in this snippet modeled after 
>> the exampled offered in
>> http://postgis.refractions.net/pipermail/postgis-users/2007-March/014971.html 
>> ?
>>
>> translate(rotate(scale(buffer(GeomFromText('Point(0.0 0.0)', 4326), 
>> r),a,b),alpha),cx,cy)
>>
>> What are the units of r, a, and b?  I think they are radians.
>>
>> My intention is to create a 1-meter radius circle at the Equator, 
>> morph  it into an ellipse and then move the ellipse to an area of 
>> interest on the globe to select the things that are in that area.  All 
>> of my things are stored with coordinates in WGS84.
>>
>> I know the values of a and b in meters.  I have a function that 
>> computes the number of meters per degree on an ellipsoid at any 
>> latitude, so I can convert a and b to radians.
>>
>> I'm pretty sure that cx and cy are the longitude and latitude, 
>> respectively, in degrees of the center of the ellipse.
>>
>> alpha is in radians.
>>
>> I've looked all over for an answer to this question and either I can't 
>> find it or don't understand what I've found. :-(