[postgis-users] Ellipses and WGS84

[Bruce Rindahl]

Your parameter r is in map units.  For your example r = 1 or a 1 degree 
circle.  The center section should be :
buffer(GeomFromText('Point(0.0 0.0)', 4326), 1)
Now you have a circle (polygon actually) with a radius of 1 unit 
(degrees) centered at point 0,0.
Next you need to scale the circle to create an ellipse.  You do this 
with the scale factors a and b.  These parameters are dimensionless 
(ratios).
The query becomes scale(buffer(GeomFromText('Point(0.0 0.0)', 4326), 
1),a,b)  where a is the meters per degree of longitude and b is the 
meters per degree of latitude at your final point cx,cy.  This gives you 
an ellipse at point(0 0) with the correct shape. 
Next is a rotation function to create an arbitrary ellipse.  You do not 
need this so eliminate it.
Finally translate the ellipse to your desired location - cx cy.
translate(scale(buffer(GeomFromText('Point(0.0 0.0)', 4326), 1),a,b),cx,cy)

Read it from the inside out. Make a unit circle (buffer) then scale it 
differently in x and y (scale (...),a,b) then move it to the location 
you want (translate(...),cx,cy)

Bruce

 David W Talmage - CONTRACTOR wrote:
> 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. :-(
>
>
>
> _______________________________________________
> postgis-users mailing list
> postgis-users at postgis.refractions.net
> http://postgis.refractions.net/mailman/listinfo/postgis-users
>
>